Semidefinite Programming (SDP) in Julia

We will show how to solve the Basic QP example problem both natively in Clarabel.jl and also by solving with Clarabel.jl within either JuMP or Convex.jl.

Clarabel.jl native interface

using Clarabel, SparseArrays, LinearAlgebra

n    = 3
nvec = Int(n*(n+1)/2)

P = spzeros(nvec,nvec)

q = [1.,0.,1.,0.,0.,1.]

A1 = -Diagonal([                         #<-- LHS of SDP constraint
1., sqrt(2), 1., sqrt(2), sqrt(2), 1.
])
A2 = [1. 2(2.) 3. 2(4.) 2(5.) 6.]        #<-- LHS of equality constraint
A  = sparse([A1;A2]);

b = [zeros(nvec);                       #<-- RHS of SDP constraint
1.]                                #<-- RHS of equality constraint

cones =
[Clarabel.PSDTriangleConeT(n),      #<--- for the SDP constraint
Clarabel.ZeroConeT(1)]             #<--- for the equality constraints

settings = Clarabel.Settings()

solver   = Clarabel.Solver()

Clarabel.setup!(solver, P, q, A, b, cones, settings)

result = Clarabel.solve!(solver)
-------------------------------------------------------------
Clarabel.jl v0.5.1  -  Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------

problem:
variables     = 6
constraints   = 7
nnz(P)        = 0
nnz(A)        = 12
cones (total) = 2
: Zero        = 1,  numel = 1
: PSDTriangle = 1,  numel = 6

settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf,  max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
0   7.3529e-02   7.3529e-02  2.78e-17  4.43e-01  7.71e-01  1.00e+00  1.40e+00   ------
1   1.3358e-01   4.3640e-02  8.99e-02  7.15e-02  1.05e-01  5.81e-03  2.31e-01  8.84e-01
2   8.4656e-02   8.2519e-02  2.14e-03  2.02e-03  2.60e-03  4.71e-04  7.55e-03  9.71e-01
3   8.6457e-02   8.6436e-02  2.12e-05  2.01e-05  2.55e-05  4.71e-06  7.56e-05  9.90e-01
4   8.6475e-02   8.6475e-02  2.12e-07  2.01e-07  2.55e-07  4.71e-08  7.56e-07  9.90e-01
5   8.6475e-02   8.6475e-02  2.12e-09  2.01e-09  2.55e-09  4.71e-10  7.56e-09  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 1.23ms

Recover the symmetric matrix X

X = zeros(n,n)
X[triu(ones(Bool,3,3))] .= result.x
X = Symmetric(X)
3×3 Symmetric{Float64, Matrix{Float64}}:
0.0128947  0.0177209  0.0251946
0.0177209  0.0243534  0.0346243
0.0251946  0.0346243  0.0492269

Using JuMP

We can solve the same problem a little more easily by using Clarabel.jl as the backend solver within JuMP. Here is the same problem again. Note that we are not required to model our SDP constraint in triangular form if we are using JuMP.

using Clarabel, JuMP

n = 3
A = [1 2 4;
2 3 5;
4 5 6]

model = JuMP.Model(Clarabel.Optimizer)
set_optimizer_attribute(model, "verbose", true)
set_optimizer_attribute(model, "equilibrate_enable",false)

@variable(model, X[1:n,1:n],PSD)
@constraint(model, tr(A*X) == 1)
@objective(model, Min, tr(X))

optimize!(model)
-------------------------------------------------------------
Clarabel.jl v0.5.1  -  Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------

problem:
variables     = 6
constraints   = 7
nnz(P)        = 0
nnz(A)        = 12
cones (total) = 2
: Zero        = 1,  numel = 1
: PSDTriangle = 1,  numel = 6

settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf,  max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: off, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
0   7.3529e-02   7.3529e-02  4.16e-17  6.19e-01  3.73e-01  1.00e+00  1.54e+00   ------
1   9.0585e-02   4.5965e-02  4.46e-02  7.46e-02  2.80e-02  6.35e-03  1.24e-01  9.56e-01
2   8.4412e-02   8.3119e-02  1.29e-03  2.28e-03  8.35e-04  1.73e-04  4.04e-03  9.68e-01
3   8.6455e-02   8.6442e-02  1.28e-05  2.26e-05  8.29e-06  1.73e-06  4.05e-05  9.90e-01
4   8.6475e-02   8.6475e-02  1.28e-07  2.26e-07  8.29e-08  1.73e-08  4.05e-07  9.90e-01
5   8.6475e-02   8.6475e-02  1.28e-09  2.26e-09  8.29e-10  1.73e-10  4.05e-09  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 1.22ms

Here is the solution again

JuMP.value.(X)
3×3 Matrix{Float64}:
0.0128947  0.0177209  0.0251946
0.0177209  0.0243534  0.0346243
0.0251946  0.0346243  0.0492269