Solving SDP problem is 100x slower than Matlab CVX

Specific Domains Optimization (Mathematical)
1

I’m trying to solve the following SDP problem
\begin{align*} maximize \ & \sum_j s_j\\ \text{subject to } \ & 2\Sigma - diag(s) \succeq 0,\\ & 0 \le s_j \le 1 \end{align*}

Julia’s SCS solver with Convex.jl is 100x slower than Matlab’s CVX, and it seems like this is because the latter solves the dual problem (which requires a lot less iteration?). How can I get Julia’s solver to be as fast as Matlab?

MWE:

using LinearAlgebra
using Random
using ToeplitzMatrices
using MATLAB
using Convex
using SCS

# simulate data
Random.seed!(2022)
n = 100
p = 200
ρ = 0.4
Sigma = Matrix(SymmetricToeplitz(ρ.^(0:(p-1))))

Solve in Julia (109 seconds):

function julia_sdp(Sigma::Matrix)
    p = size(Sigma, 1)
    svar = Variable(p)
    problem = maximize(sum(svar), svar ≥ 0, 1 ≥ svar, 2*Sigma - Diagonal(svar) == Semidefinite(p))
    solve!(problem, SCS.Optimizer())
    return evaluate(svar)
end
@time s_julia = julia_sdp(Sigma) # compile
@time s_julia = julia_sdp(Sigma) 

------------------------------------------------------------------
	       SCS v3.0.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 40201, constraints m: 80401
cones: 	  z: primal zero / dual free vars: 59901
	  l: linear vars: 400
	  s: psd vars: 20100, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, warm_start: 0
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
	  nnz(A): 100701, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 8.69e+00  9.98e-01  1.69e+03 -1.09e+03  1.00e-01  2.40e-02 
   250| 6.13e-03  4.90e-04  2.35e-01 -1.75e+02  5.22e-01  4.40e+00 
   500| 2.62e-03  1.79e-04  1.02e-01 -1.73e+02  5.22e-01  9.49e+00 
   750| 1.72e-03  9.00e-05  6.30e-02 -1.73e+02  5.22e-01  1.51e+01 
  1000| 1.36e-03  6.81e-05  5.18e-02 -1.73e+02  5.22e-01  2.17e+01 
  1250| 1.02e-03  5.15e-05  3.82e-02 -1.73e+02  5.22e-01  2.85e+01 
  1500| 8.86e-04  4.53e-05  3.45e-02 -1.73e+02  5.22e-01  3.61e+01 
  1750| 8.61e+00  4.50e+00  4.15e-01 -1.75e+02  5.22e-01  4.42e+01 
  2000| 6.42e-04  3.26e-05  2.48e-02 -1.72e+02  5.22e-01  5.25e+01 
  2250| 5.36e-04  2.66e-05  1.99e-02 -1.72e+02  5.22e-01  6.06e+01 
  2500| 5.53e-04  2.89e-05  2.20e-02 -1.72e+02  5.22e-01  6.82e+01 
  2750| 4.47e-04  2.23e-05  1.70e-02 -1.72e+02  5.22e-01  7.48e+01 
  3000| 4.60e-04  2.41e-05  1.84e-02 -1.72e+02  5.22e-01  8.10e+01 
  3250| 3.65e-04  1.81e-05  1.37e-02 -1.72e+02  5.22e-01  8.70e+01 
  3500| 3.71e-04  1.92e-05  1.47e-02 -1.72e+02  5.22e-01  9.29e+01 
  3700| 2.99e-04  8.99e-05  9.86e-03 -1.72e+02  5.22e-01  9.76e+01 
------------------------------------------------------------------
status:  solved
timings: total: 9.77e+01s = setup: 8.60e-02s + solve: 9.76e+01s
	 lin-sys: 1.09e+01s, cones: 8.10e+01s, accel: 7.57e-01s
------------------------------------------------------------------
objective = -172.217851
------------------------------------------------------------------
109.175072 seconds (891.39 k allocations: 88.280 MiB, 0.05% gc time, 0.03% compilation time)

Solve with Matlab’s CVX (1.1 seconds)

function matlab_sdp(Sigma::Matrix)
    p = Float64(size(Sigma, 1))
    @mput Sigma p
    mat"""
    cvx_begin
        variable s(p);
        maximize(sum(s));
        2*Sigma-diag(s) == semidefinite(p);
        0 <= s <= 1;
    cvx_end
    """
    @mget s
    return s
