Iteration limit with a (relatively) small size problem

[HubertVilluendas]

Hello! Working with Clarabel, I'm currently investigating a variant of the knapsack problem for which I have a semidefinite approach. Namely the problem is expressed as follow:

$$\left[\begin{array}{rl} \text{minimize} & c^\top \mathrm{diag}(X) + \varphi\left(\lambda,X\right) \\ \text{subject to} & w^\top \mathrm{diag}(X) \geq q\\ & \begin{pmatrix} 1 & \mathrm{diag}(X)^\top\\ \mathrm{diag}(X) & X \end{pmatrix}\succcurlyeq 0\\ &\text{ additionnal constraints} \end{array}\right.$$

where our variable is $X$ a symetric matrix of size $n\times n$. I have an heuristic where I solve this problem using a subproblem for which I add cuts separating fractionnal points.

The problem I encounter is that even for relatively small instances (n = 76) the solver doesn't converge and goes to iteration limit, leaving me with a quite bad solution.

For investigation purpose, I've registered the model for several problematic instances in a .cbf (compressed):

Difficult instances (cbf).zip

Could anyone help me with this problem? Otherwise i just report the issue and give you the data to help the developpers solve for this kind of convergence problems.

Here is the JSON corresponding to the problematic instances, in which there is all the informations about the instances: the number of items $n$, the costs and weights, also the time spent to solve the problem and the solution returned by Clarabel.
ProblematicModel_Clarabel_01.json
ProblematicModel_Clarabel_02.json
ProblematicModel_Clarabel_03.json
ProblematicModel_Clarabel_04.json

and also this is my code building the model (written in julia):

using JuMP
import HiGHS
import Clarabel

#= 
    𝐷𝑒𝑠𝑐𝑟𝑖𝑝𝑡𝑖𝑜𝑛 𝑜𝑓 𝑡ℎ𝑒 𝑝𝑟𝑜𝑏𝑙𝑒𝑚:
    
Version of the Knapsack problem where we add "compactness contraints":
we have n items i ∈ {1,…,n} with costs cᵢ and weights wᵢ. The goal is
to find the selection of items S ⊂ {1,…,n} that 𝐦𝐢𝐧𝐢𝐦𝐢𝐳𝐞𝐬 the cost of
the selection while ensuring that the total weight of the items in S
exceed a parameter q > 0.
The selection must also be 𝐜𝐨𝐦𝐩𝐚𝐜𝐭: given Δ ∈ 𝐍, there may be no gap 
of more than Δ between two selected items.
With a {0,1} formulation, (xᵢ = 1 if and only if i ∈ S), this gives
the 𝐜𝐨𝐦𝐩𝐚𝐜𝐢𝐭𝐲 𝐜𝐨𝐧𝐬𝐭𝐫𝐚𝐢𝐧𝐭:
    ⌊(j-i-1)/Δ⌋(xᵢ + xⱼ - 1) ≤ xᵢ₊₁ + ⋯ + xⱼ₋₁      ∀ i,j, j-i > Δ 

    𝐴𝑝𝑝𝑟𝑜𝑎𝑐ℎ:

We consider a semidefinite relaxation of this problem with penality on
the compacity constraint:

    minimize    c^⊤ diag(X) + φ(λ, X)
    subject to  •   w^⊤ diag(X) ≥ q
                •   [   1   diag(X)^⊤]
                    [diag(X)     X   ] ∈ PSDCone
                • Additionnal strengthening constraints

with φ(λ,⋅) → 𝐑 a linear map.

When dealing with a Knapsack problem, one can add cuts to separate an
optimal fractionnal point x∗ from the {0,1}-Knapsack polytope by adding
"cover inequalities". In the context of the min-Knapsack problem, these
inequalities becomes "𝐦𝐚𝐱𝐢𝐦𝐚𝐥 𝐢𝐧𝐬𝐮𝐟𝐟𝐢𝐜𝐢𝐞𝐧𝐭 𝐬𝐮𝐛𝐬𝐞𝐭 𝐜𝐮𝐭𝐬": if S ⊂ {1,…,n}
is such that ∑_{i∈S} wᵢ < q and forall j ∈ {1,…,n}∖S

    wⱼ + ∑_{i∈S} wᵢ ≥ q

then we say that S is a maximal insufficient subset (MIS) and the follo-
-wing inequality 

    ∑_{i∉S} xᵢ ≥ 1

is valid for the min-Knapsack polytope. Hence we add this cuts to sepa-
-rate optimal fractionnal points.
=#


# _______________________________________________________________________
#
#                   𝐔 𝐒 𝐄 𝐅 𝐔 𝐋   𝐅 𝐔 𝐍 𝐂 𝐓 𝐈 𝐎 𝐍 𝐒
# _______________________________________________________________________

# A function that takes a subset MIS, a vector Weights ∈ 𝐑ⁿ and q ∈ 𝐑 and
# returns true if and only if MIS is indeed a maximal insufficient subset. 
function check_maximality(MIS,n,Weights,q)
    global Maximality = (total_weight(MIS,Weights) < q)
    for j in 1:n
        if j ∉ MIS
            global Maximality = Maximality && (Weights[j] + total_weight(MIS,Weights) >= q)
        end
    end
    return(Maximality)
end

# For practicity reasons, I might need to sum weights over a collection S
# that could be empty in some cases. To avoid errors, we need a function
# that returns ∑_{i∈S} wᵢ with the convention ∑_{i∈∅} wᵢ = 0.
function total_weight(subset,Weights)
    if subset==Set()
        return(0)
    else
        return(sum(Weights[i] for i ∈ subset))
    end
end

# A fonction that takes a 𝑠𝑢𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 subset and build a maximal insufficient
# subset by greedly adding items of maximal weights from {1,…,n}∖S.
function build_MaximalInsufficientSubset(subset,n,q,Weights)
    global count_for_precaution = 0
    while (total_weight(subset,Weights) < q) && (count_for_precaution < n)
        global count_for_precaution += 1
        min_weight = Inf
        min_index = -1
        for i in 1:n
            if i ∉ subset && Weights[i] < min_weight
                min_weight = Weights[i]
                min_index = i
            end
        end
        global LastAdded = min_index
        push!(subset,LastAdded)
    end
    delete!(subset,LastAdded)
    return(subset)
end

# The function that comes into play when considering the penalized version.
function payoff(n,i,j,Delta,lambda)
    if abs(j-i) > Delta
        return(lambda)
    else
        return(0)
    end
end

# Function that builds the semidefinite model of the comapct Knapsack
# with the appropriated penality. The model can also take a subset S
# and add the maximal insufficient subset cut to the model.
function build_model_SDP_penality(n::Int,
    Delta::Int,
    Proportion::Float64,
    Weights::Vector{Float64},
    Costs::Vector{Float64},
    Penalite::String="Lagrangian",lambda=0,
    MaxInsSubInequality::Bool=false, subset=Set())

    MaxQ = sum(Weights[i] for i in 1:n)
    q = Proportion*MaxQ

    model = Model(Clarabel.Optimizer)

    @variable(model, X[1:n, 1:n] in SymmetricMatrixSpace())

    if Penalite == "Lagrangian"
        @objective(model, Min, sum(X[i,i]*Costs[i] for i in 1:n)
        + sum(sum(payoff(n,i,j,Delta,lambda)*(floor((j-i-1)/Delta)*X[i,j] - sum(X[k,k] for k in (i+1):(j-1))) for j in (i+Delta+1):n) for i in 1:(n-Delta-1)))
    elseif Penalite == "MaximalGap"
        @variable(model, tau, lower_bound = 0)
        @objective(model, Min, sum(X[i,i]*Costs[i] for i in 1:n) + tau)
        for i in 1:(n-Delta-1)
            for j in (i+Delta+1):n
                @constraint(model, tau >= payoff(n,i,j,Delta,lambda)*(floor((j-i-1)/Delta)*X[i,j] - sum(X[k,k] for k in (i+1):(j-1))))
            end
        end
    elseif Penalite == "PayEachGap"
        @objective(model, Min, sum(X[i,i]*Costs[i] for i in 1:n)
        + sum(sum(payoff(n,i,j,Delta,lambda)*X[i,j] for j in (i+Delta+1):n) for i in 1:(n-Delta-1)))
        for i in 1:(n-Delta-1)
            for j in (i+Delta+1):n
                @constraint(model, floor((j-i-1)/Delta)*X[i,j] <= sum(X[k,k] for k in (i+1):(j-1)))
            end
        end
    else
        error(
            """
            Encountered an error in building the penalized model: type of penality ($(Penalite)) not recognized.
            Penalities implemented: 
                - "Lagrangian" : min            tr(diag(c)X) + ∑ᵢⱼ λᵢⱼ × (⌊(i-j-1)/Δ⌋ Xᵢⱼ - ∑ₖ Xₖₖ )
                                 subject to     • tr(diag(w)X) ≥ q
                                                • Quadratic inequalities
                                                •   [   1    diag(X)ᵀ]
                                                    [diag(X)    X    ] ≽ 0

                - "MaximalGap" : min            tr(diag(c)X) + λ τ
                                 subject to     • tr(diag(w)X) ≥ q
                                                • Quadratic inequalities
                                                •   [   1    diag(X)ᵀ]
                                                    [diag(X)    X    ] ≽ 0
                                                • τ ≥ ⌊(i-j-1)/Δ⌋ Xᵢⱼ - ∑ₖ Xₖₖ      ∀ i,j ∈ [n]

                - "PayEachGap" : min            tr(diag(c)X) + ∑ᵢⱼ λᵢⱼ × Xᵢⱼ
                                 subject to     • tr(diag(w)X) ≥ q
                                                • Quadratic inequalities
                                                •   [   1    diag(X)ᵀ]
                                                    [diag(X)    X    ] ≽ 0
                                                • ⌊(i-j-1)/Δ⌋ Xᵢⱼ ≤ ∑ₖ Xₖₖ          ∀ i,j ∈ [n] 
            """,
        )
    end

    @constraint(model, sum(X[i,i]*Weights[i] for i in 1:n) >= q)

    Diag_SDP = [X[i,i] for i in 1:n]
    global MyMatrix_SDP= vcat(hcat(1,Diag_SDP'),hcat(Diag_SDP,X))
    @constraint(model, MyMatrix_SDP in PSDCone())

    @constraint(model, sum(Weights[i]*Weights[i]*X[i,i] for i in 1:n)+2*sum(sum(Weights[i]*Weights[k]*X[i,k] for k in (i+1):n) for i in 1:(n-1)) <= sum(Weights[i]*Weights[i] for i in 1:n)*sum(sum(X[i,k] for k in 1:n) for i in 1:n))
    for i in 1:n
        @constraint(model, q*X[i,i] <= sum(Weights[j]*X[i,j] for j in 1:n))
        @constraint(model, sum(Weights[j]*X[i,j] for j in 1:n) <= q*(X[i,i]-1) + sum(X[j,j]*Weights[j] for j in 1:n))
        for j in i:n
            @constraint(model, X[i,j] >= 0)
            @constraint(model, X[i,i] >= X[i,j])
            @constraint(model, X[i,j] >= X[i,i] + X[j,j] -1)
            for k in j:n
                @constraint(model, X[i,k] + X[j,k] <= X[k,k] + X[i,j])
                @constraint(model, X[i,j] + X[i,k] + X[j,k] +1 >= X[i,i] + X[j,j] + X[k,k])
            end
        end
    end

    if MaxInsSubInequality
        if subset==Set()
            global MaxInsufficientSubset = Set([i for i in 1:n])
            while total_weight(MaxInsufficientSubset,Weights) >= q
                global random_idx = rand([i for i in 1:n])
                delete!(MaxInsufficientSubset,random_idx)
            end
            global MaxInsufficientSubset = build_MaximalInsufficientSubset(MaxInsufficientSubset,n,q,Weights)        
        else
            global MaxInsufficientSubset = subset
        end
        if check_maximality(MaxInsufficientSubset,n,Weights,q)
            MISComp = setdiff([i for i in 1:n],MaxInsufficientSubset)
            @constraint(model, sum(X[i,i] for i ∈ MISComp) >= 1)
        end
    end
    
    return(model=model, X=X)
end


# _______________________________________________________________________
#
#   𝐅 𝐔 𝐍 𝐂 𝐓 𝐈 𝐎 𝐍   𝐓 𝐎   𝐁 𝐔 𝐈 𝐋 𝐃   𝐓 𝐇 𝐄   𝐌 𝐀 𝐈 𝐍   𝐌 𝐎 𝐃 𝐄 𝐋
# _______________________________________________________________________

# This function build the main model 
function build_model_SDP_penality_and_cover(n::Int,
    Delta::Int,
    Proportion::Float64,
    Weights::Vector{Float64},
    Costs::Vector{Float64},
    Penalite::String="Lagrangian",lambda=0)

    MaxQ = sum(Weights[i] for i in 1:n)
    q = Proportion*MaxQ

    fractionnalmodel = build_model_SDP_penality(n,Delta,Proportion,Weights,Costs,Penalite,lambda)
    set_optimizer(fractionnalmodel.model, Clarabel.Optimizer)
    optimize!(fractionnalmodel.model)

    FractionnalSolution = [value(fractionnalmodel.X[i,i]) for i in 1:n]

    fractionnality_score = 2*sqrt((sum((FractionnalSolution[i]-floor(FractionnalSolution[i]+0.5))^2 for i in 1:n))/n) 

    if fractionnality_score < 0.01
        return(model=fractionnalmodel.model,
                X=fractionnalmodel.X,
                InternalStatus="Already_Solved",
                Objective=objective_value(fractionnalmodel.model),
                Value=value.(fractionnalmodel.X),
                Time=solve_time(fractionnalmodel.model),
                Status=string(termination_status(fractionnalmodel.model)))
    else
        SeparationProblem = Model(HiGHS.Optimizer)
        @variable(SeparationProblem, alpha_bin[1:n], Bin)
        @objective(SeparationProblem, Min, sum(FractionnalSolution[i]*(1-alpha_bin[i]) for i in 1:n))
        @constraint(SeparationProblem, sum(Weights[i]*alpha_bin[i] for i in 1:n) <= q-1e-5)

        undo_relax = relax_integrality(SeparationProblem)

        optimize!(SeparationProblem)

        global WeCanFindAMaxInsufficientSubset = false

        if objective_value(SeparationProblem)<1    # La valeur de l'objectif relaxé est supérieur à 1, donc l'objectif entier est supérieur à 1 : il n'y a pas de MIS qui sépare x_sep
            undo_relax()
            set_silent(SeparationProblem)
            optimize!(SeparationProblem)
            if objective_value(SeparationProblem)<1    # L'objectif entier est supérieur à 1 : il n'y a pas de MIS qui sépare x_sep
                global WeCanFindAMaxInsufficientSubset = true
                MIS = Set()
                for i in 1:n
                    if value(alpha_bin[i])==1
                        push!(MIS,i)
                    end
                end
            end
        end

        if WeCanFindAMaxInsufficientSubset
            MIS = build_MaximalInsufficientSubset(MIS,n,q,Weights)
            ReturnedModel = build_model_SDP_penality(n,Delta,Proportion,Weights,Costs,"Lagrangian",lambda,true,MIS)
            return(model=ReturnedModel.model, X=ReturnedModel.X, InternalStatus="Not_Solved_Yet")
        else
            return(model=fractionnalmodel.model,
                X=fractionnalmodel.X,
                InternalStatus="Already_Solved",
                Objective=objective_value(fractionnalmodel.model),
                Value=value.(fractionnalmodel.X),
                Time=solve_time(fractionnalmodel.model),
                Status=string(termination_status(fractionnalmodel.model)))
        end
    end
end


# _______________________________________________________________________
#
#  𝐄 𝐗 𝐀 𝐌 𝐏 𝐋 𝐄   𝐖 𝐈 𝐓 𝐇   𝐀   𝐃 𝐈 𝐅 𝐅 𝐈 𝐂 𝐔 𝐋 𝐓   𝐈 𝐍 𝐒 𝐓 𝐀 𝐍 𝐂 𝐄
# _______________________________________________________________________

n = 76
Delta = 1
Proportion = 0.61
Weights = [0.18713555395798373,
            0.17814921319594204,
            0.10146577202651962,
            0.02757808131639905,
            0.004213595335090659,
            0.0006190590302739831,
            0.00021966610751657464,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            0.00021966610751657464,
            0.00021966610751657464,
            0.002616023644061025,
            0.020588705168144403,
            0.0814961258886492,
            0.15079079798705955,
            0.14460020768431972,
            0.0771028037383177,
            0.019190829938493473,
            0.002616023644061025,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5,
            1.9969646137870423e-5]

Costs = [1.8388833910079883,
        1.7172987958747963,
        1.4667085321806315,
        1.8613485745343212,
        1.7224790447290017,
        1.3792444942246807,
        1.9640110973938827,
        1.5822743432404311,
        1.4260750996836382,
        1.6271907453147336,
        1.3158402475803694,
        1.5335647625904443,
        1.4707734686990828,
        1.6474593513023494,
        1.144753335413823,
        1.4756716371605125,
        1.2834789792096073,
        1.1359682628321588,
        1.9033748049752246,
        1.3319994948840725,
        1.6582592449483418,
        1.201199086270169,
        1.1656552150326713,
        1.8911754931060474,
        1.0397609995885824,
        1.6076692995511288,
        1.5705190810941874,
        1.8251845890478742,
        1.3431855035703422,
        1.8447413292922465,
        1.5145139477375658,
        1.3127601255884005,
        1.021067231020306,
        1.9385151352187926,
        1.7753463500166764,
        1.2135684085193286,
        1.6819331065241299,
        1.8477834394667791,
        1.290120961621613,
        1.352759588962427,
        1.9685809242336059,
        1.848623643458187,
        1.7616618887389854,
        1.2324330301416218,
        1.7005521374883386,
        1.2403806651816862,
        1.6698159298245212,
        1.5361148146006798,
        1.3829895629282252,
        1.734908322746779,
        1.0619914955187562,
        1.8828405463470403,
        1.2802851819043495,
        1.523922254391386,
        1.1531689078476064,
        1.613463206214997,
        1.9157429879887,
        1.8140834797015406,
        1.3973268072224116,
        1.0494523024429996,
        1.0913508882671499,
        1.9411005509737749,
        1.7925225538933933,
        1.1440036730841436,
        1.21327435698343,
        1.1594541138860361,
        1.1861005896223145,
        1.4062311188437104,
        1.8075507938766027,
        1.5984865648420516,
        1.965146332258116,
        1.95125077887132,
        1.2211837803842451,
        1.69242078803077,
        1.7836232762765452,
        1.1375915232457083]

lambda = 0.01*sum(Costs[i] for i in 1:n)

MyModel = build_model_SDP_penality_and_cover(n,
        Delta,
        Proportion,
        Weights,
        Costs,
        "Lagrangian",lambda)

if MyModel.InternalStatus == "Not_Solved_Yet"
    optimize!(MyModel.model)
end

[goulart-paul]

When I run this example on my own machine it converges in ~70 iterations, but I am working on the v0.10 code (which should be released soon). I don't think anything really relevant to SDPs changed though.

If you are using this in Julia, then you might have better luck (or at least much shorter runtimes) if you use a different linear solver than the default. See here : https://clarabel.org/stable/julia/linear_solvers/. The Panua Pardiso one is probably fastest but requires a license. MKL and HSL not quite as fast, but still a lot faster than default. You need to be using HSL before using Clarabel with these solvers since they are optional dependencies.

[HubertVilluendas]

Hello!

Thank you very much for your response and for testing my example! It's reassuring to know that with your version of Clarabel, the problem converges in approximately 70 iterations, this gives me confidence that the interior point algorithm provides a reasonable heuristic for solving this type of problem ;)
Thaks you also for your advices, I'll investigate using HSL or get a licence for Pardiso. Thank you again for your support!