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

reumle · December 1, 2019, 6:17pm

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)

# 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

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.

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

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!

reumle · December 1, 2019, 6:22pm

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

ericphanson · December 1, 2019, 7:31pm

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> 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), and got improved timings:

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:

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

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

abulak · December 1, 2019, 7:50pm

Just out of curiosity: a JuMP solution

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

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

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

odow · December 1, 2019, 8:25pm

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

    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

abulak · December 1, 2019, 8:36pm

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)

odow · December 1, 2019, 8:46pm

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.

reumle · December 2, 2019, 12:32am

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:

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

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

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.

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> 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"

ericphanson · December 2, 2019, 1:15am

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

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> 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    │

reumle · December 2, 2019, 9:07pm

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

timings now

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

│ 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…

reumle · December 2, 2019, 10:58pm

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)

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

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

# 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


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.

abulak · December 3, 2019, 11:02pm

I see

│ 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> 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

abulak · December 3, 2019, 11:48pm

after a while of investigation:

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…

Topic		Replies	Views
Convex.jl : remarkable increase in perf in dev version (0.13.0) compared to v0.12.6, (now as fast as R/CVXR) Optimization (Mathematical) performance	1	642	December 7, 2019
Solving SDP problem is 100x slower than Matlab CVX Optimization (Mathematical) convex-optimization	5	1185	April 19, 2022
Solving large(?) problem with Convex.jl and SCS Optimization (Mathematical) question , first-steps	13	2525	February 13, 2020
Understanding DCP in Convex.jl and the difference of elapsed time Optimization (Mathematical) convex-optimization , disciplined-convex-p	8	633	December 6, 2021
Improve speed SDP variable entries continuous 0-1 in Convex.jl Optimization (Mathematical)	4	167	February 6, 2024