Model building is slow when adding a large number of constraints, but the model necessitates a certain type of refrencing of constraints

Specific Domains Optimization (Mathematical)
,
1

Hi there JuMP-Crew
Given a symmetric matrix Q and the all ones matrix E and a set of vectors v_1,…,v_n called “Vertex_Set” and a set of two-tuples of these vectors (or rather their indices) called “Inclusion_Set”, I want to build the linear programming model

max z
subject to
Q-zE == S,
v_i’Sv_j >= 0, forall (i,j) in Inclusion_Set
S is a symmetric matrix variable
z is a scalar variable

This is a subproblem of an iterative algorithm (taken from An Adaptive Linear Approximation Algorithm for Copositive Programs – Optimization Online), where the Vertex_Set and Inclusion_set are updated in every iteration in a way so that the feasible region increases. I omit the details because they are not important to the question.

My problem is that much of the algorithm’s running time goes into the model creation. My implementation of the model building looks as follows

    function Build_LP_Model( )
    #1) Initialization and problem data-----------------------------------------------------------
    # This would actually be function input, but I put it here for the sake 
    # of self-containment
    n = 300
    Q = randn(n,n)
    Q = Q+Q'
    
    E = ones(n,n)

    Vertex_Set_Dict = Dict{Int64, @NamedTuple{vector::Vector{Float64}, parent1::Int64, parent2::Int64}}(
    i => (
          vector  = LinearAlgebra.Diagonal(ones(n))[i,:], 
          parent1 = i, 
          parent2 = i
          ) 
        for i = 1:n
    );
    Inclusion_Set = [ (i,j) for i=1:n for j = i:n ]
#2) Construct Model-------------------------------------------------------------------      
    println("Constructing LP model")
    @time begin
        model = Model(Gurobi.Optimizer)
        set_silent(model)
        @variable(model, z)
        @variable(model, S[i = 1:n, j = 1:n], Symmetric)
        
        @objective(model, Max, z)
        
        @constraint(model, Q-z*E == S)
        println("Easy Part ends here ")
    end    
    @time begin
        Edge_con = Dict{Tuple{Int64, Int64},ConstraintRef}(
            (i,j) => @constraint(
                        model,
                        Vertex_Set_Dict[i].vector'*S*Vertex_Set_Dict[j].vector >= 0
                        ) for (i,j) in Inclusion_Set
        )
        println("Hard Part (Edge_con) ends here ")   
    end
    return (model, Edge_con, S)

    end

This will produce the output

Constructing LP model
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value ******
Academic license ***** - for non-commercial use only - registered to ******
Easy Part ends here
0.640964 seconds (2.30 M allocations: 159.583 MiB, 75.93% gc time, 8.38% compilation time)
Hard Part (Edge_con) ends here
43.991088 seconds (110.03 M allocations: 7.041 GiB, 2.21% gc
time)
(A JuMP Model
├ solver: Gurobi …

The problem is that building the inequality constraints takes almost a minute. Actually, solving the problem will take less than a fraction of a second. In the original paper, the authors have an implementation in C++ where the model building takes half a second (in 2009) for similarly sized problems (see Table 1 in the paper)

A couple of important notes:

  • In every iteration there is an update on the Vertex_Set (only ever increases) and the Inclusion_Set (may lose and gain elements). Constraints will be added and deleted accordingly.
  • Thus, it is important to keep track of the present constraints in what I call the Edge_Set_Dict.
  • Data from the optimal solution as well as from all the sets mentioned above is used to update the respective sets, so all of it must be tracked as well.

My question is now the following

  • Is there an obvious way in JuMP to speed up the model creation, in other words: Am I doing something stupid here?
  • If there is no way to do it in JuMP, what would be an alternative in Julia and where can I read up on such an alternative?

Thank you for your help!

2

You’re recomputing in the constraint expressions over and over again:
the constraint you compute is essentially tr(S, vi'vj) which means it will need quadratic computation with the size of vs. But You can precompute vi'*S for all necessary is to cut compute to linear time per (i,j) pair.
Your vertex set here is fully dense, but I assume it’s only to make MWE simpler.

Try something like this:

@time begin
        VS = map(1:n) do i # assumes Vertex_Set_Dict is dense
            vi = Vertex_Set_Dict[i].vector
            @expression(model, [dot(vi, c) for c in eachcol(S)])
        end

        Edge_con = Dict(
            (i, j) => let
                vj = Vertex_Set_Dict[j].vector
                @constraint(model, dot(VS[i], vj) >= 0)
            end for (i, j) in Inclusion_Set
        )

problem formulation with this is already below 0.5s even for n = 300 on my computer.

EDIT:
On a freshly updated JuMP just adding those ~45_000 constraints to the model results in 0.172214 seconds (1.40 M allocations: 66.347 MiB, 10.63% gc time)
That’s ~31 allocations per constraint, @blegat is this a number you expect? looks too large to me?

[4076af6c] JuMP v1.23.6
1 Like
3

Is Vertex_Set_Dict[i].vector always a {0,1} entry vector? If so, just keep the non-zero entries:

using JuMP, Gurobi
function Build_LP_Model()
    n = 300
    Q = randn(n, n)
    Q = Q + Q'
    E = ones(n, n)
    vertex_set = Dict(
        i => (vertices = Int[i], parent1 = i,  parent2 = i) for i in 1:n
    )
    inclusion_set = [(i, j) for i in 1:n for j in i:n]
    model = Model(Gurobi.Optimizer)
    set_silent(model)
    @variable(model, z)
    @variable(model, S[i in 1:n, j in 1:n], Symmetric)        
    @objective(model, Max, z)
    @constraint(model, Q - z * E == S)
    @constraint(
        model,
        Edge_con[(i, j) in inclusion_set],
        sum(S[ii, jj] for ii in vertex_set[i].vertices,
                          jj in vertex_set[j].vertices) >= 0
    )
    return model
end
@time model = Build_LP_Model()
model[:Edge_con]
model[:S]
4

@odow: Thank you. No its not, only in the first iteration. But they will be sparse through most of the algorithm.

1 Like
5

@Abulak: Thank you! I think this solves the problem.

1 Like