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

Specific Domains Optimization (Mathematical)
8

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"