Tips and Tricks

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial was originally contributed by Arpit Bhatia.

This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP.

It uses the following packages:

using JuMP
import SCS
import LinearAlgebra

import MathOptInterface as MOI

Background theory

A subset $C$ of a vector space $V$ is a cone if $\forall x \in C$ and positive scalars $\lambda > 0$, the product $\lambda x \in C$.

A cone $C$ is a convex cone if $\lambda x + (1 - \lambda) y \in C$, for any $\lambda \in [0, 1]$, and any $x, y \in C$.

Conic programming problems are convex optimization problems in which a convex function is minimized over the intersection of an affine subspace and a convex cone. An example of a conic-form minimization problems, in the primal form is:

\[\begin{aligned} & \min_{x \in \mathbb{R}^n} & a_0^T x + b_0 \\ & \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m \end{aligned}\]

The corresponding dual problem is:

\[\begin{aligned} & \max_{y_1, \ldots, y_m} & -\sum_{i=1}^m b_i^T y_i + b_0 \\ & \;\;\text{s.t.} & a_0 - \sum_{i=1}^m A_i^T y_i & = 0 \\ & & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m \end{aligned}\]

where each $\mathcal{C}_i$ is a closed convex cone and $\mathcal{C}_i^*$ is its dual cone.

Second-Order Cone

The SecondOrderCone (or Lorentz Cone) of dimension $n$ is a cone of the form:

\[K_{soc} = \{ (t, x) \in \mathbb{R}^n : t \ge ||x||_2 \}\]

It is most commonly used to represent the L2-norm of the vector $x$:

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:3])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in SecondOrderCone())
@objective(model, Min, t)
optimize!(model)
value(t), value.(x)

(0.5773503148525415, [0.3333333333638126, 0.3333333333638124, 0.3333333333638124])

Rotated Second-Order Cone

A Second-Order Cone rotated by $\pi/4$ in the $(x_1,x_2)$ plane is called a RotatedSecondOrderCone. It is a cone of the form:

\[K_{rsoc} = \{ (t,u,x) \in \mathbb{R}^n : 2tu \ge ||x||_2^2, t,u \ge 0 \}\]

When u = 0.5, it represents the sum of squares of a vector $x$:

data = [1.0, 2.0, 3.0, 4.0]
target = [0.45, 1.04, 1.51, 1.97]
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, θ)
@variable(model, t)
@expression(model, residuals, θ * data .- target)
@constraint(model, [t; 0.5; residuals] in RotatedSecondOrderCone())
@objective(model, Min, t)
optimize!(model)
value(θ), value(t)

(0.49800000000000544, 0.004979422159618532)

Exponential Cone

The MOI.ExponentialCone is a set of the form:

\[K_{exp} = \{ (x,y,z) \in \mathbb{R}^3 : y \exp (x/y) \le z, y > 0 \}\]

It can be used to model problems involving log and exp.

Exponential

To model $\exp(x) \le z$, use (x, 1, z):

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x == 1.5)
@variable(model, z)
@objective(model, Min, z)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
value(z), exp(1.5)

(4.481654619603659, 4.4816890703380645)

Logarithm

To model $x \le \log(z)$, use (x, 1, z):

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x)
@variable(model, z == 1.5)
@objective(model, Max, x)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
value(x), log(1.5)

(0.4055956586655673, 0.4054651081081644)

Log-sum-exp

To model $t \ge \log\left(\sum e^{x_i}\right)$, use:

N = 3
x0 = rand(N)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[i = 1:N] == x0[i])
@variable(model, t)
@objective(model, Min, t)
@variable(model, u[1:N])
@constraint(model, sum(u) <= 1)
@constraint(model, [i = 1:N], [x[i] - t, 1, u[i]] in MOI.ExponentialCone())
optimize!(model)
value(t), log(sum(exp.(x0)))

(1.4726849472138905, 1.472772274699489)

Entropy

The entropy maximization problem consists of maximizing the entropy function, $H(x) = -x\log{x}$ subject to linear inequality constraints.

\[\begin{aligned} & \max & - \sum_{i=1}^n x_i \log x_i \\ & \;\;\text{s.t.} & \mathbf{1}^\top x = 1 \\ & & Ax \leq b \end{aligned}\]

We can model this problem using an exponential cone by using the following transformation:

\[t\leq -x\log{x} \iff t\leq x\log(1/x) \iff (t, x, 1) \in K_{exp}\]

Thus, our problem becomes,

\[\begin{aligned} & \max & 1^Tt \\ & \;\;\text{s.t.} & Ax \leq b \\ & & 1^T x = 1 \\ & & (t_i, x_i, 1) \in K_{exp} && \forall i = 1 \ldots n \\ \end{aligned}\]

m, n = 10, 15
A, b = randn(m, n), rand(m, 1)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b)
@constraint(model, [i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone())
optimize!(model)
objective_value(model)

2.6632061721520914

The MOI.ExponentialCone has a dual, the MOI.DualExponentialCone, that offers an alternative formulation that can be more efficient for some formulations.

There is also the MOI.RelativeEntropyCone for explicitly encoding the relative entropy function

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x[1:n])
@objective(model, Max, -t)
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b)
@constraint(model, [t; ones(n); x] in MOI.RelativeEntropyCone(2n + 1))
optimize!(model)
objective_value(model)

