# Optimal Experimental Design¶

Optimal experimental design is a theory at the interface of statistics and optimization, which studies how to allocate some statistical trials within a set of available design points. The goal is to allow for the best possible estimation of an unknown parameter . In what follows, we assume the standard linear model with multiresponse experiments: a trial in the design point gives a multidimensional observation that can be written as , where is of dimension , is a matrix, and the error vectors are i.i.d. with a unit variance.

Several optimization criteria exist, leading to different SDP, SOCP and LP formulations. As such, optimal experimental design problens are natural examples for problems in conic optimization. For a review of the different formulations and more references, see [1].

The code below initializes the data used in all the examples of this page. It should be run prior to any of the codes presented in this page.

```
>>> import cvxopt as cvx
>>> import picos
>>> #---------------------------------#
>>> # First generate some data : #
>>> # _ a list of 8 matrices A #
>>> # _ a vector c #
>>> #---------------------------------#
>>> A = [cvx.matrix([[1,0,0,0,0],
... [0,3,0,0,0],
... [0,0,1,0,0]]),
... cvx.matrix([[0,0,2,0,0],
... [0,1,0,0,0],
... [0,0,0,1,0]]),
... cvx.matrix([[0,0,0,2,0],
... [4,0,0,0,0],
... [0,0,1,0,0]]),
... cvx.matrix([[1,0,0,0,0],
... [0,0,2,0,0],
... [0,0,0,0,4]]),
... cvx.matrix([[1,0,2,0,0],
... [0,3,0,1,2],
... [0,0,1,2,0]]),
... cvx.matrix([[0,1,1,1,0],
... [0,3,0,1,0],
... [0,0,2,2,0]]),
... cvx.matrix([[1,2,0,0,0],
... [0,3,3,0,5],
... [1,0,0,2,0]]),
... cvx.matrix([[1,0,3,0,1],
... [0,3,2,0,0],
... [1,0,0,2,0]])
... ]
>>> c = cvx.matrix([1,2,3,4,5])
```

## Multi-response c-optimal design (SOCP)¶

We compute the c-optimal design (`c=[1,2,3,4,5]`

)
for the observation matrices `A[i].T`

from the variable `A`

defined above.
The results below suggest that we should allocate 12.8% of the
experimental effort on design point #5, and 87.2% on the design point #7.

Primal problem

The SOCP for multiresponse c-optimal design is:

```
>>> # create the problem, variables and params
>>> c_primal_SOCP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> cc = picos.Constant('c', c)
>>> z = [picos.RealVariable('z[{0}]'.format(i), AA[i].size[1])
... for i in range(s)]
>>> mu = picos.RealVariable('mu', s)
```

```
>>> # define the constraints and objective function
>>> cones = c_primal_SOCP.add_list_of_constraints([abs(z[i]) <= mu[i] for i in range(s)])
>>> lin = c_primal_SOCP.add_constraint(picos.sum([AA[i] * z[i] for i in range(s)]) == cc)
>>> c_primal_SOCP.set_objective('min', (1|mu) )
>>> print(c_primal_SOCP)
Second Order Cone Program
minimize ∑(mu)
over
3×1 real variable z[i] ∀ i ∈ [0…7]
8×1 real variable mu
subject to
‖z[i]‖ ≤ mu[i] ∀ i ∈ [0…7]
∑(A[i]·z[i] : i ∈ [0…7]) = c
```

```
>>> #solve the problem and retrieve the optimal weights of the optimal design.
>>> solution = c_primal_SOCP.solve(solver='cvxopt')
>>> mu = mu.value
>>> w = mu / sum(mu) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[...]
[...]
[ 1.28e-01]
[...]
[ 8.72e-01]
[...]
```

The `[...]`

above indicate a numerical zero entry
(*i.e., which can be something like 2.84e-10*).
We use the ellipsis `...`

instead for clarity and compatibility with **doctest**.

Dual problem

This is only to check that we obtain the same solution with the dual problem, and to provide one additional example in this tutorial:

