# Solving SDP problem is 100x slower than Matlab CVX

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)


[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.

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
solve!(problem, () -> Dualization.dual_optimizer(SCS.Optimizer))

2 Likes

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())
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

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

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