end
@time s_matlab = matlab_sdp(Sigma) # compile
@time s_matlab = matlab_sdp(Sigma)

>> >> >> >> >> >> >> >>  
Calling SDPT3 4.0: 20500 variables, 200 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 200
 dim. of sdp    var  = 200,   num. of sdp  blk  =  1
 dim. of linear var  = 400
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|1.9e+02|1.6e+01|1.5e+06| 8.400000e+04  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|0.524|4.4e-05|7.8e+00|6.2e+05| 4.042192e+04 -1.699825e+03| 0:0:00| chol  1  1 
 2|1.000|1.000|5.2e-05|5.9e-02|3.7e+04| 3.332107e+04 -6.961080e-01| 0:0:00| chol  1  1 
 3|0.990|1.000|1.2e-06|1.8e-02|5.9e+02| 5.836344e+02  1.060232e+01| 0:0:00| chol  1  1 
 4|1.000|1.000|9.2e-08|1.8e-03|2.4e+02| 3.333237e+02  9.407462e+01| 0:0:00| chol  1  1 
 5|1.000|0.807|5.8e-10|4.9e-04|7.2e+01| 2.405317e+02  1.682483e+02| 0:0:00| chol  1  1 
 6|0.798|1.000|1.0e-10|1.8e-05|3.0e+01| 1.946140e+02  1.646145e+02| 0:0:00| chol  1  1 
 7|0.888|0.866|2.1e-11|3.9e-06|3.8e+00| 1.754671e+02  1.716357e+02| 0:0:00| chol  1  1 
 8|0.894|0.907|3.7e-12|5.3e-07|1.2e+00| 1.730688e+02  1.718583e+02| 0:0:00| chol  1  1 
 9|0.845|1.000|1.2e-12|1.8e-08|3.3e-01| 1.722023e+02  1.718724e+02| 0:0:00| chol  1  1 
10|0.933|0.938|8.0e-14|2.8e-09|2.5e-02| 1.719021e+02  1.718774e+02| 0:0:00| chol  1  1 
11|0.940|0.993|1.2e-14|2.0e-11|1.7e-03| 1.718792e+02  1.718775e+02| 0:0:00| chol  1  1 
12|0.964|0.988|1.0e-14|1.2e-12|5.9e-05| 1.718776e+02  1.718776e+02| 0:0:01| chol  1  1 
13|1.000|1.000|9.6e-12|1.0e-12|7.5e-06| 1.718776e+02  1.718776e+02| 0:0:01| chol  1  1 
14|1.000|1.000|3.8e-13|1.5e-12|1.3e-07| 1.718776e+02  1.718776e+02| 0:0:01|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 14
 primal objective value =  1.71877551e+02
 dual   objective value =  1.71877551e+02
 gap := trace(XZ)       = 1.33e-07
 relative gap           = 3.87e-10
 actual relative gap    = 3.86e-10
 rel. primal infeas (scaled problem)   = 3.78e-13
 rel. dual     "        "       "      = 1.50e-12
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 2.0e+02, 1.2e+01, 2.7e+01
 norm(A), norm(b), norm(C) = 2.5e+01, 1.5e+01, 3.7e+01
 Total CPU time (secs)  = 0.56  
 CPU time per iteration = 0.04  
 termination code       =  0
 DIMACS: 2.9e-12  0.0e+00  1.9e-11  0.0e+00  3.9e-10  3.9e-10
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +171.878

  1.157232 seconds (17.41 k allocations: 957.003 KiB, 1.95% compilation time)

Check answer is the same:

[s_julia s_matlab]
200×2 Matrix{Float64}:
 0.999983  1.0
 0.97131   0.971428
 0.811161  0.811429
 0.875385  0.875428
 0.849962  0.849829
 0.86042   0.860068
 0.85656   0.855973
 0.858421  0.857611
 0.857984  0.856956
 0.858443  0.857218
 0.858523  0.857113
 0.858736  0.857155
 0.858881  0.857138
 ⋮

Any help in speeding up the Julia code would be greatly appreciated.

2

I’d try the Hypatia.jl solver, I usually have much better luck with SDPs using that, in terms of speed, accuracy and the number of problems it solves successfully. I have a feeling that SCS doesn’t work particularly well with some of the transformations MathOptInterface is doing because it fails so often when used with Convex.jl

2 Likes
3