```
>>> # create the problem, variables and params
>>> c_dual_SOCP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> cc = picos.Constant('c',c)
>>> u = picos.RealVariable('u',c.size)
>>> # define the constraints and objective function
>>> cones = c_dual_SOCP.add_list_of_constraints(
... [abs(AA[i].T*u)<=1 for i in range(s)])
>>> c_dual_SOCP.set_objective('max', (cc|u) )
>>> print(c_dual_SOCP)#
Second Order Cone Program
maximize ⟨c, u⟩
over
5×1 real variable u
subject to
‖A[i]ᵀ·u‖ ≤ 1 ∀ i ∈ [0…7]
>>> #solve the problem and retrieve the weights of the optimal design
>>> solution = c_dual_SOCP.solve(solver='cvxopt')
>>> mu = [cons.dual[0] for cons in cones] #Lagrangian duals of the SOC constraints
>>> mu = cvx.matrix(mu)
>>> w=mu/sum(mu) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[...]
[...]
[ 1.28e-01]
[...]
[ 8.72e-01]
[...]
```

## Single-response c-optimal design (LP)¶

When the observation matrices are row vectors (single-response framework),
the SOCP above reduces to a simple LP, because the variables
are scalar.
We solve below the LP for the case where there are 11
available design points, corresponding to the columns of the matrices
`A[4]`

, `A[5]`

, `A[6]`

, and `A[7][:,:-1]`

defined in the preambule.

The optimal design allocates 3.37% to point #5 (2nd column of `A[5]`

),
27.9% to point #7 (1st column of `A[6]`

),
11.8% to point #8 (2nd column of `A[6]`

),
27.6% to point #9 (3rd column of `A[6]`

),
and 29.3% to point #11 (2nd column of `A[7]`

).

```
>>> # create the problem, variables and params
>>> c_primal_LP = picos.Problem()
>>> A1 = [cvx.sparse(a[:,i],tc='d') for i in range(3) for a in A[4:]] #12 column vectors
>>> A1 = A1[:-1] # remove the last design point (it is the same as the last-but-one)
>>> s = len(A1)
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A1)] # each AA[i].T is a 1 x 5 observation matrix
>>> cc = picos.Constant('c', c)
>>> z = [picos.RealVariable('z[{0}]'.format(i), 1) for i in range(s)]
>>> mu = picos.RealVariable('mu', s)
```

```
>>> #define the constraints and objective function
>>> abs_con = c_primal_LP.add_list_of_constraints([abs(z[i]) <= mu[i] for i in range(s)])
>>> lin_con = c_primal_LP.add_constraint(picos.sum([AA[i]*z[i] for i in range(s)]) == cc)
>>> c_primal_LP.set_objective('min', (1|mu))
```

Note that there are no cone constraints, because the constraints of the form are handled as two inequalities when is scalar, so the problem is a LP indeed:

```
>>> print(c_primal_LP)
Linear Program
minimize ∑(mu)
over
1×1 real variable z[i] ∀ i ∈ [0…10]
11×1 real variable mu
subject to
|z[i]| ≤ mu[i] ∀ i ∈ [0…10]
∑(A[i]·z[i] : i ∈ [0…10]) = c
```

```
>>> #solve the problem and retrieve the weights of the optimal design
>>> solution = c_primal_LP.solve(solver='cvxopt')
>>> mu = mu.value
>>> w = mu / sum(mu) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[...]
[...]
[ 3.37e-02]
[...]
[ 2.79e-01]
[ 1.18e-01]
[ 2.76e-01]
[...]
[ 2.93e-01]
```

## SDP formulation of c-optimal design¶

We give below the SDP for c-optimality, in primal and dual form. You can observe that we obtain the same results as with the SOCP presented earlier: 12.8% on design point #5, and 87.2% on design point #7.

Primal problem

The SDP formulation of the c-optimal design problem is:

```
>>> # create the problem, variables and params
>>> c_primal_SDP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> cc = picos.Constant('c', c)
>>> mu = picos.RealVariable('mu',s)
>>> # define the constraints and objective function
>>> lmi = c_primal_SDP.add_constraint(
... picos.sum([mu[i] * AA[i] * AA[i].T for i in range(s)]) >> cc*cc.T)
>>> lin_cons = c_primal_SDP.add_constraint(mu >= 0)
>>> c_primal_SDP.set_objective('min', (1|mu) )
>>> print(c_primal_SDP)
Semidefinite Program
minimize ∑(mu)
over
8×1 real variable mu
subject to
∑(mu[i]·A[i]·A[i]ᵀ : i ∈ [0…7]) ≽ c·cᵀ
mu ≥ 0
```

