Experiment design

Originally Contributed by: Arpit Bhatia, Chris Coey

This tutorial covers experiment design examples (D-optimal, A-optimal, and E-optimal) from section 7.5 of the book Convex Optimization by Boyd and Vandenberghe.

The tutorial uses the following packages

using JuMP
import SCS
import LinearAlgebra
import Random
Info

This tutorial uses sets from MathOptInterface. By default, JuMP exports the MOI symbol as an alias for the MathOptInterface.jl package. We recommend making this more explicit in your code by adding the following lines:

import MathOptInterface
const MOI = MathOptInterface

We set a seed so the random numbers are repeatable:

Random.seed!(1234)
MersenneTwister(1234)

The relaxed experiment design problem

The basic experiment design problem is as follows.

Given the menu of possible choices for experiments, $v_{1}, \ldots, v_{p}$, and the total number $m$ of experiments to be carried out, choose the numbers of each type of experiment, i.e., $m_{1}, \ldots, m_{p}$ to make the error covariance $E$ small (in some sense).

The variables $m_{1}, \ldots, m_{p}$ must, of course, be integers and sum to $m$ the given total number of experiments. This leads to the optimization problem:

\[\begin{aligned} \min\left(\mathrm{w.r.t.} \mathbf{S}_{+}^{n}\right) & E=\left(\sum_{j=1}^{p} m_{j} v_{j} v_{j}^{T}\right)^{-1} \\ \text{subject to} & m_{i} \geq 0 \\ & \sum\limits_{i=1}^p m_{i} = m \\ & m_{i} \in \mathbb{Z},\quad i=1,\ldots,p \end{aligned}\]

The basic experiment design problem can be a hard combinatorial problem when $m$, the total number of experiments, is comparable to $n$, since in this case the $m_{i}$ are all small integers.

In the case when $m$ is large compared to $n$, however, a good approximate solution can be found by ignoring, or relaxing, the constraint that the $m_{i}$ are integers.

Let $\lambda_{i}=m_{i} / m,$ which is the fraction of the total number of experiments for which $a_{j}=v_{i},$ or the relative frequency of experiment $i$. We can express the error covariance in terms of $\lambda_{i}$ as:

\[E=\frac{1}{m}\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1}\]

The vector $\lambda \in \mathbf{R}^{p}$ satisfies $\lambda \succeq 0, \mathbf{1}^{T} \lambda=1,$ and also, each $\lambda_{i}$ is an integer multiple of $1 / m$. By ignoring this last constraint, we arrive at the problem:

\[\begin{aligned} \min\left(\mathrm{w.r.t.} \mathbf{S}_{+}^{n}\right) & E=(1 / m)\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1} \\ \text{subject to:}\quad & \lambda \succeq 0 \\ & \mathbf{1}^{T} \lambda=1 \end{aligned}\]

Several scalarizations have been proposed for the experiment design problem, which is a vector optimization problem over the positive semidefinite cone.

q = 4 # dimension of estimate space
p = 8 # number of experimental vectors
nmax = 3 # upper bound on lambda
n = 12

V = randn(q, p)

eye = Matrix{Float64}(LinearAlgebra.I, q, q);

A-optimal design

In A-optimal experiment design, we minimize tr $E$, the trace of the covariance matrix. This objective is simply the mean of the norm of the error squared:

\[\mathbf{E}\|e\|_{2}^{2}=\mathbf{E} \operatorname{tr}\left(e e^{T}\right)=\operatorname{tr} E\]

The A-optimal experiment design problem in SDP form is

\[\begin{aligned} \min & \mathbf{1}^{T} u \\ \text{subject to} & \left[\begin{aligned}{\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}} & {e_{k}} \\ {e_{k}^{T}} & {u_{k}}\end{aligned}\right] \succeq 0, \quad k=1, \ldots, n \\ & \lambda \succeq 0 \\ & \mathbf{1}^{T} \lambda=1 \end{aligned}\]