Try https://github.com/jump-dev/Dualization.jl:

solve!(problem, () -> Dualization.dual_optimizer(SCS.Optimizer))
2 Likes
4

Yep. Something like 20 times faster. With Hypatia it takes some 9 seconds compared to some 200 seconds using SCS on my old laptop. Here is the complete plug’n’play code, just in case (but I could not resist and I had to rename and reorganize things a bit, but just related to style):

using Convex
using Hypatia
using Random
using SCS
using ToeplitzMatrices

Random.seed!(2022)
n = 100
p = 200
ρ = 0.4
Σ = Matrix(SymmetricToeplitz(ρ.^(0:(p-1))))

function julia_sdp(Σ::Matrix)
    s = Variable(size(Σ,1),Positive())
    add_constraint!(s, s ≤ 1)
    constraint = 2*Σ - diagm(s) ⪰ 0
    problem = maximize(sum(s),constraint)
    #solve!(problem, SCS.Optimizer; silent_solver = true)
    #solve!(problem, () -> Dualization.dual_optimizer(SCS.Optimizer))
    solve!(problem, Hypatia.Optimizer; silent_solver = true)
    return evaluate(s)  
end

@time s_julia = julia_sdp(Σ) # compile
@time s_julia = julia_sdp(Σ)
4 Likes
5

On my mac, Hypatia.jl code above solves the original problem in ~3 seconds which is good enough for my purposes. Thank you very much!

For the approach using Dualization.jl, I’m getting roughly 9 seconds, but it seems to be solving to greater accuracy than Matlab’s CVX.

function julia_dual(Sigma::Matrix)
    p = size(Sigma, 1)
    svar = Variable(p)
    problem = maximize(sum(svar), svar ≥ 0, 1 ≥ svar, 2*Sigma - Diagonal(svar) == Semidefinite(p))
    solve!(problem, Dualization.dual_optimizer(SCS.Optimizer))
    return evaluate(svar)
end
@time s_dual = julia_dual(Sigma) # compile
@time s_dual = julia_dual(Sigma)

  9.215603 seconds (2.61 M allocations: 222.878 MiB, 0.49% gc time, 0.32% compilation time)
 11.798738 seconds (2.53 M allocations: 218.782 MiB, 0.38% gc time)
------------------------------------------------------------------
	       SCS v3.0.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 80401, constraints m: 60701
cones: 	  z: primal zero / dual free vars: 40201
	  l: linear vars: 400
	  s: psd vars: 20100, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, warm_start: 0
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
	  nnz(A): 121201, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 2.25e+01  2.00e+00  1.50e+04 -7.12e+03  1.00e-01  2.68e-02 
   250| 4.05e-03  2.54e-04  1.95e-02  1.73e+02  1.00e-01  4.25e+00 
   475| 1.72e-04  2.62e-04  8.63e-03  1.72e+02  1.02e+00  8.69e+00 
------------------------------------------------------------------
status:  solved
timings: total: 8.82e+00s = setup: 1.28e-01s + solve: 8.69e+00s
	 lin-sys: 1.07e+00s, cones: 7.06e+00s, accel: 8.35e-02s
------------------------------------------------------------------
objective = 172.107787
------------------------------------------------------------------
------------------------------------------------------------------
	       SCS v3.0.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 80401, constraints m: 60701
cones: 	  z: primal zero / dual free vars: 40201
	  l: linear vars: 400
	  s: psd vars: 20100, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, warm_start: 0
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
	  nnz(A): 121201, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 2.25e+01  2.00e+00  1.50e+04 -7.12e+03  1.00e-01  3.85e-02 
   250| 4.05e-03  2.54e-04  1.95e-02  1.73e+02  1.00e-01  5.77e+00 
   475| 1.72e-04  2.62e-04  8.63e-03  1.72e+02  1.02e+00  1.11e+01 
------------------------------------------------------------------
status:  solved
timings: total: 1.13e+01s = setup: 1.88e-01s + solve: 1.11e+01s
	 lin-sys: 1.38e+00s, cones: 9.03e+00s, accel: 1.04e-01s
------------------------------------------------------------------
objective = 172.107787
------------------------------------------------------------------

I will try to see how to adjust convergence tolerance

2 Likes
6

By the way, I think you can use SCS with CVX if you want a comparison of just the modeling layer (Convex vs CVX) or call SDPT3 from Julia.

4 Likes