Adding second order cone is slow

Specific Domains Optimization (Mathematical)
1

Hi,

I have a conic model that i want to solve. The model looks like this

using MathOptInterface
const MOI = MathOptInterface
using Gurobi
using JuMP
n=200
model = Model(Gurobi.Optimizer)
@variable(model, t[i=1:n, j=1:n] >= 0)
@variable(model, T[i=1:100000])
@constraint(model, [i=eachindex(T)], T[i] == sum([t[rand(UInt)%n + 1, rand(UInt)%n + 1] for _ in 1:100]))
TT = rand(Float64, length(T)) * 100
@variable(model, ratio[eachindex(T)])
@constraint(model, [i=eachindex(T)], ratio[i] >= T[i]/TT[i])
@constraint(model, [i=eachindex(T)], [ratio[i] + T[i], ratio[i] - T[i], 2 *sqrt(TT[i])] in SecondOrderCone())
@time optimize!(model)

If i remove the second order cone constraint the optimize part ends pretty fast

5.020953 seconds (6.52 M allocations: 884.520 MiB, 3.63% gc time)

However when using Gurobi or Cplex converting the second order cone to quadratics constraints takes too long. what’s the best way to speed this up?

P.s using Mosek the whole process doesn’t take a lot of time.
P.P.s ratio is \max\left( \frac{T}{TT}, \frac{TT}{T}\right)

2

Just try

@constraint(model, [i=eachindex(T)], ratio[i] >= TT[i] / T[i])

The SOC formulation with Gurobi is slow because we need to add slack variables for each term, then convert the SOC to quadratic, and we have to do this 100,000 times.

Mosek is fast because it has native support for SOC.

1 Like
3

i decided to run some benchmark with the direct mode

# model = Model(Gurobi.Optimizer)
function f(n, m; op = Gurobi.Optimizer, direct=false)
    op = optimizer_with_attributes(op, MOI.Silent() => true)
    if !direct
        model = JuMP.Model(op)
    else
        model = JuMP.direct_model(MOI.instantiate(op))
    end
    Random.seed!(1984)
    TT = rand(Float64, m)
    @variable(model, t[i=1:n, j=1:n] >= 0)
    @variable(model, T[i=1:m])
    @constraint(model, [i=eachindex(T)], T[i] == sum([t[rand(UInt)%n + 1, rand(UInt)%n + 1] for _ in 1:100]))
    
    @variable(model, ratio[eachindex(T)])
    @constraint(model, [i=eachindex(T)], ratio[i] >= T[i]/TT[i])
    
    if !direct
        @constraint(model, [i=eachindex(T)], [ratio[i] + T[i], ratio[i] - T[i], 2 * sqrt(TT[i])] in SecondOrderCone())
    else
        @variable(model, rp[eachindex(T)])
        @constraint(model, [i=eachindex(T)], rp[i] == ratio[i] + T[i])
        @variable(model, rm[eachindex(T)])
        @constraint(model, [i=eachindex(T)], rm[i] == ratio[i] - T[i])
        @variable(model, sq[eachindex(T)])
        for i in eachindex(T)
            fix(sq[i], 2 * sqrt(TT[i]))
        end
        @constraint(model, [i=eachindex(T)], [rp[i], rm[i], sq[i]] in SecondOrderCone())
    end
    @objective(model, Min, sum(ratio))
    if !direct
        MOI.Utilities.reset_optimizer(model.moi_backend)
        MOI.Utilities.attach_optimizer(model.moi_backend)
    end
    nothing
end

# score = Dict()
for d in [true, false]
    for S in [ Mosek, CPLEX, Gurobi]
        for i in [10, 100, 1000, 10000, 100000]
            if haskey(score, (i, S, d))
                continue
            end
            GC.gc()
            @show i, S, d
            s = time()
            @time f(200, i, op = S.Optimizer, direct=d)
            score[(i, S, d)] = time() - s
            flush(stdout)
        end
    end
end

gives

(i, S, d) = (10, Mosek, true)
  0.044993 seconds (1.08 M allocations: 86.453 MiB)
(i, S, d) = (100, Mosek, true)
  0.058396 seconds (1.11 M allocations: 89.859 MiB, 20.76% gc time)
(i, S, d) = (1000, Mosek, true)
  0.101414 seconds (1.37 M allocations: 123.827 MiB, 11.99% gc time)
(i, S, d) = (10000, Mosek, true)
  0.573692 seconds (4.02 M allocations: 473.221 MiB, 10.35% gc time)
(i, S, d) = (100000, Mosek, true)
  5.419253 seconds (30.48 M allocations: 3.762 GiB, 11.09% gc time)
(i, S, d) = (10, CPLEX, true)
  0.055432 seconds (722.39 k allocations: 52.289 MiB)
(i, S, d) = (100, CPLEX, true)
  0.062187 seconds (742.48 k allocations: 55.140 MiB)
(i, S, d) = (1000, CPLEX, true)
  0.131571 seconds (943.21 k allocations: 83.800 MiB, 11.73% gc time)
