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