```
>>> #solve the problem and retrieve the weights of the optimal design
>>> solution = c_primal_SDP.solve(solver='cvxopt')
>>> w = mu.value
>>> w = w / sum(w) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[...]
[...]
[ 1.28e-01]
[...]
[ 8.72e-01]
[...]
```

Dual problem

This is only to check that we obtain the same solution with the dual problem, and to provide one additional example in this tutorial:

```
>>> #create the problem, variables and params
>>> c_dual_SDP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> cc = picos.Constant('c', c)
>>> m = c.size[0]
>>> X = picos.SymmetricVariable('X',(m,m))
```

```
>>> #define the constraints and objective function
>>> lin_cons = c_dual_SDP.add_list_of_constraints(
... [(AA[i]*AA[i].T | X ) <= 1 for i in range(s)])
>>> psd = c_dual_SDP.add_constraint(X>>0)
>>> c_dual_SDP.set_objective('max', cc.T*X*cc)
```

```
>>> print(c_dual_SDP)
Semidefinite Program
maximize cᵀ·X·c
over
5×5 symmetric variable X
subject to
⟨A[i]·A[i]ᵀ, X⟩ ≤ 1 ∀ i ∈ [0…7]
X ≽ 0
```

```
>>> # solve the problem and retrieve the weights of the optimal design
>>> solution = c_dual_SDP.solve(solver='cvxopt')
>>> mu = [cons.dual for cons in lin_cons] #Lagrangian duals of the linear constraints
>>> mu = cvx.matrix(mu)
>>> w = mu / sum(mu) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[...]
[...]
[ 1.28e-01]
[...]
[ 8.72e-01]
[...]
```

And the optimal positive semidefinite matrix X is:

```
>>> print(X)
[ 5.92e-03 8.98e-03 2.82e-03 -3.48e-02 -1.43e-02]
[ 8.98e-03 1.36e-02 4.27e-03 -5.28e-02 -2.17e-02]
[ 2.82e-03 4.27e-03 1.34e-03 -1.66e-02 -6.79e-03]
[-3.48e-02 -5.28e-02 -1.66e-02 2.05e-01 8.39e-02]
[-1.43e-02 -2.17e-02 -6.79e-03 8.39e-02 3.44e-02]
```

## A-optimality (SOCP)¶

We compute the A-optimal design
for the observation matrices `A[i].T`

defined in the preambule.
The optimal design allocates
24.9% on design point #3,
14.2% on point #4,
8.51% on point #5,
12.1% on point #6,
13.2% on point #7,
and 27.0% on point #8.

Primal problem

The SOCP for the A-optimal design problem is:

```
>>> # create the problem, variables and params
>>> A_primal_SOCP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> Z = [picos.RealVariable('Z[{0}]'.format(i), AA[i].T.size) for i in range(s)]
>>> mu = picos.RealVariable('mu', s)
```

```
>>> #define the constraints and objective function
>>> cone_cons = A_primal_SOCP.add_list_of_constraints(
... [abs(Z[i]) <= mu[i] for i in range(s)])
>>> lin_cons = A_primal_SOCP.add_constraint(
... picos.sum([AA[i] * Z[i] for i in range(s)]) == 'I')
>>> A_primal_SOCP.set_objective('min', (1|mu) )
>>> print(A_primal_SOCP)
Second Order Cone Program
minimize ∑(mu)
over
3×5 real variable Z[i] ∀ i ∈ [0…7]
8×1 real variable mu
subject to
‖Z[i]‖ ≤ mu[i] ∀ i ∈ [0…7]
∑(A[i]·Z[i] : i ∈ [0…7]) = I
```

