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

Specific Domains Optimization (Mathematical)
1

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!