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.8.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.30e-01  7.71e-01  1.00e+00  1.33e+00   ------
  1   1.3789e-01   4.0894e-02  9.70e-02  7.31e-02  1.14e-01  5.93e-03  2.32e-01  8.70e-01
  2   8.3910e-02   8.1077e-02  2.83e-03  2.57e-03  3.46e-03  6.34e-04  9.67e-03  9.64e-01
  3   8.6450e-02   8.6422e-02  2.79e-05  2.55e-05  3.38e-05  6.39e-06  9.73e-05  9.90e-01
  4   8.6475e-02   8.6475e-02  2.79e-07  2.55e-07  3.38e-07  6.39e-08  9.73e-07  9.90e-01
  5   8.6475e-02   8.6475e-02  2.79e-09  2.55e-09  3.38e-09  6.39e-10  9.73e-09  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  806μs

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

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.8.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.11e-01  3.73e-01  1.00e+00  1.48e+00   ------
  1   9.5000e-02   4.5267e-02  4.97e-02  7.92e-02  3.07e-02  6.41e-03  1.30e-01  9.46e-01
  2   8.3793e-02   8.2081e-02  1.71e-03  2.90e-03  1.09e-03  2.10e-04  5.15e-03  9.62e-01
  3   8.6449e-02   8.6432e-02  1.69e-05  2.88e-05  1.08e-05  2.12e-06  5.18e-05  9.90e-01
  4   8.6475e-02   8.6475e-02  1.69e-07  2.88e-07  1.08e-07  2.12e-08  5.18e-07  9.90e-01
  5   8.6475e-02   8.6475e-02  1.69e-09  2.88e-09  1.08e-09  2.12e-10  5.18e-09  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  781μs

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

This page was generated using Literate.jl.

Semidefinite Programming (SDP) in Julia

Clarabel.jl native interface

Using JuMP