```
>>> # solve the problem and retrieve the weights of the optimal design
>>> solution = A_primal_SOCP.solve(solver='cvxopt')
>>> w = mu.value
>>> w = w / sum(w) #normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[ 2.49e-01]
[ 1.42e-01]
[ 8.51e-02]
[ 1.21e-01]
[ 1.32e-01]
[ 2.70e-01]
```

Dual problem

This is only to check that we obtain the same solution with the dual problem, and to provide one additional example in this tutorial:

```
>>> #create the problem, variables and params
>>> D_SOCPual_A=picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> m = AA[0].size[0]
>>> U = picos.RealVariable('U',(m,m))
>>> #define the constraints and objective function
>>> cone_cons = D_SOCPual_A.add_list_of_constraints(
... [abs(AA[i].T*U) <= 1 for i in range(s)])
>>> D_SOCPual_A.set_objective('max', 'I'|U)
>>> print(D_SOCPual_A)
Second Order Cone Program
maximize tr(U)
over
5×5 real variable U
subject to
‖A[i]ᵀ·U‖ ≤ 1 ∀ i ∈ [0…7]
```

```
>>> # solve the problem and retrieve the weights of the optimal design
>>> solution = D_SOCPual_A.solve(solver='cvxopt')
>>> mu = [cons.dual[0] for cons in cone_cons] # Lagrangian duals of the SOC constraints
>>> mu = cvx.matrix(mu)
>>> w = mu / sum(mu) # normalize mu to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[ 2.49e-01]
[ 1.42e-01]
[ 8.51e-02]
[ 1.21e-01]
[ 1.32e-01]
[ 2.70e-01]
```

## A-optimality with multiple constraints (SOCP)¶

A-optimal designs can also be computed by SOCP when the vector of weights is subject to several linear constraints. To give an example, we compute the A-optimal design for the observation matrices given in the preambule, when the weights must satisfy: and . This problem has the following SOCP formulation:

The optimal solution allocates 29.7% and 20.3% to the design points #3 and #4, and respectively 6.54%, 11.9%, 9.02% and 22.5% to the design points #5 to #8:

```
>>> # create the problem, variables and params
>>> A_multiconstraints = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> mu = picos.RealVariable('mu',s)
>>> w = picos.RealVariable('w',s)
>>> Z = [picos.RealVariable('Z[{0}]'.format(i), AA[i].T.size)
... for i in range(s)]
>>> # define the constraints and objective function
>>> lin_cons0 = A_multiconstraints.add_constraint(
... picos.sum([AA[i] * Z[i] for i in range(s)]) == 'I')
>>> lin_cons1 = A_multiconstraints.add_constraint( (1|w[:4]) <= 0.5)
>>> lin_cons2 = A_multiconstraints.add_constraint( (1|w[4:]) <= 0.5)
>>> cone_cons = A_multiconstraints.add_list_of_constraints(
... [ abs(Z[i]) **2 <= mu[i] * w[i] for i in range(s)])
>>> A_multiconstraints.set_objective('min', (1|mu) )
>>> print(A_multiconstraints)
Quadratically Constrained Program
minimize ∑(mu)
over
3×5 real variable Z[i] ∀ i ∈ [0…7]
8×1 real variable mu, w
subject to
∑(A[i]·Z[i] : i ∈ [0…7]) = I
∑(w[:4]) ≤ 0.5
∑(w[4:]) ≤ 0.5
‖Z[i]‖² ≤ mu[i]·w[i] ∀ i ∈ [0…7]
```

```
>>> # solve the problem and retrieve the weights of the optimal design
>>> solution = A_multiconstraints.solve(solver='cvxopt')
>>> w = w.value
>>> w = w / sum(w) # normalize w to get the optimal weights
```

The optimal design is:

```
>>> print(w)
[...]
[...]
[ 2.97e-01]
[ 2.03e-01]
[ 6.54e-02]
[ 1.19e-01]
[ 9.02e-02]
[ 2.25e-01]
```

## Exact A-optimal design (MISOCP)¶

In the exact version of A-optimality, a number of trials is given, and the goal is to find the optimal number of times that a trial on design point #i should be performed, with .

The SOCP formulation of A-optimality for constrained designs also accept integer constraints, which results in a MISOCP for exact A-optimality:

The exact optimal design is :