aOpt = Model(SCS.Optimizer)
set_silent(aOpt)
@variable(aOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(aOpt, u[1:q], lower_bound = 0)
@constraint(aOpt, sum(np) <= n)
for i in 1:q
    matrix = [
        V*LinearAlgebra.diagm(0 => np ./ n)*V' eye[:, i]
        eye[i, :]' u[i]
    ]
    @constraint(aOpt, matrix >= 0, PSDCone())
end
@objective(aOpt, Min, sum(u))
optimize!(aOpt)
objective_value(aOpt)
5.041251236485684
value.(np)
8-element Vector{Float64}:
  1.7479403176790234
  1.1153179223786454
 -1.921535647576537e-6
  1.6619528795282619
  3.0000000606181523
  0.8414252145434773
  1.3825638730556045
  2.250800850048898

E-optimal design

In $E$ -optimal design, we minimize the norm of the error covariance matrix, i.e. the maximum eigenvalue of $E$.

Since the diameter (twice the longest semi-axis) of the confidence ellipsoid $\mathcal{E}$ is proportional to $\|E\|_{2}^{1 / 2}$, minimizing $\|E\|_{2}$ can be interpreted geometrically as minimizing the diameter of the confidence ellipsoid.

E-optimal design can also be interpreted as minimizing the maximum variance of $q^{T} e$, over all $q$ with $\|q\|_{2}=1$. The E-optimal experiment design problem in SDP form is:

\[\begin{aligned} \min & t \\ \text{subject to} & \sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T} \succeq t I \\ & \lambda \succeq 0 \\ & \mathbf{1}^{T} \lambda=1 \end{aligned}\]

eOpt = Model(SCS.Optimizer)
set_silent(eOpt)
@variable(eOpt, 0 <= np[1:p] <= nmax)
@variable(eOpt, t)
@constraint(
    eOpt,
    V * LinearAlgebra.diagm(0 => np ./ n) * V' - (t .* eye) >= 0,
    PSDCone(),
)
@constraint(eOpt, sum(np) <= n)
@objective(eOpt, Max, t)
optimize!(eOpt)
objective_value(eOpt)
0.44894117767003006
value.(np)
8-element Vector{Float64}:
  2.999999556115164
  0.6751761431273157
 -1.6787116162002895e-6
  1.045839200265888
  3.0000003496920167
  1.787000623804308
  0.3015143744944239
  2.1904663003978775

D-optimal design

The most widely used scalarization is called $D$ -optimal design, in which we minimize the determinant of the error covariance matrix $E$. This corresponds to designing the experiment to minimize the volume of the resulting confidence ellipsoid (for a fixed confidence level). Ignoring the constant factor $1 / m$ in $E$, and taking the logarithm of the objective, we can pose this problem as convex optimization problem:

\[\begin{aligned} \min & \log \operatorname{det}\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1} \\ \text{subject to} & \lambda \succeq 0 \\ & \mathbf{1}^{T} \lambda=1 \end{aligned}\]

dOpt = Model(SCS.Optimizer)
set_silent(dOpt)
@variable(dOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(dOpt, t)
@objective(dOpt, Max, t)
@constraint(dOpt, sum(np) <= n)
E = V * LinearAlgebra.diagm(0 => np ./ n) * V'
@constraint(
    dOpt,
    [t, 1, (E[i, j] for i in 1:q for j in 1:i)...] in MOI.LogDetConeTriangle(q)
)
optimize!(dOpt)
objective_value(dOpt)
0.19015532184794745
value.(np)
8-element Vector{Float64}:
 -5.076115905111846e-7
  2.567412611814496
  7.077829975296941e-9
  0.2627487622133296
  2.9435559860747222
  2.3925058575379854
  2.8369409784128434
  0.9968446438338621

Tip

This tutorial was generated using Literate.jl. View the source .jl file on GitHub.