2.663173643175568

PowerCone

The MOI.PowerCone is a three-dimensional set parameterized by a scalar value α. It has the form:

\[K_{p} = \{ (x,y,z) \in \mathbb{R}^3 : x^{\alpha} y^{1-\alpha} \ge |z|, x \ge 0, y \ge 0 \}\]

The power cone permits a number of reformulations. For example, when $p > 1$, we can model $t \ge x^p$ using the power cone $(t, 1, x)$ with $\alpha = 1 / p$. Thus, to model $t \ge x^3$ with $x \ge 0$

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x >= 1.5)
@constraint(model, [t, 1, x] in MOI.PowerCone(1 / 3))
@objective(model, Min, t)
optimize!(model)
value(t), value(x)

(3.3747548357745134, 1.4999552263154623)

The MOI.PowerCone has a dual, the MOI.DualPowerCone, that offers an alternative formulation that can be more efficient for some formulations.

P-Norm

The p-norm $||x||_p = \left(\sum\limits_{i} |x_i|^p\right)^{\frac{1}{p}}$ can be modeled using MOI.PowerCones. See the Mosek Modeling Cookbook for the derivation.

function p_norm(x::Vector, p)
    N = length(x)
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, r[1:N])
    @variable(model, t)
    @constraint(model, [i = 1:N], [r[i], t, x[i]] in MOI.PowerCone(1 / p))
    @constraint(model, sum(r) == t)
    @objective(model, Min, t)
    optimize!(model)
    return value(t)
end

x = rand(5);
LinearAlgebra.norm(x, 4), p_norm(x, 4)

(0.9316922467512209, 0.9316875994951573)

Positive Semidefinite Cone

The set of positive semidefinite matrices (PSD) of dimension $n$ form a cone in $\mathbb{R}^n$. We write this set mathematically as:

\[\mathcal{S}_{+}^n = \{ X \in \mathcal{S}^n \mid z^T X z \geq 0, \: \forall z\in \mathbb{R}^n \}.\]

A PSD cone is represented in JuMP using the MOI sets PositiveSemidefiniteConeTriangle (for upper triangle of a PSD matrix) and PositiveSemidefiniteConeSquare (for a complete PSD matrix). However, it is preferable to use the PSDCone shortcut as illustrated below.

Example: largest eigenvalue of a symmetric matrix

Suppose $A$ has eigenvalues $\lambda_{1} \geq \lambda_{2} \ldots \geq \lambda_{n}$. Then the matrix $t I-A$ has eigenvalues $t-\lambda_{1}, t-\lambda_{2}, \ldots, t-\lambda_{n}$. Note that $t I-A$ is PSD exactly when all these eigenvalues are non-negative, and this happens for values $t \geq \lambda_{1}$. Thus, we can model the problem of finding the largest eigenvalue of a symmetric matrix as:

\[\begin{aligned} \lambda_{1} = \min t \\ \text { s.t. } t I-A \succeq 0 \end{aligned}\]

A = [3 2 4; 2 0 2; 4 2 3]
I = Matrix{Float64}(LinearAlgebra.I, 3, 3)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* I - A in PSDCone())
optimize!(model)
objective_value(model)

8.000003377698668

GeometricMeanCone

The MOI.GeometricMeanCone is a cone of the form:

\[K_{geo} = \{ (t, x) \in \mathbb{R}^n : x \ge 0, t \le \sqrt[n-1]{x_1 x_2 \cdots x_{n-1}} \}\]

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:4])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in MOI.GeometricMeanCone(5))
optimize!(model)
value(t), value.(x)

(0.0, [0.25000000110304893, 0.250000001103052, 0.25000000110305065, 0.2500000011030498])

RootDetCone

The MOI.RootDetConeSquare is a cone of the form:

\[K = \{ (t, X) \in \mathbb{R}^{1+d^2} : t \le \det(X)^{\frac{1}{d}} \}\]

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; vec(X)] in MOI.RootDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))

(2.236116364743911, 2.2360679574027222)

If X is symmetric, then you can use MOI.RootDetConeTriangle instead. This can be more efficient because the solver does not need to add additional constraints to ensure X is symmetric.

When forming the function, use triangle_vec to obtain the column-wise upper triangle of the matrix as a vector in the order that JuMP requires.

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; triangle_vec(X)] in MOI.RootDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))

(2.2361792479821463, 2.2360679863316992)

LogDetCone

The MOI.LogDetConeSquare is a cone of the form:

\[K = \{ (t, u, X) \in \mathbb{R}^{2+d^2} : t \le u \log(\det(X / u)) \}\]

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; u; vec(X)] in MOI.LogDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))

(1.497920445822069, 1.4978661006407137)

If X is symmetric, then you can use MOI.LogDetConeTriangle instead. This can be more efficient because the solver does not need to add additional constraints to ensure X is symmetric.

When forming the function, use triangle_vec to obtain the column-wise upper triangle of the matrix as a vector in the order that JuMP requires.

model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; u; triangle_vec(X)] in MOI.LogDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))

(1.4979204481049417, 1.497866101098124)

Other Cones and Functions

For other cones supported by JuMP, check out the MathOptInterface Manual.