```
>>> # create the problem, variables and params
>>> A_exact = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> m = AA[0].size[0]
>>> N = picos.Constant('N', 20) # number of trials allowed
>>> I = picos.Constant('I', cvx.spmatrix([1]*m,range(m),range(m),(m,m))) #identity matrix
>>> Z = [picos.RealVariable('Z[{0}]'.format(i), AA[i].T.size) for i in range(s)]
>>> n = picos.IntegerVariable('n', s)
>>> t = picos.RealVariable('t', s)
```

```
>>> # define the constraints and objective function
>>> cone_cons = A_exact.add_list_of_constraints(
... [ abs(Z[i])**2 <= n[i] * t[i] for i in range(s)])
>>> lin_cons = A_exact.add_constraint(
... picos.sum([AA[i]*Z[i] for i in range(s)]) == I)
>>> wgt_cons = A_exact.add_constraint( (1|n) <= N )
>>> A_exact.set_objective('min',1|t)
>>> print(A_exact)
Mixed-Integer Quadratically Constrained Program
minimize ∑(t)
over
8×1 integer variable n
3×5 real variable Z[i] ∀ i ∈ [0…7]
8×1 real variable t
subject to
‖Z[i]‖² ≤ n[i]·t[i] ∀ i ∈ [0…7]
∑(A[i]·Z[i] : i ∈ [0…7]) = I
∑(n) ≤ N
```

```
>>> #solve the problem and display the optimal design
>>> solution = A_exact.solve()
>>> print(n)
[...]
[...]
[ 5.00e+00]
[ 3.00e+00]
[ 2.00e+00]
[ 2.00e+00]
[ 3.00e+00]
[ 5.00e+00]
```

Note

The above output is not validated as we lack an appropriate solver on the build server.

## Approximate and exact D-optimal design ((MI)SOCP)¶

The D-optimal design problem has a SOCP formulation involving a geometric mean in the objective function:

By introducing a new variable such that
, we can pass
this problem to PICOS with the function `geomean`

,
which reformulates the geometric mean inequality as a set of equivalent second order cone
constraints.
The example below allocates respectively 22.7%, 3.38%, 1.65%, 5.44%, 31.8% and 35.1%
to the design points #3 to #8.

```
>>> #create the problem, variables and params
>>> D_SOCP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> m = AA[0].size[0]
>>> mm = picos.Constant('m', m)
>>> L = picos.RealVariable('L', (m,m))
>>> V = [picos.RealVariable('V['+str(i)+']', AA[i].T.size) for i in range(s)]
>>> w = picos.RealVariable('w',s)
>>> # additional variable to handle the geometric mean in the objective function
>>> t = picos.RealVariable('t',1)
```

```
>>> # define the constraints and objective function
>>> lin_cons = D_SOCP.add_constraint(picos.sum([AA[i]*V[i] for i in range(s)]) == L)
>>> # L is lower triangular
>>> lowtri_cons = D_SOCP.add_list_of_constraints( [L[i,j] == 0
... for i in range(m)
... for j in range(i+1,m)])
>>> cone_cons = D_SOCP.add_list_of_constraints([abs(V[i]) <= (mm**0.5)*w[i]
... for i in range(s)])
>>> wgt_cons = D_SOCP.add_constraint(1|w <= 1)
>>> geomean_cons = D_SOCP.add_constraint(t <= picos.geomean(picos.maindiag(L)))
>>> D_SOCP.set_objective('max',t)
```

```
>>> #solve the problem and display the optimal design
>>> print(D_SOCP)
Optimization Problem
maximize t
over
1×1 real variable t
3×5 real variable V[i] ∀ i ∈ [0…7]
5×5 real variable L
8×1 real variable w
subject to
L = ∑(A[i]·V[i] : i ∈ [0…7])
L[i,j] = 0 ∀ (i,j) ∈ zip([0,0,…,2,3],[1,2,…,4,4])
‖V[i]‖ ≤ m^(1/2)·w[i] ∀ i ∈ [0…7]
∑(w) ≤ 1
geomean(maindiag(L)) ≥ t
```

