Hello, julia community. Recently, I aim to solve a large scale optimization problem, which contains 10M exp cone. When I try to model it by using JUMP, it is too slow.
@constraint(jump_model, con[i=1:m], [x[m*n + 2*i - 1], 1.0, x[m*n + 2*i]] in MOI.ExponentialCone())
Hence, I try to use scs inner api to solve it directly.
I try to use the following code:
function fisher_market_old(m, n)
rng = Random.MersenneTwister(1)
w = rand(rng, m)
u = rand(rng, m, n)
b = 0.25*m*ones(n)
m, n = size(u)
model = Model(SCS.Optimizer)
@variable(model, x[1:n, 1:m] >= 0)
@variable(model, p[1:2 * m])
@objective(model, Min, -sum(w[i] * p[2(i-1)+1] for i in 1:m))
for j in 1:n
@constraint(model, sum(x[j, i] for i in 1:m) == b[j])
end
for i in 1:m
@constraint(model, p[2(i-1)+1] == sum(u[i, j] * x[j, i] for j in 1:n))
@constraint(model, [p[2(i-1)+1], 1.0, p[2(i-1)+2]] in MOI.ExponentialCone())
end
JLD2.save("./fisher_market_data/fisher_market_m_$(m)_n_$(n)_pure_random.jld2", "u", u, "w", w, "b", b)
optimize!(model)
println("Optimal solution found.")
optimal_value = -objective_value(model)
println("Optimal value: ", optimal_value)
end
function fisher_market_load_data_scs_api(m, n)
data = JLD2.load("./fisher_market_data/fisher_market_m_$(m)_n_$(n)_pure_random.jld2")
u = data["u"]
w = data["w"]
btilde = data["b"]
b = zeros(m * n + m + n + 3m)
b[1:n] .= btilde
for i in 1:m
# b[m * n + n + 3(i - 1) + 2] = -1.0
b[m * n + n + m + 3(i - 1) + 2] = -1.0
end
u = vec(u')
sparse_u = sparse(u)
nnz_u = nnz(sparse_u)
println(nnz_u)
row_indices = zeros(Int, nnz_u)
for i = 1:nnz_u
row_indices[i] = (sparse_u.nzind[i] - 1) ÷ n + 1 + n
end
c = zeros(m*n + 2*m)
@views c[m*n+1:2:m*n+2*m-1] .= -w
row1 = 1:(m * n)
row1 = row1 .+ n .+ m
col1 = 1:(m * n)
val1 = fill(1.0, m * n)
rows = repeat(1:n, m)
cols = repeat((0:(m-1)) * n, inner=n) .+ repeat(1:n, m)
values = ones(n * m)
rows = vcat([rows, row1]...)
cols = vcat([cols, col1]...)
val = vcat([values, val1]...)
row2 = row_indices
col2 = sparse_u.nzind
val2 = -sparse_u.nzval
row = vcat([rows, row2]...)
col = vcat([cols, col2]...)
val = vcat([val, val2]...)
row3 = (n+1) : (n+m)
row3 = row3
col3 = m * n .+ 2 .* (1:m) .- 1
val3 = fill(1.0, m)
row = vcat([row, row3]...)
col = vcat([col, col3]...)
val = vcat([val, val3]...)
row4 = vcat([[3*i + 1, 3i+3] for i in 0:(m-1)]...) .+ (m + n)
row4 = row4 .+ (n * m)
col4 = Vector{Integer}(1:2m)
col4 = col4 .+ (m * n)
val4 = fill(1.0, 2m)
row = vcat([row, row4]...)
col = vcat([col, col4]...)
val = vcat([val, val4]...)
A = sparse(row, col, val, m + n + 3m + m *n, m*n + m + m)
A = SparseMatrixCSC(A)
mA, nA = size(A)
println("mA: ", mA)
println("nA: ", nA)
println("A.nzval[end-10:end]: ", A.nzval)
println("A.rowval[end-10:end]: ", A.rowval)
println("A.colptr[end-10:end]: ", A.colptr)
println("c: ", c)
println("b: ", -b)
# dataA = JLD2.load("./A.jld2")
# A = dataA["A"]
sol_res = SCS.scs_solve(
SCS.DirectSolver,
mA,
nA,
A,
sparse(zeros(nA, nA)),
-b,
c,
m + n, # z
m * n, # l
Float64[], # bu
Float64[], # bl
Vector{Integer}([]), # q
Vector{Integer}([]), # s
m,
0,
Float64[],
zeros(nA),
zeros(mA),
zeros(mA)
)
end
fisher_market_old(3, 1)
fisher_market_load_data_scs_api(3, 1)
Now, I observe that the matrix that I constructed is the same as jump model:
In jump interface:
A.nzval: [1.0, -0.00790928339056074, 1.0, 1.0, -0.4886128300795012, 1.0, 1.0, -0.21096820215853596, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
A.rowval: [0, 1, 4, 0, 2, 5, 0, 3, 6, 1, 7, 9, 2, 10, 12, 3, 13, 15]
A.colptr: [0, 3, 6, 9, 11, 12, 14, 15, 17, 18]
In my constructation:
A.nzval: [1.0, -0.00790928339056074, 1.0, 1.0, -0.4886128300795012, 1.0, 1.0, -0.21096820215853596, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
A.rowval: [1, 2, 5, 1, 3, 6, 1, 4, 7, 2, 8, 10, 3, 11, 13, 4, 14, 16]
A.colptr: [1, 4, 7, 10, 12, 13, 15, 16, 18, 19]
The only difference is zeroindex and oneindex. And if I save the data from jump and load it by using the following code, then the result is the same.
in scs moi wrapper:
JLD2.save("A.jld2", "A", A)
in my code load:
dataA = JLD2.load("./A.jld2")
A = dataA["A"]
Hence, my question is how to accelerated building a model with 10M cone or why the scs result can not be reproduce, since zero-index or one-index?