# Solving SDP problem is 100x slower than Matlab CVX

**URL:** <https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704>\
**Category:** Optimization (Mathematical)\
**Tags:** convex-optimization\
**Created:** [April 19, 2022, 7:31pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704 "2022-04-19T19:31:36Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![biona001](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/biona001/32/16497_2.png) [@biona001](https://discourse.julialang.org/u/biona001)\
**Post date:** [April 19, 2022, 7:31pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/1 "2022-04-19T19:31:37Z")

</div>

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:

```julia
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):

```julia
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)

```julia
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:

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 19, 2022, 8:19pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/2 "2022-04-19T20:19:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 19, 2022, 8:37pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/3 "2022-04-19T20:37:38Z")

</div>

Try [https://github.com/jump-dev/Dualization.jl:](https://github.com/jump-dev/Dualization.jl:)

```julia
solve!(problem, () -> Dualization.dual_optimizer(SCS.Optimizer))

```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [April 19, 2022, 8:50pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/4 "2022-04-19T20:50:16Z")

</div>

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

```julia
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(Σ) 

```

---

<div class="post-metadata">

**Author:** ![biona001](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/biona001/32/16497_2.png) [@biona001](https://discourse.julialang.org/u/biona001)\
**Post date:** [April 19, 2022, 9:03pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/5 "2022-04-19T21:03:48Z")

</div>

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.

```julia
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

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [April 19, 2022, 9:28pm UTC](https://discourse.julialang.org/t/solving-sdp-problem-is-100x-slower-than-matlab-cvx/79704/6 "2022-04-19T21:28:35Z")

</div>

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.