Excessively slow on large, sparse LP

[mlubin]

Here's code to generate an LP with 6M nonzeros, 2M variables and constraints. When the solver is Clp, it takes about a second to generate the model and ~10 seconds to load it into the solver. With SCS it takes ~15 seconds to generate the model (why does this part depend on the solver?) and it never finishes loading the problem into the solver.

When I interrupt the code, it's usually on this line:

SCS.jl/src/MOI_wrapper.jl

Line 317 in d6cf3f4

A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))

# Copyright 2020 Google LLC.
# SPDX-License-Identifier: Apache-2.0
import MathOptInterface
import Clp
import SCS
import Random

const MOI = MathOptInterface

struct PMedianData
    num_facilities::Int
    num_customers::Int
    num_locations::Int
    customer_locations::Vector{Float64}
end

function generate_moi_problem(model, data::PMedianData)
    facility_variables = MOI.add_variables(model, data.num_locations)
    for v in facility_variables
        MOI.add_constraint(model, v, MOI.GreaterThan(0.0))
        MOI.add_constraint(model, v, MOI.LessThan(1.0))
    end
    assignment_variables = [MOI.add_variable(model) for i in 1:data.num_customers, j in 1:data.num_locations]
    for v in assignment_variables
        MOI.add_constraint(model, v, MOI.GreaterThan(0.0))
        # "Less than 1.0" constraint is redundant.
    end
    objective = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(abs(data.customer_locations[i] - j), assignment_variables[i,j])
        for i in 1:data.num_customers for j in 1:data.num_locations],
            0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    for i in 1:data.num_customers, j in 1:data.num_locations
        # assignment_variables[i,j] <= facility_variables[j]
        MOI.add_constraint(model, MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(1.0, assignment_variables[i,j]),
            MOI.ScalarAffineTerm(-1.0, facility_variables[j])], 0.0),
            MOI.LessThan(0.0))
    end
    for i in 1:data.num_customers
        # sum_j assignment_variables[i,j] = 1
        MOI.add_constraint(model, MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(1.0, assignment_variables[i,j]) for j in 1:data.num_locations],
            0.0),
            MOI.EqualTo(1.0)
        )
    end
    # sum_j facility_variables[j] = M
    MOI.add_constraint(model, MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(1.0, facility_variables), 0.0),
        MOI.EqualTo{Float64}(data.num_facilities)
    )
    return assignment_variables, facility_variables
end

function solve_clp(data::PMedianData; time_limit_sec=Inf)
    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        Clp.Optimizer()), Float64)
    if isfinite(time_limit_sec)
        MOI.set(model, MOI.TimeLimitSec(), time_limit_sec)
    end
    @time x, y = generate_moi_problem(model, data)
    @time MOI.optimize!(model)
    @show MOI.get(model, MOI.ObjectiveValue())
end

function solve_scs(data::PMedianData)
    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        SCS.Optimizer()), Float64)
    @time x, y = generate_moi_problem(model, data)
    @time MOI.optimize!(model)
    @show MOI.get(model, MOI.ObjectiveValue())
end

Random.seed!(10)

num_facilities = 10
num_customers = 2000
num_locations = 1000
data = PMedianData(num_facilities, num_customers, num_locations, rand(num_customers) .* num_locations)
solve_clp(data; time_limit_sec=5)
solve_scs(data)

Note: It looks like ECOS has a very similar code structure and therefore hits the same bottleneck in this example.

[blegat]

Most of the time seems to be taken here:
https://github.com/JuliaLang/julia/blob/b79a6d11a771297839c403b358880fbd54c9f285/stdlib/SparseArrays/src/sparsematrix.jl#L747-L748
allocating the colptr which is huge here as the number of variable is large even if the constraints has a few terms.
In fact calling sparse and then findnz just for removing redundant entries is a bit silly as it's going from I, J, V representation to colptr, rowval, V and then back to I, J, V. I think this was written before we had MOI.Utilities.canonical.

why does this part depend on the solver?

For SCS, the scalar constraints need to be bridged by the Bridges.Constraint.VectorizeBridge into vector constraints. This takes some time as we need to convert scalar functions to vector functions. This time could be decreased with jump-dev/MathOptInterface.jl#523

[kalmarek]

@mlubin can you check the timings after #182 is applied? does this solve/alleviate the issue?
(I've hit this more than once, but was a bit terrified to venture into bridges...)

[mlubin]

I didn't do precise timings, but the model loads in ~10 seconds now whereas previously it never finished.