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 variables 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.