How to efficiently generate constraints


#1

I have a rather large SDP problem which I deconstruct (block diagonalize) to smaller ones; A baby example (~5000×5000) in JuMP-0.18.01 took

215.068401 seconds (124.47 M allocations: 231.684 GiB, 16.09% gc time)

(there are 7229 dense constraints, at some point I spent a lot of time trying to optimize this)

Much to my surprise, upgrading to JuMP-0.18.04 caused this:

905.073732 seconds (893.81 M allocations: 1.343 TiB, 38.38% gc time)

below You can find some code but let me vaguely describe it. The problem itself is very simple:

  • P is a vector of SD blocks (JuMP.Variables)
  • M is a Vector{Array{Float64, 2}} of the same length and same dimensions as P
  • all linear constraints look like lhs == x - λ*u,
  • where
    • lhs = sum(vecdot(M[π], P[π]) for π in 1:length(P)) (its a sum of traces of products)
    • x and u are Float64s

I suspect that lhs thing is not the right way, but I couldn’t find anything in the docs.

So here is the code:

Ns = size.(data.Uπs, 2) #Uπs: 5000×n projection matrix
m = JuMP.Model();

P = Vector{Matrix{JuMP.Variable}}(length(Ns))

for (k,n) in enumerate(Ns)
    P[k] = JuMP.@variable(m, [i=1:n, j=1:n]) # small semidefinite blocks
    JuMP.@SDconstraint(m, P[k] >= 0.0)
end

λ = JuMP.@variable(m, λ) 
JuMP.@objective(m, Max, λ)

info("Adding $(length(data.orbits)) constraints... ")
@time addconstraints!(m,P,λ,X,orderunit, data) # the timing above is from here
function addconstraints!(m::JuMP.Model,
    P::Vector{Matrix{JuMP.Variable}}, λ::JuMP.Variable,
    X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData)

    orderunit_orb = orbit_spvector(orderunit.coeffs, data.orbits)
    X_orb = orbit_spvector(X.coeffs, data.orbits)
    UπsT = [U' for U in data.Uπs]

    cnstrs = constraints(parent(X).pm)
    orb_cnstr = spzeros(Float64, size(parent(X).pm)...)

    M = [Array{Float64}(n,n) for n in size.(UπsT,1)]

    for (t, orbit) in enumerate(data.orbits)
        orbit_constraint!(orb_cnstr, cnstrs, orbit)
        constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims)

        lhs = @expression(m, sum(vecdot(M[π], P[π]) for π in 1:length(P))
        x, u = X_orb[t], orderunit_orb[t]
        JuMP.@constraint(m, lhs == x - λ*u)
    end
end

function constraintLHS!(M, cnstr, Us, Ust, dims)
    for π in 1:endof(Us)
        M[π] = dims[π].*Ust[π]*cnstr*Us[π]
    end
end

#2

Which line got slower when updating JuMP ? It would help if you could give a minimal example so that we could fix the regression.
EDIT: You could time each line with TimerOutputs and post the timing for Julia v0.18.1 and JuMP v0.18.4.


#3

I noticed at one point that there were slowdowns on Julia 0.6 for JuMP 0.18.4 because of changes that were made for compatibility with Julia 1.0. We unfortunately don’t have good tools to help control for this. Have you tried running the code on Julia 1.0?