Example: experiment design
This tutorial was generated using Literate.jl. Download the source as a .jl file.
This tutorial was originally contributed by Arpit Bhatia and Chris Coey.
This tutorial covers experiment design examples (D-optimal, A-optimal, and E-optimal) from section 7.5 of (Boyd and Vandenberghe, 2004).
Required packages
This tutorial uses the following packages:
using JuMP
import SCS
import LinearAlgebra
import MathOptInterface as MOI
import RandomWe set a seed so the random numbers are repeatable:
Random.seed!(1234)Random.TaskLocalRNG()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, that is, $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
n_max = 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 = n_max)
@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)
@assert is_solved_and_feasible(aOpt)
objective_value(aOpt)5.103223362935127value.(np)8-element Vector{Float64}:
2.9495211350240993
1.7790921964770618
8.104871341338243e-6
4.706163785716405e-6
2.108266984184275
1.6983463798099494
1.2635010842547665
2.201229233017492E-optimal design
In $E$ -optimal design, we minimize the norm of the error covariance matrix, that is, 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] <= n_max)
@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)
@assert is_solved_and_feasible(eOpt)
objective_value(eOpt)0.43538638563660814value.(np)8-element Vector{Float64}:
2.9999986784944226
2.75281016389529
-3.6095749146944787e-7
-1.487244304015834e-6
2.181846279181994
2.325263546495149
0.21896764657670037
1.5211163050178793D-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 = n_max)
@variable(dOpt, t)
@objective(dOpt, Max, t)
@constraint(dOpt, sum(np) <= n)
E = V * LinearAlgebra.diagm(0 => np ./ n) * V'
@constraint(dOpt, [t; 1; triangle_vec(E)] in MOI.LogDetConeTriangle(q))
optimize!(dOpt)
@assert is_solved_and_feasible(dOpt)
objective_value(dOpt)0.30845393467249815value.(np)8-element Vector{Float64}:
0.42758242662215623
2.9100163168940334
2.9168724590080685e-6
2.629197980393792e-6
2.9157806299837166
2.67337566459234
2.735395012219622
0.3378388086258122