# Convex.jl (+SCS) more than 100\* slower than R counterpart, CVXR + SCS

**URL:** <https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717>\
**Category:** Optimization (Mathematical)\
**Tags:** performance\
**Created:** [December 1, 2019, 6:17pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717 "2019-12-01T18:17:53Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![reumle](https://avatars.discourse-cdn.com/v4/letter/r/5daacb/32.png) [@reumle](https://discourse.julialang.org/u/reumle)\
**Post date:** [December 1, 2019, 6:17pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/1 "2019-12-01T18:17:53Z")

</div>

Hello, I see significantly worse perf in the julia implementation. Also scaling seems linear on R side and quadratic on Julia side…

[MWE below at the bottom, in R and Julia]

### timings

- in R/CVXR/scs, (in seconds)

```julia-auto
# A tibble: 4 x 5
   nCol nRow moment setupTm solveTm
  <dbl> <dbl> <dttm> <dbl> <dbl>
1 36 200 2019-12-01 10:09:56 0.0868 0.0389
2 36 1000 2019-12-01 10:09:56 0.0319 0.0658
3 36 5000 2019-12-01 10:09:56 0.0409 0.265 
4 36 10000 2019-12-01 10:09:56 0.0529 0.676 

```

- in julia/convex.jl/scs.jl

```julia-auto
4×5 DataFrame
│ Row │ nCol │ nRow │ moment1 │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ DateTime │ Float64 │ Float64 │
├─────┼───────┼───────┼─────────────────────────┼─────────┼─────────┤
│ 1 │ 36 │ 200 │ 2019-12-01T10:26:16.791 │ 0.001 │ 0.096 │
│ 2 │ 36 │ 1000 │ 2019-12-01T10:26:16.888 │ 0.0 │ 1.285 │
│ 3 │ 36 │ 5000 │ 2019-12-01T10:26:18.173 │ 0.002 │ 49.869 │
│ 4 │ 36 │ 10000 │ 2019-12-01T10:27:08.044 │ 0.014 │ 225.797 │

```

- Timing is around 300 times longer for the julia package for the largest(but still fairly small) run.

### MWE’s

- Julia
  - This is the code from function “test4”, but `x` is now a variable to optimise, instead of a parameter.
  - I specify below `SCSSolver( linear_solver= SCS.Direct…)` so that both (julia and R) programs run the same algo: `Direct` is the default setting in CVXR. “Indirect” seems faster by a factor of 2, this does not change the larger story..
  - SCS version under R is 2.1.1 ,under Julia 2.0.2.

```julia-auto
df5= DataFrame(nCol = Int64[], nRow= Int64[], moment1=DateTime[], setupTm=Float64[], solveTm=Float64[])

for nRow in [200,1000,5000,10000] #different problem sizes

  m=nRow
  n=36
  moment1=Dates.now()

  s = Variable(m)

  x = Variable(n); # TODO: This is also more interesting with x = Variable(n)

  # We can compare it to minimize norm(A*x-b, 1) in Matlab with variable x

  A = randn(m,n)

  b = randn(m)

  p = Problem(:minimize,sum(s), [A*x - b <= s, A*x - b >= -s, x>-10, x< 10, sum(x) >10])
  moment2= Dates.now()
  solve!(p,SCSSolver( linear_solver= SCS.Direct, max_iters= 10))

  #----------------------------------------- 3 solve problem
  moment3= Dates.now()

  #----------------------------------------- 4 collect timings.
  cc=[n,m,moment1,Dates.value(moment2-moment1)/1000,Dates.value(moment3-moment2)/1000]
  push!(df5,cc)
end

```

- R

```julia-auto
library("scs")	
library("CVXR")	
library("tibble")
library("dplyr")

df3=tibble(nCol=numeric(),nRow=numeric(),moment=as.POSIXct(character()),setupTm=numeric(),solveTm=numeric())
#for (bscale in c(0.1,0.3,1,2))
for (nRow in c(200,1000,5000,10000)){
moment1=Sys.time()

nCol=36

AA = matrix(rnorm(nRow*nCol),ncol=nCol)
b = rnorm(nRow)
s <- Variable( nrow( AA)) # les coeffs d'investissement
x <- Variable( ncol( AA))

objective1 = Minimize( sum(s) )
constraints = list( x< 10, x > -10,
	                     AA %*% x - b <= s,
	                     AA %*% x - b >= -s,
	                     sum( x) > 10)

prob1 = Problem( objective1, constraints)	
	

moment2 = Sys.time()
solut2	=	psolve( prob1,verbose=TRUE,solver="SCS",max_iters=10,warm_start=FALSE)	
#	sc_a =	solut1$getValue(a)
moment3 = Sys.time()

df3=df3 %>% add_row(nCol=nCol,nRow=nRow,moment=moment1,setupTm=moment2-moment1,solveTm=moment3-moment2)
}

```

### session info

- FWIW I have run a chunk where A is non random (hence reproducible) with A[i,j] = 1/(i+j) and CVXR/convex.jl give essentially the same answers.

Any feedback useful, including pointing out errors in my test. Should I open an issue on GitHub? I would like to try and solve this problem! Cheers!

---

<div class="post-metadata">

**Author:** ![reumle](https://avatars.discourse-cdn.com/v4/letter/r/5daacb/32.png) [@reumle](https://discourse.julialang.org/u/reumle)\
**Post date:** [December 1, 2019, 6:22pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/2 "2019-12-01T18:22:37Z")

</div>

the test code is derived from the github repo for convex.jl: [https://github.com/JuliaOpt/Convex.jl/blob/master/test/oldperf.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/test/oldperf.jl), function `test4`

---

<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:** [December 1, 2019, 7:31pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/3 "2019-12-01T19:31:14Z")

</div>

Thanks for the report. With Convex.jl v0.12.6 and SCS.jl v0.6.2, I get the following timings (after running the MWE once to compile, though that should only affect row 1 anyway):

```julia
julia> df5                                                                         
4×5 DataFrame
│ Row │ nCol │ nRow │ moment1 │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ DateTime │ Float64 │ Float64 │
───────┼───────┼────────────
│ 1 │ 36 │ 200 │ 2019-12-01T18:51:14.496 │ 0.0 │ 0.45 │
│ 2 │ 36 │ 1000 │ 2019-12-01T18:51:14.946 │ 0.001 │ 0.189 │
│ 3 │ 36 │ 5000 │ 2019-12-01T18:51:15.136 │ 0.078 │ 2.31 │
│ 4 │ 36 │ 10000 │ 2019-12-01T18:51:17.524 │ 0.01 │ 8.136 │

```

That seems much worse than the R timings you got, but much better than the Convex timings you got.

I also tried on the upcoming MathOptInterface branch ([https://github.com/JuliaOpt/Convex.jl/pull/330](https://github.com/JuliaOpt/Convex.jl/pull/330)), and got improved timings:

```julia
julia> df5                                                                
4×5 DataFrame
│ Row │ nCol │ nRow │ moment1 │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ DateTime │ Float64 │ Float64 │
├─────┼───────┼───────┼─────────────────────────┼─────────┼─────────┤
│ 1 │ 36 │ 200 │ 2019-12-01T19:13:03.385 │ 0.0 │ 0.011 │
│ 2 │ 36 │ 1000 │ 2019-12-01T19:13:03.396 │ 0.001 │ 0.165 │
│ 3 │ 36 │ 5000 │ 2019-12-01T19:13:03.562 │ 0.004 │ 0.418 │
│ 4 │ 36 │ 10000 │ 2019-12-01T19:13:03.984 │ 0.004 │ 1.515 │

```

However, I think that the solve time might depend quite a bit on the particular instance you are solving, so I don’t think it’s a good idea to benchmark by generating new random data between R and Julia. It’s also usually a good idea to benchmark by running the problem more than once and to take the minimum time, in case there are other things happening on your computer during some of the runs that slow it down artificially.

With that in mind, I rewrote your MWE in Julia as the following:

```julia
using DelimitedFiles, DataFrames, BenchmarkTools

n = 36
for m in [200,1000,5000,10000]
	A = randn(m,n)
	b = randn(m)
	writedlm("n=$(n)_m=$(m)_A.txt", A)  
	writedlm("n=$(n)_m=$(m)_b.txt", b)  
end

function setup_problem(A, b)
	m, n = size(A)
	s = Variable(m)
  	x = Variable(n)
  	p = minimize(sum(s), [A*x - b <= s, A*x - b >= -s, x>-10, x< 10, sum(x) >10])
end

function solve_problem(p)
	solver = SCSSolver(linear_solver= SCS.Direct, max_iters= 10, verbose=false)
	solve!(p, solver; verbose = false)
end

df= DataFrame(n = Int64[], m = Int64[], setupTm=Float64[], solveTm=Float64[])

for m in [200,1000,5000,10000]
 	A = readdlm("n=$(n)_m=$(m)_A.txt")
	b = readdlm("n=$(n)_m=$(m)_b.txt")  

	setup_time = @belapsed setup_problem($A, $b)
	p = setup_problem(A, b)
	solve_time = @belapsed solve_problem($p)
	push!(df, [n, m, setup_time, solve_time])
end

```

which gave me

```julia
julia> df                                                           
4×4 DataFrame
│ Row │ n │ m │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │
├─────┼───────┼───────┼───────────┼───────────┤
│ 1 │ 36 │ 200 │ 8.411e-6 │ 0.0106288 │
│ 2 │ 36 │ 1000 │ 1.0203e-5 │ 0.0791753 │
│ 3 │ 36 │ 5000 │ 1.7449e-5 │ 1.52196 │
│ 4 │ 36 │ 10000 │ 2.6018e-5 │ 6.14177 │

```

On the upcoming MathOptInterface branch, I got

```julia
julia> df                                                      
4×4 DataFrame
│ Row │ n │ m │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │
──────
│ 1 │ 36 │ 200 │ 7.61e-6 │ 0.00876672 │
│ 2 │ 36 │ 1000 │ 9.47e-6 │ 0.0417753 │
│ 3 │ 36 │ 5000 │ 1.7649e-5 │ 0.314997 │
│ 4 │ 36 │ 10000 │ 2.5681e-5 │ 0.741291 │

```

Note, I used the exact same data between the release version of Convex and the MathOptInterface (MOI) version, to be sure they were solving the same problem. It looks like MOI gives us a good speed up!

To use that yourself, simply add the package by typing `]` to enter package mode, and

```julia
] add https://github.com/ericphanson/Convex.jl#MathOptInterface

```

then switch `SCSSolver` to `SCS.Optimizer` in your code. (This branch will become part of the v0.13 release of Convex.jl).

I suggest you compare between Convex and R using the same matrices `A` and `b` to make sure that the difference in problem instances isn’t a factor, and also to try the MOI branch. Please feel free to open an issue if you still see a big performance gap, especially a scaling one.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 1, 2019, 7:50pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/4 "2019-12-01T19:50:37Z")

</div>

Just out of curiosity: a JuMP solution

```julia
using JuMP
using JuMP.MOIU
using SCS

function test_problem(m, n)
    model = Model()
    @variable(model, s[1:m])
    @variable(model, x[1:n])
    
    A = rand(m,n)
    b = rand(m)
    @constraints model begin
        A*x - b .<= s
        A*x - b .>= -s
        x .<= 10
        x .>= -10
        sum(x) >= 10
    end
    
    @objective(model, Min, sum(s))
    
    return model
end

with_SCS() = with_optimizer(SCS.Optimizer, 
    linear_solver=SCS.Direct,
    max_iters = 10,
    verbose = 0,
    )

let n = 36
    for nRow in [200,1000,5000,10000] #different problem sizes
        @info "" nRow
        @time m = test_problem(nRow, n)
        JuMP.set_optimizer(m, with_SCS())
        @time MOIU.attach_optimizer(m)
        @time optimize!(m)
    end
end

```

gives me

```julia
┌ Info: 
└ nRow = 200
  0.002295 seconds (19.66 k allocations: 2.593 MiB)
  0.012505 seconds (42.31 k allocations: 3.551 MiB)
  0.008903 seconds (25.83 k allocations: 6.964 MiB)
┌ Info: 
└ nRow = 1000
  0.022204 seconds (89.47 k allocations: 12.378 MiB, 36.56% gc time)
  0.035574 seconds (175.88 k allocations: 15.779 MiB, 21.48% gc time)
  0.065067 seconds (117.24 k allocations: 57.169 MiB, 19.22% gc time)
┌ Info: 
└ nRow = 5000
  0.090468 seconds (441.52 k allocations: 61.093 MiB, 23.42% gc time)
  0.142002 seconds (843.78 k allocations: 76.152 MiB, 31.80% gc time)
  0.708513 seconds (589.23 k allocations: 883.075 MiB, 35.21% gc time)
┌ Info: 
└ nRow = 10000
  0.235614 seconds (881.52 k allocations: 121.542 MiB, 31.78% gc time)
  0.364521 seconds (1.68 M allocations: 152.897 MiB, 51.57% gc time)
  1.993709 seconds (1.17 M allocations: 3.210 GiB, 24.94% gc time)

```

but only a fraction of the time is spent in `SCS`:

```julia
┌ Info: 
└ nRow = 10000
  0.328984 seconds (881.52 k allocations: 121.542 MiB, 53.73% gc time)
  0.359167 seconds (1.68 M allocations: 152.897 MiB, 51.34% gc time)
----------------------------------------------------------------------------
        SCS v2.1.1 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 740108
eps = 1.00e-05, alpha = 1.50, max_iters = 10, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 10036, constraints m = 20073
Cones: linear vars: 20073
Setup time: 1.64e-01s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 5.05e+20 3.91e+20 1.00e+00 -3.27e+24 2.83e+23 2.78e+24 7.91e-03 
    10| 1.83e+00 2.41e+00 3.39e-01 3.53e+03 7.15e+03 3.35e-13 5.65e-02 
----------------------------------------------------------------------------
Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate, returning best found solution.
Timing: Solve time: 5.67e-02s
        Lin-sys: nnz in L factor: 1130847, avg solve time: 2.72e-03s
        Cones: avg projection time: 1.33e-05s
        Acceleration: avg step time: 1.66e-03s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.0343e-15, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 1.3550e-17
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 1.8296e+00
dual res: |A'y + c|_2 / (1 + |c|_2) = 2.4124e+00
rel gap: |c'x + b'y| / (1 + |c'x| + |b'y|) = 3.3934e-01
----------------------------------------------------------------------------
c'x = 3527.5601, -b'y = 7151.7760
============================================================================
  1.866334 seconds (1.17 M allocations: 3.210 GiB, 25.50% gc time)

```

---

<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:** [December 1, 2019, 8:25pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/5 "2019-12-01T20:25:24Z")

</div>

> [@abulak](#):
>
> model = Model() @variable(model, s[1:m]) @variable(model, x[1:n]) A = rand(m,n) b = rand(m) @constraints model begin A_x - b .\<= s A_x - b .\>= -s x .\<= 10 x .\>= -10 sum(x) \>= 10 end

You can do slightly better by passing the bounds on `x` as proper variable bounds:

```Julia
    model = Model()
    @variable(model, s[1:m])
    @variable(model, -10 <= x[1:n] <= -10)
    
    A = rand(m,n)
    b = rand(m)
    @constraints model begin
        A*x - b .<= s
        A*x - b .>= -s
        sum(x) >= 10
    end

```

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 1, 2019, 8:36pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/6 "2019-12-01T20:36:23Z")

</div>

> [@odow](#):
>
> @variable(model, s[1:m]) @variable(model, -10 \<= x[1:n] \<= -10)

indeed, that cuts the time of the `optimize!` call by a third!  
Is it SCS, or MOI related? (I never understood what is the real difference between variable bounds and constraints)

---

<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:** [December 1, 2019, 8:46pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/7 "2019-12-01T20:46:15Z")

</div>

> Is it SCS, or MOI related? (I never understood what is the real difference between variable bounds and constraints)

SCS related. Most solvers have a special treatment for variable bounds, so there can be a big performance difference between `l <= x <= u` as bounds, and two linear constraints `1.0 * x >= l` and `1.0 * x <= u`.

As a rule of thumb, you should always specify these types of constraints as variable bounds. If the solver does not support variable bounds, JuMP will convert them to constraints. But JuMP will **not** convert linear constraints `1.0 * x <= u` into variable bounds! Why? Because you we don’t know if you want to modify the `1.0` coefficient in future, or add additional terms to the constraint.

---

<div class="post-metadata">

**Author:** ![reumle](https://avatars.discourse-cdn.com/v4/letter/r/5daacb/32.png) [@reumle](https://discourse.julialang.org/u/reumle)\
**Post date:** [December 2, 2019, 12:32am UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/8 "2019-12-02T00:32:32Z")

</div>

Eric, thanks a lot for the detailed and useful feedback. I ran your code, which ensured exactly the same data were being used Julia-side and R-side. Tonight I just wanted to establish a Baseline, and I will run newer versions of convex tomorrow.

Essentially the timing difference [convex vs CVXR] remains (more or less) unchanged, and also scaling still looks quadratic.

### timings (all in secs)

- Julia  
You see a first run with 200 rows in order to compile all functions, and then 3 independent runs:

```julia-auto
13×4 DataFrame
│ Row │ n │ m │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │
├─────┼───────┼───────┼─────────┼─────────┤
│ 1 │ 36 │ 200 │ 0.001 │ 0.059 │
│ 2 │ 36 │ 200 │ 0.0 │ 0.058 │
│ 3 │ 36 │ 1000 │ 0.0 │ 1.256 │
│ 4 │ 36 │ 5000 │ 0.0 │ 45.586 │
│ 5 │ 36 │ 10000 │ 0.0 │ 218.275 │
│ 6 │ 36 │ 200 │ 0.0 │ 0.062 │
│ 7 │ 36 │ 1000 │ 0.0 │ 1.257 │
│ 8 │ 36 │ 5000 │ 0.0 │ 46.429 │
│ 9 │ 36 │ 10000 │ 0.0 │ 210.443 │
│ 10 │ 36 │ 200 │ 0.001 │ 0.061 │
│ 11 │ 36 │ 1000 │ 0.0 │ 1.27 │
│ 12 │ 36 │ 5000 │ 0.0 │ 46.037 │
│ 13 │ 36 │ 10000 │ 0.0 │ 267.79 │

```

- R/CVXR

```julia-auto
> df3
# A tibble: 12 x 5
    nCol nRow moment setupTm solveTm
   <dbl> <dbl> <dttm> <dbl> <dbl>
 1 36 200 2019-12-02 00:16:45 0.0509 0.0429
 2 36 1000 2019-12-02 00:16:45 0.0329 0.0738
 3 36 5000 2019-12-02 00:16:45 0.0299 0.311 
 4 36 10000 2019-12-02 00:16:46 0.0309 0.923 
 5 36 200 2019-12-02 00:16:47 0.0299 0.0319
 6 36 1000 2019-12-02 00:16:47 0.0409 0.0768
 7 36 5000 2019-12-02 00:16:47 0.0329 0.292 
 8 36 10000 2019-12-02 00:16:48 0.0299 0.845 
 9 36 200 2019-12-02 00:16:49 0.0270 0.0309
10 36 1000 2019-12-02 00:16:49 0.0278 0.0628
11 36 5000 2019-12-02 00:16:49 0.0339 0.295 
12 36 10000 2019-12-02 00:16:50 0.0299 0.653 

```

### code

- Julia  
I ran your code but with a change. In order to keep it comparable to the previous version and also to R I did not use `@belapsed`, and kept using the `Dates.now()` function.

The part which I modified is the final loop. The rest is unchanged. Please find it below.

```julia-auto
using DelimitedFiles, DataFrames, BenchmarkTools

n = 36
for m in [200,1000,5000,10000]
	A = randn(m,n)
	b = randn(m)
	writedlm("n=$(n)_m=$(m)_A.txt", A)  
	writedlm("n=$(n)_m=$(m)_b.txt", b)  
end

function setup_problem(A, b)
	m, n = size(A)
	s = Variable(m)
  	x = Variable(n)
  	p = minimize(sum(s), [A*x - b <= s, A*x - b >= -s, x>-10, x< 10, sum(x) >10])
end

function solve_problem(p)
	solver = SCSSolver(linear_solver= SCS.Direct, max_iters= 10, verbose=false)
	solve!(p, solver; verbose = false)
end

df= DataFrame(n = Int64[], m = Int64[], setupTm=Float64[], solveTm=Float64[])
df2= DataFrame(n = Int64[], m = Int64[], setupTm=Float64[], solveTm=Float64[])

for m in [200,1000,5000,10000]
#	for m in [200]
 	A = readdlm("n=$(n)_m=$(m)_A.txt")
	b = readdlm("n=$(n)_m=$(m)_b.txt")  

	moment1=Dates.now()
	p=setup_problem(A, b)
	moment2= Dates.now()
	
	q = solve_problem(p)
	moment3=Dates.now()
	push!(df2, [n, m, Dates.value(moment2-moment1)/1000, Dates.value(moment3-moment2)/1000])

end

```

- R/CVXR code. Reads Julia txt output for A and b.

```julia-auto
library("scs")	
library("CVXR")	
library("tibble")
library("dplyr")

df3=tibble(nCol=numeric(),nRow=numeric(),moment=as.POSIXct(character()),setupTm=numeric(),solveTm=numeric())
#for (bscale in c(0.1,0.3,1,2))
for (nRow in c(200,1000,5000,10000)){

nCol=36

fDir="C:\\Users\\elm\\.atom\\"

fName=paste0(fDir,"n=",nCol,"_m=",nRow,"_A.txt")
AA = as.matrix(read.table(file=fName,header=FALSE))

fName=paste0(fDir,"n=",nCol,"_m=",nRow,"_b.txt")
b = as.matrix(read.table(file=fName,header=FALSE))

moment1=Sys.time()

s <- Variable( nrow( AA)) 
x <- Variable( ncol( AA))

objective1 = Minimize( sum(s) )
constraints = list( x< 10, x > -10,
	                     AA %*% x - b <= s,
	                     AA %*% x - b >= -s,
	                     sum( x) > 10)

prob1 = Problem( objective1, constraints)	
	

moment2 = Sys.time()
solut2	=	psolve( prob1,verbose=TRUE,solver="SCS",max_iters=10,warm_start=FALSE)	
#	sc_a =	solut1$getValue(a)
moment3 = Sys.time()

df3=df3 %>% add_row(nCol=nCol,nRow=nRow,moment=moment1,setupTm=moment2-moment1,solveTm=moment3-moment2)
}

```

Thanks again.  
Emmanuel.  
PS your code calls `setup_problem` twice? Is it intended? (I am a beginner in Julia)  
PPS current package versions, hardware info.

```julia-auto
julia> pkgs
"Convex" => v"0.12.5"
…
"SCS" => v"0.6.2"
…

julia> VERSION
v"1.2.0"

julia> Sys.cpu_summary()
Intel(R) Core(TM) i5-8400 CPU @ 2.80GHz:

julia> Sys.windows_version()
v"10.0.17763"

```

---

<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:** [December 2, 2019, 1:15am UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/9 "2019-12-02T01:15:13Z")

</div>

Thanks for the update!

I’m was a bit puzzled by the differences in performance we’re seeing, but I think it’s actually coming from using Julia 1.2 (I’m on 1.3). The difference between Convex 0.12.5 and 0.12.6 is minimal and should not affect performance (it was just printing fixes and documentation changes), so that shouldn’t be a problem. And we’re on the same version of SCS.jl.  
I get

```julia
julia> Sys.cpu_summary()
Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz

```

so your computer should be strictly faster (my CPU is several generations behind and a mobile chip, since I’m on a laptop, compared to your desktop CPU). I’m also on MacOS instead of Windows, but I don’t think that should matter much.

**But** when I run the code on Julia 1.2, instead of 1.3, I also get very long solve times for `m = 5000` and `m = 10000`. I don’t really know why the performance would be so much worse on 1.2, but it seems like it is, on my computer at least. Luckily 1.3 was just officially released, so it’s a good time to upgrade!

> PS your code calls `setup_problem` twice? Is it intended? (I am a beginner in Julia)

Yes, so `@belapsed` is a macro from BenchmarkTools.jl that runs the code many times and takes the minimum, returning that time (mimicking Base Julia’s `@elapsed` but for benchmarking, hence the `b`). So from that call, I don’t actually get the problem `p`. So then I call `setup_problem` one more time to actually get a problem instance `p` to use for the next step. It only takes `1e-5` seconds anyway to call `setup_problem` according to my benchmarking 🙂. But that actually makes me realize one difference between the two benchmarking methods we’ve used: in the one I suggested, the same problem `p` is reused in the loop for the solve time benchmark. It’s probably a better idea to use a fresh problem each time like you do (even though we aren’t warmstarting or such, so it shouldn’t really matter). To do that with BenchmarkTools, we can change the command to

```julia
solve_time = @belapsed solve_problem(p) setup=(p=setup_problem($A, $b)) evals=1

```

which tells BenchmarkTools to do 1 “eval” per “sample”, and for each sample, set it up by getting a fresh problem `p`. (We do the `evals=1` so that we really are calling `setup` for each new solve). That doesn’t change much for me though; on Convex 0.12.6 and Julia 1.3, I get

```julia
julia> df                                               
4×4 DataFrame
│ Row │ n │ m │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │
──┼──┼────┼
│ 1 │ 36 │ 200 │ 7.712e-6 │ 0.00848731 │
│ 2 │ 36 │ 1000 │ 9.643e-6 │ 0.0768669 │
│ 3 │ 36 │ 5000 │ 1.7546e-5 │ 1.71201 │
│ 4 │ 36 │ 10000 │ 2.5349e-5 │ 6.58762 │

```

---

<div class="post-metadata">

**Author:** ![reumle](https://avatars.discourse-cdn.com/v4/letter/r/5daacb/32.png) [@reumle](https://discourse.julialang.org/u/reumle)\
**Post date:** [December 2, 2019, 9:07pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/10 "2019-12-02T21:07:54Z")

</div>

I just did what you said(Julia 1.3), and the result is pretty spectacular!

- timings now

```julia
4×4 DataFrame
│ Row │ n │ m │ setupTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │
├─────┼───────┼───────┼─────────┼─────────┤
│ 1 │ 36 │ 200 │ 0.364 │ 7.979 │
│ 2 │ 36 │ 1000 │ 0.0 │ 0.061 │
│ 3 │ 36 │ 5000 │ 0.0 │ 1.005 │
│ 4 │ 36 │ 10000 │ 0.0 │ 3.772 │

```

timings last night

```julia
│ 2 │ 36 │ 200 │ 0.0 │ 0.058 │
│ 3 │ 36 │ 1000 │ 0.0 │ 1.256 │
│ 4 │ 36 │ 5000 │ 0.0 │ 45.586 │
│ 5 │ 36 │ 10000 │ 0.0 │ 218.275 │

```

so a **40 to 60 times speedup**! This is remarkable.😃

Thanks a lot for your help and advice.

I am not done though, I would like to use Julia/convex for such programs, but with `m` (ie #Rows) above 1e6…  
This is where my i5 4210u laptop (+R/CVXR) would start choking . Let me try and see how it goes…

---

<div class="post-metadata">

**Author:** ![reumle](https://avatars.discourse-cdn.com/v4/letter/r/5daacb/32.png) [@reumle](https://discourse.julialang.org/u/reumle)\
**Post date:** [December 2, 2019, 10:58pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/11 "2019-12-02T22:58:07Z")

</div>

Hi, Abulak,Odow.  
Thanks a lot for taking the time to rewrite the problem under JuMP. Here are the timings I get with your code. It sems JuMP scales quadratically (ish) whereas CVXR scales linearly (at least at these sizes). I am going to try to understand why…

- timings, JuMP (2 runs)

```julia
│ Row │ n │ m │ setPbmTm │ setOptTm │ attachOptTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼───────┼────────┼──────────┼──────────┼─────────────┼─────────┤
│ 1 │ 36 │ 300 │ 0.045 │ 0.0 │ 1.225 │ 0.013 │
│ 2 │ 36 │ 1000 │ 0.012 │ 0.0 │ 0.023 │ 0.04 │
│ 3 │ 36 │ 3000 │ 0.035 │ 0.0 │ 0.501 │ 0.182 │
│ 4 │ 36 │ 10000 │ 0.105 │ 0.0 │ 0.134 │ 1.165 │
│ 5 │ 36 │ 30000 │ 0.423 │ 0.0 │ 0.449 │ 8.381 │
│ 6 │ 36 │ 100000 │ 1.611 │ 0.001 │ 1.93 │ 108.514 │
│ 7 │ 36 │ 300 │ 0.004 │ 0.001 │ 0.021 │ 0.018 │
│ 8 │ 36 │ 1000 │ 0.014 │ 0.0 │ 0.031 │ 0.052 │
│ 9 │ 36 │ 3000 │ 0.035 │ 0.001 │ 0.05 │ 0.198 │
│ 10 │ 36 │ 10000 │ 0.11 │ 0.0 │ 0.487 │ 1.307 │
│ 11 │ 36 │ 30000 │ 0.397 │ 0.001 │ 0.481 │ 8.43 │
│ 12 │ 36 │ 100000 │ 1.492 │ 0.001 │ 2.061 │ 104.602 │

```

- Julia code(I followed Odow’s hint, changed the timing collection, and not much else)

```julia
using JuMP
using JuMP.MOIU
using SCS

df2= DataFrame(n = Int64[], m = Int64[], setPbmTm=Float64[], setOptTm=Float64[], attachOptTm=Float64[], solveTm=Float64[])

function test_problem(m, n)
	model = Model()
    @variable(model, s[1:m])
    @variable(model, -10 <= x[1:n] <= -10)
    
    A = rand(m,n)
    b = rand(m)
    @constraints model begin
        A*x - b .<= s
        A*x - b .>= -s
        sum(x) >= 10
    end
    
    @objective(model, Min, sum(s))
    
    return model
end

with_SCS() = with_optimizer(SCS.Optimizer, 
    linear_solver=SCS.Direct,
    max_iters = 10,
    verbose = 0,
    )

let n = 36
    for nRow in [300,1000,3000,10000,30000,100000] #different problem sizes
		moment1= Dates.now()
        m = test_problem(nRow, n)
		moment2= Dates.now()
        JuMP.set_optimizer(m, with_SCS())
		moment3= Dates.now()
        MOIU.attach_optimizer(m)
		moment4= Dates.now()
        optimize!(m)
		moment5= Dates.now()
		push!(df2, [n, nRow, Dates.value(moment2-moment1)/1000, Dates.value(moment3-moment2)/1000, Dates.value(moment4-moment3)/1000, Dates.value(moment5-moment4)/1000])
    end
end

```

- R CVXR timings. They depend fairly linearly from # rows ,CVXR is significantly faster at larger sizes

```julia
# A tibble: 6 x 5
   nCol nRow moment setupTm solveTm
  <dbl> <dbl> <dttm> <dbl> <dbl>
1 36 300 2019-12-02 23:41:26 0.284 0.0349
2 36 1000 2019-12-02 23:41:26 0.0329 0.0868
3 36 3000 2019-12-02 23:41:26 0.0459 0.146 
4 36 10000 2019-12-02 23:41:26 0.0798 0.635 
5 36 30000 2019-12-02 23:41:27 0.112 1.59  
6 36 100000 2019-12-02 23:41:29 0.271 5.53  

```

- R code: Essentially same as my first post

```julia

library("scs")	
library("CVXR")	
library("tibble")
library("dplyr")

df3=tibble(nCol=numeric(),nRow=numeric(),moment=as.POSIXct(character()),setupTm=numeric(),solveTm=numeric())
#for (bscale in c(0.1,0.3,1,2))
for (nRow in c(300,1000,3000,10000,30000,100000)){
moment1=Sys.time()

nCol=36

AA = matrix(rnorm(nRow*nCol),ncol=nCol)
b = rnorm(nRow)
s <- Variable( nrow( AA)) # les coeffs d'investissement
x <- Variable( ncol( AA))

objective1 = Minimize( sum(s) )
constraints = list( x< 10, x > -10,
	                     AA %*% x - b <= s,
	                     AA %*% x - b >= -s,
	                     sum( x) > 10)

prob1 = Problem( objective1, constraints)	
	

moment2 = Sys.time()
solut2	=	psolve( prob1,verbose=TRUE,solver="SCS",max_iters=10,warm_start=FALSE)	
#	sc_a =	solut1$getValue(a)
moment3 = Sys.time()

df3=df3 %>% add_row(nCol=nCol,nRow=nRow,moment=moment1,setupTm=moment2-moment1,solveTm=moment3-moment2)
}

```

So thanks for your help.  
Cheers, Emmanuel.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 3, 2019, 11:02pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/12 "2019-12-03T23:02:22Z")

</div>

I see

```julia
│ Row │ n │ m │ setPbmTm │ setOptTm │ attachOptTm │ solveTm │
│ │ Int64 │ Int64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼───────┼────────┼──────────┼──────────┼─────────────┼─────────┤
│ 1 │ 36 │ 300 │ 1.251 │ 0.751 │ 19.433 │ 1.316 │
│ 2 │ 36 │ 1000 │ 0.012 │ 0.001 │ 0.029 │ 0.063 │
│ 3 │ 36 │ 3000 │ 0.039 │ 0.0 │ 0.076 │ 0.206 │
│ 4 │ 36 │ 10000 │ 0.202 │ 0.0 │ 0.17 │ 1.446 │
│ 5 │ 36 │ 30000 │ 0.525 │ 0.0 │ 0.661 │ 6.919 │
│ 6 │ 36 │ 100000 │ 2.04 │ 0.0 │ 2.503 │ 61.686 │
│ 7 │ 36 │ 300 │ 0.004 │ 0.0 │ 0.02 │ 0.013 │
│ 8 │ 36 │ 1000 │ 0.012 │ 0.0 │ 0.029 │ 0.048 │
│ 9 │ 36 │ 3000 │ 0.04 │ 0.001 │ 0.061 │ 0.196 │
│ 10 │ 36 │ 10000 │ 0.136 │ 0.0 │ 0.164 │ 1.572 │
│ 11 │ 36 │ 30000 │ 0.592 │ 0.001 │ 0.712 │ 6.809 │
│ 12 │ 36 │ 100000 │ 2.016 │ 0.001 │ 2.052 │ 59.533 │

```

on my

```julia
julia> versioninfo()
Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3820QM CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)
Environment:
  JULIA_NUM_THREADS = 4

```

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [December 3, 2019, 11:48pm UTC](https://discourse.julialang.org/t/convex-jl-scs-more-than-100-slower-than-r-counterpart-cvxr-scs/31717/13 "2019-12-03T23:48:23Z")

</div>

after a while of investigation:

```julia
const n=36;
const nRow=30_000;
bm = let m = test_problem(nRow, n)
    JuMP.set_optimizer(m, with_SCS())
    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    @assert bm.state == MOIU.ATTACHED_OPTIMIZER
    bm
end

model = bm.optimizer.model

@time MOIU.allocate_load(model.optimizer, model.model_cache, false);
# 6.330737 seconds (2.27 M allocations: 13.898 GiB, 15.13% gc time)

```

is the function that takes the bulk of time, but I’m frankly too lazy to look into this at the moment… 😉