```
>>> solution = D_SOCP.solve(solver='cvxopt')
>>> print(w)
[...]
[...]
[ 2.27e-01]
[ 3.38e-02]
[ 1.65e-02]
[ 5.44e-02]
[ 3.18e-01]
[ 3.51e-01]
```

As for the A-optimal problem, there is an alternative SOCP formulation of D-optimality [2], in which integer constraints may be added. This allows us to formulate the exact D-optimal problem as a MISOCP. For , we obtain the following N-exact D-optimal design: :

```
>>> # create the problem, variables and params
>>> D_exact = picos.Problem()
>>> L = picos.RealVariable('L',(m,m))
>>> V = [picos.RealVariable('V['+str(i)+']',AA[i].T.size) for i in range(s)]
>>> T = picos.RealVariable('T', (s,m))
>>> n = picos.IntegerVariable('n', s)
>>> N = picos.Constant('N', 20)
>>> # additional variable to handle the geomean inequality
>>> t = picos.RealVariable('t',1)
```

```
>>> # define the constraints and objective function
>>> lin_cons = D_exact.add_constraint(
... picos.sum([AA[i]*V[i] for i in range(s)]) == L)
>>> # L is lower triangular
>>> lowtri_cons = D_exact.add_list_of_constraints( [L[i,j] == 0
... for i in range(m)
... for j in range(i+1,m)])
>>> cone_cons = D_exact.add_list_of_constraints([ abs(V[i][:,k])**2 <= n[i]/N*T[i,k]
... for i in range(s) for k in range(m)])
>>> lin_cons2 = D_exact.add_list_of_constraints([(1|T[:,k]) <= 1
... for k in range(m)])
>>> wgt_cons = D_exact.add_constraint(1|n <= N)
>>> geomean_cons = D_exact.add_constraint(t <= picos.geomean( picos.maindiag(L)))
>>> D_exact.set_objective('max',t)
>>> print(D_exact)
Mixed-Integer Optimization Problem
maximize t
over
8×1 integer variable n
1×1 real variable t
3×5 real variable V[i] ∀ i ∈ [0…7]
5×5 real variable L
8×5 real variable T
subject to
L = ∑(A[i]·V[i] : i ∈ [0…7])
L[i,j] = 0 ∀ (i,j) ∈ zip([0,0,…,2,3],[1,2,…,4,4])
‖V[i][:,j]‖² ≤ n[i]/N·T[i,j] ∀ (i,j) ∈
zip([0,0,…,7,7],[0,1,…,3,4])
∑(T[:,i]) ≤ 1 ∀ i ∈ [0…4]
∑(n) ≤ N
geomean(maindiag(L)) ≥ t
```

```
>>> #solve the problem and display the optimal design
>>> solution = D_exact.solve()
>>> print(n)
[...]
[...]
[ 5.00e+00]
[ 1.00e+00]
[...]
[ 1.00e+00]
[ 6.00e+00]
[ 7.00e+00]
```

Note

The above output is not validated as we lack an appropriate solver on the build server.

## Former MAXDET formulation of the D-optimal design (SDP)¶

A so-called MAXDET Programming formulation of the D-optimal design
has been known since the late 90’s [3], and
can be reformulated as a SDP thanks to the `detrootn`

function.
The following code finds the same design as the SOCP approach presented above.

```
>>> # problem, variables and parameters
>>> D_MAXDET = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> m = AA[0].size[0]
>>> w = picos.RealVariable('w', s, lower=0)
>>> t = picos.RealVariable('t', 1)
>>> # constraint and objective
>>> wgt_cons = D_MAXDET.add_constraint(1|w <= 1)
>>> Mw = picos.sum([w[i] * AA[i] * AA[i].T for i in range(s)])
>>> detrootn_cons = D_MAXDET.add_constraint(t <= picos.DetRootN(Mw))
>>> D_MAXDET.set_objective('max', t)
```

```
>>> print(D_MAXDET)
Optimization Problem
maximize t
over
1×1 real variable t
8×1 real variable w (bounded below)
subject to
∑(w) ≤ 1
det(∑(w[i]·A[i]·A[i]ᵀ : i ∈ [0…7]))^(1/5) ≥ t
```