(i, S, d) = (10000, CPLEX, true)
  0.718920 seconds (2.95 M allocations: 368.506 MiB, 12.04% gc time)
(i, S, d) = (100000, CPLEX, true)
  7.030515 seconds (23.02 M allocations: 3.158 GiB, 13.35% gc time)
(i, S, d) = (10, Gurobi, true)
  0.074310 seconds (802.43 k allocations: 33.653 MiB)
(i, S, d) = (100, Gurobi, true)
  0.160088 seconds (4.47 M allocations: 147.561 MiB, 19.40% gc time)
(i, S, d) = (1000, Gurobi, true)
  0.947570 seconds (45.58 M allocations: 1.389 GiB, 24.04% gc time)
(i, S, d) = (10000, Gurobi, true)
 21.800440 seconds (902.18 M allocations: 27.117 GiB, 22.70% gc time)
(i, S, d) = (10, Mosek, false)
  0.212848 seconds (1.99 M allocations: 126.750 MiB, 7.86% gc time)
(i, S, d) = (100, Mosek, false)
  0.128150 seconds (1.80 M allocations: 119.043 MiB, 12.32% gc time)
(i, S, d) = (1000, Mosek, false)
  0.199612 seconds (2.18 M allocations: 162.147 MiB, 15.54% gc time)
(i, S, d) = (10000, Mosek, false)
  0.845246 seconds (5.92 M allocations: 629.632 MiB, 14.35% gc time)
(i, S, d) = (10, CPLEX, false)
  0.090568 seconds (1.41 M allocations: 80.558 MiB)
(i, S, d) = (100, CPLEX, false)
  0.145304 seconds (1.43 M allocations: 83.932 MiB, 13.03% gc time)
(i, S, d) = (1000, CPLEX, false)
  0.792863 seconds (1.70 M allocations: 117.808 MiB, 2.05% gc time)
(i, S, d) = (10000, CPLEX, false)
 17.449910 seconds (4.37 M allocations: 489.156 MiB, 0.74% gc time)
(i, S, d) = (10, Gurobi, false)
  0.124415 seconds (1.49 M allocations: 61.948 MiB)
(i, S, d) = (100, Gurobi, false)
  0.250013 seconds (5.15 M allocations: 176.211 MiB, 15.26% gc time)
(i, S, d) = (1000, Gurobi, false)
  1.679482 seconds (44.87 M allocations: 1.381 GiB, 14.56% gc time)
(i, S, d) = (10000, Gurobi, false)
 38.727610 seconds (753.94 M allocations: 22.796 GiB, 11.13% gc time)

image
image722×480 21.5 KB

4

/(::Int64,::VariableRef) is not defined. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.

5

I guess I meant ratio * T >= TT.

6

Continuous model is non-convex -- solving as a MIP.

i don’t think this sounds like a good idea

7

I’ve observed the same behavior, where adding SecondOrderCone constraints takes a (very) large proportion of model build time.

The code for CPLEX and Gurobi’s MOI wrappers are pretty similar.
The culprit, I think, might be this line in Gurobi’s MOI wrapper: the internal Gurobi model will be updated every time a SOC constraint is added. That can add up to a lot of time.

8

switching to a direct model in cplex kinda solves the problem

9

you’re right it calls ret = @grb_ccall(updatemodel, Cint, (Ptr{Cvoid},), model.ptr_model) every update

10

As for the ratio * T >= TT, it’s a rotated second-order cone constraint, under the assumption that ratio and T are non-negative.
If either can take negative values, the constraint is indeed non-convex.

You’ll have to formulate it explicitly as an RSOC constraint though, I don’t think JuMP will do the quadratic → RSOC conversion for you.
In addition, CPLEX’s MOI wrapper does not support RSOC, so you can;t use direct mode there. One way around that is to reformulate the RSOC as a SOC yourself, and use direct mode.

11

is there any benefit in using RSOC instead of a SOC like i do?

12

Sorry, no, your original formulation was fine. :man_shrugging:
I got confused by the ratio * T >= TT later on.

13

The real issue is that Gurobi and CPLEX don’t support SOC and RSOC explicitly. They infer them from quadratic constraints.

Add lower bounds of 0 to T and ratio.

This worked fine for me:

using Gurobi
using JuMP

n = 200
N = 100_000
TT = rand(Float64, N) * 100

model = Model(Gurobi.Optimizer)
@variables(model, begin
    t[1:n, 1:n] >= 0
    T[1:N] >= 0
    ratio[1:N] >= 0
end)
@constraints(model, begin
    [i = 1:N], T[i] == sum([t[rand(UInt)%n + 1, rand(UInt)%n + 1] for _ in 1:100])
    [i = 1:N], ratio[i] >= T[i] / TT[i]
    [i = 1:N], ratio[i] * T[i] >= TT[i]
end)

optimize!(model)
1 Like
14

Interesting. i had no idea! thanks

1 Like