```
>>> #solve and display
>>> solution = D_MAXDET.solve(solver='cvxopt')
>>> print(w)
[ ...]
[ ...]
[ 2.27e-01]
[ 3.38e-02]
[ 1.65e-02]
[ 5.44e-02]
[ 3.18e-01]
[ 3.51e-01]
```

## General Phi_p optimal design (SDP)¶

The A- and D-optimal design problems presented above can be obtained as special cases of the general Kiefer optimal design problem, where is a real in :

These problems are easy to enter in PICOS, thanks to the `tracepow`

function,
that automatically replaces inequalities involving trace of matrix powers as a set of equivalent linear matrix
inequalities (SDP) (cf. [4] ). Below are two examples with and ,
allocating respectively (20.6%, 0.0%, 0.0%, 0.92%, 40.8%, 37.7%), and
(24.8%, 16.6%, 10.8%, 14.1%, 7.84%, 26.0%) of the trials to the design points 3 to 8.

```
>>> #problems, variables and parameters
>>> P0dot2_SDP = picos.Problem()
>>> Pminus3_SDP = picos.Problem()
>>> AA = [picos.Constant('A[{0}]'.format(i), Ai)
... for i, Ai in enumerate(A)] # each AA[i].T is a 3 x 5 observation matrix
>>> s = len(AA)
>>> m = AA[0].size[0]
>>> w = picos.RealVariable('w', s, lower=0)
>>> t = picos.RealVariable('t', 1)
```

```
>>> # constraint and objective
>>> wgt02_cons = P0dot2_SDP.add_constraint(1|w <= 1)
>>> wgtm3_cons = Pminus3_SDP.add_constraint(1|w <= 1)
```

```
>>> Mw = picos.sum([w[i]*AA[i]*AA[i].T for i in range(s)])
```

```
>>> tracep02_cons = P0dot2_SDP.add_constraint(t <= picos.PowerTrace(Mw, 0.2))
>>> P0dot2_SDP.set_objective('max', t)
```

```
>>> tracepm3_cons = Pminus3_SDP.add_constraint(t >= picos.PowerTrace(Mw, -3))
>>> Pminus3_SDP.set_objective('min', t)
```

```
>>> # p=0.2
>>> print(P0dot2_SDP)
Optimization Problem
maximize t
over
1×1 real variable t
8×1 real variable w (bounded below)
subject to
∑(w) ≤ 1
tr(∑(w[i]·A[i]·A[i]ᵀ : i ∈ [0…7])^(1/5)) ≥ t
```

```
>>> #solve and display
>>> solution = P0dot2_SDP.solve(solver='cvxopt')
>>> print(w)
[ ...]
[ ...]
[ 2.06e-01]
[ ...]
[ ...]
[ 9.20e-03]
[ 4.08e-01]
[ 3.77e-01]
```

```
>>> # p=-3
>>> print(Pminus3_SDP)
Optimization Problem
minimize t
over
1×1 real variable t
8×1 real variable w (bounded below)
subject to
∑(w) ≤ 1
tr(∑(w[i]·A[i]·A[i]ᵀ : i ∈ [0…7])^(-3)) ≤ t
>>> solution = Pminus3_SDP.solve(solver='cvxopt')
>>> print(w)
[ ...]
[ ...]
[ 2.48e-01]
[ 1.66e-01]
[ 1.08e-01]
[ 1.41e-01]
[ 7.83e-02]
[ 2.60e-01]
```

## References¶

“Computing Optimal Designs of multiresponse Experiments reduces to Second-Order Cone Programming”, G. Sagnol,

Journal of Statistical Planning and Inference, 141(5), p.1684-1708, 2011.“Computing exact D-optimal designs by mixed integer second order cone programming”, G. Sagnol and R. Harman, Submitted: arXiv:1307.4953.

“Determinant maximization with linear matrix inequality constraints”, L. Vandenberghe, S. Boyd and S.P. Wu,

SIAM journal on matrix analysis and applications, 19(2), 499-533, 1998.“On the semidefinite representations of real functions applied to symmetric matrices”, G. Sagnol,

Linear Algebra and its Applications, 439(10), p.2829-2843, 2013.