Convergence to a suboptimal solution

[dhendryc]

Hi,

I want to solve a max mineig problem and Pajarito is not always converging to the correct solution as I get a better solution with a different solver. It seems to me that Pajarito is cutting off solution since its reported optimality gap is 0.0.

I also tried replacing HiGHS with SCIP and Hypatia with Clarabel and while this set-up solves the instance below I can generate other instances there this set-up also converges to a suboptimal solution.

Below is a MWE that reproduces the issue. (I used Julia 1.11, Pajarito v0.8.2, Hypatia v0.7.4., HiGHS v1.19.0)

MWE:

using JuMP
using Pajarito
using HiGHS
using Hypatia
using LinearAlgebra

A = [0.640560125688334 0.8988924651705463 0.6026315293750658 0.09942414283811873; 0.3244420648938249 0.21960977858670605 0.8752427026103645 0.7330721267274701; 0.4284151183931746 0.06037385523635341 0.4479437123418444 0.40800090600302463; 0.3539981006442484 0.11431932418843094 0.6797418029178991 0.3676010970291791; 0.6240687276682546 0.8276760414998834 0.8951460038632557 0.9408622568595247; 0.6093405087159935 0.5475397233252607 0.3753804665727667 0.6283323077663837; 0.8148945670877987 0.8159334877118128 0.42139733688904235 0.4596363701822491; 0.8581761761029205 0.6791489434158131 0.8162784435702374 0.38054865157982964; 0.4037354612337182 0.9234073469913137 0.4925105187447588 0.4202196080406787; 0.0465632242904217 0.588608026827289 0.21073940556559967 0.9357195644138552; 0.46079375823300184 0.9655435124426841 0.33129739228282684 0.8749097390632607; 0.23653413851585414 0.9689717890441853 0.4285687872702951 0.9887116589634328; 0.5616823570640757 0.3700511410277585 0.7705341197413721 0.8850442312936958; 0.8295347567225702 0.9252731672804814 0.11200245309187484 0.654354549533481; 0.6877230841636806 0.36562522511027085 0.5380654173827643 0.33035378144723415; 0.48457722652079827 0.15673193660803542 0.20485906499945394 0.6362942989719393; 0.8574183603562568 0.5761171392212233 0.6599616466199261 0.5833842527846491; 0.5019688218144165 0.4815792427281611 0.5051499066121654 0.7951315576134046; 0.46509258174600643 0.8499021885218158 0.6089184177322454 0.1776254689978578; 0.9139243738764997 0.1512127891717162 0.8357702479228285 0.25057409400489294]
N = 6

# Build E-optimal Pajarito model (based on build_E_pajarito_model)
function build_E_pajarito_model_minimal(m, n, time_limit; verbose=true)
    # Generate data
    #A, N, ub = build_data(seed, m, n, zero_one=zero_one)
    ub = fill(1.0, m)
    
    println("Problem size: m=$m, n=$n, N=$N, sum(ub)=$(sum(ub))")

    # Setup solvers
    # MIP solver 
    oa_solver = optimizer_with_attributes(HiGHS.Optimizer,
        MOI.Silent() => !verbose,
        "mip_feasibility_tolerance" => 1e-8,
        "mip_rel_gap" => 1e-6,
    )
    
    # SDP solver
    conic_solver = optimizer_with_attributes(Hypatia.Optimizer, 
        MOI.Silent() => !verbose,
    )
    
    # Pajarito optimizer
    opt = optimizer_with_attributes(Pajarito.Optimizer,
        "time_limit" => time_limit, 
        "iteration_limit" => 100000,
        "oa_solver" => oa_solver, 
        "conic_solver" => conic_solver,
        "tol_rel_gap" => 5e-2,
        "tol_abs_gap" => 1e-6,
        MOI.Silent() => !verbose,
    )

    # Create model
    model = Model(opt)
    
    # Add variables
    @variable(model, x[1:m])
    set_integer.(x)
    @variable(model, t)
    
    # Constraints
    @constraint(model, sum(x) == N)
    @constraint(model, x in MOI.Nonnegatives(m))
    @constraint(model, x .<= ub)
    
    # Objective: maximize t (minimize negative largest eigenvalue)
    @objective(model, Max, t)

    # PSD constraint: A' * diag(x) * A - t*I ⪰ 0
    # This ensures that the smallest eigenvalue of A' * diag(x) * A is >= t
    # Information matrix: A' * diag(x) * A - t*I ⪰ 0
    info_matrix = [
        @expression(model, 
            (i == j ? -t : 0.0) + sum(A[k, i] * x[k] * A[k, j] for k in 1:m)
        ) for i in 1:n, j in 1:n
    ]
    # Add PSD constraint
    @constraint(model, info_matrix in PSDCone())


    return model, x, t, A, N, ub
end

# Main solving function
function solve_e_optimal_pajarito_minimal(;m=50, n=7, time_limit=300, verbose=true)
    println("=== E-Optimal Design with Pajarito ===")
    println("Parameters: m=$m, n=$n")
    
    # Build model
    model, x, t, A, N, ub = build_E_pajarito_model_minimal(m, n, time_limit, verbose=verbose)
    
    # Solve the model
    println("\nSolving...")
    optimize!(model)
    
    # Extract results
    status = termination_status(model)
    solution_x = value.(x)
    solution_t = value(t)
    solve_time_sec = solve_time(model)
    
    # Get additional solver info
    paja_opt = JuMP.unsafe_backend(model)
    num_iterations = paja_opt.num_iters
    num_cuts = paja_opt.num_cuts
    
    println("\n=== Results ===")
    println("Status: $status")
    println("Optimal t value: $solution_t")
    println("Solve time: $(solve_time_sec) seconds")
    println("Number of iterations: $num_iterations")
    println("Number of cuts: $num_cuts")
    println("Solution x: $(solution_x[1:min(20, length(solution_x))])")
    println("Sum of x: $(sum(solution_x))")
    
    return solution_x, solution_t, status, solve_time_sec
end

# Example usage
if abspath(PROGRAM_FILE) == @__FILE__
    # Run a small example
    println("Running minimal E-optimal design example...")
    solution_x, solution_t, status, solve_time_sec = solve_e_optimal_pajarito_minimal(
        m=size(A, 1),     # m (number of experiments)
        n=size(A, 2),      # n (number of parameters)  
        time_limit=300,    # time_limit
        verbose=true,   # verbose
    )
    
    println("\n=== Compare with a different solver ===")
    println("Solution found with a different solver: 0.373251984739451")
    println("x = [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0]")
end

[odow]

This is probably related to the tolerances of the sub solvers. Have you tried experimenting with different tolerances for Hypatia and HiGHS?

[odow]

Confirmed that changing the tolerance of HiGHS changes the solution:

julia> using JuMP, LinearAlgebra

julia> import HiGHS, Hypatia, Pajarito

julia> function main(mip_feasibility_tolerance)
           opt = optimizer_with_attributes(Pajarito.Optimizer,
               "oa_solver" => optimizer_with_attributes(HiGHS.Optimizer,
                   MOI.Silent() => true,
                   "mip_feasibility_tolerance" => mip_feasibility_tolerance,
               ), 
               "conic_solver" => optimizer_with_attributes(Hypatia.Optimizer,
                   MOI.Silent() => true,
               ),
               MOI.Silent() => false,
           )

           A = [
               0.64056012568833400 0.89889246517054630 0.60263152937506580 0.09942414283811873
               0.32444206489382490 0.21960977858670605 0.87524270261036450 0.73307212672747010
               0.42841511839317460 0.06037385523635341 0.44794371234184440 0.40800090600302463
               0.35399810064424840 0.11431932418843094 0.67974180291789910 0.36760109702917910
               0.62406872766825460 0.82767604149988340 0.89514600386325570 0.94086225685952470
               0.60934050871599350 0.54753972332526070 0.37538046657276670 0.62833230776638370
               0.81489456708779870 0.81593348771181280 0.42139733688904235 0.45963637018224910
               0.85817617610292050 0.67914894341581310 0.81627844357023740 0.38054865157982964
               0.40373546123371820 0.92340734699131370 0.49251051874475880 0.42021960804067870
               0.04656322429042170 0.58860802682728900 0.21073940556559967 0.93571956441385520
               0.46079375823300184 0.96554351244268410 0.33129739228282684 0.87490973906326070
               0.23653413851585414 0.96897178904418530 0.42856878727029510 0.98871165896343280
               0.56168235706407570 0.37005114102775850 0.77053411974137210 0.88504423129369580
               0.82953475672257020 0.92527316728048140 0.11200245309187484 0.65435454953348100
               0.68772308416368060 0.36562522511027085 0.53806541738276430 0.33035378144723415
               0.48457722652079827 0.15673193660803542 0.20485906499945394 0.63629429897193930
               0.85741836035625680 0.57611713922122330 0.65996164661992610 0.58338425278464910
               0.50196882181441650 0.48157924272816110 0.50514990661216540 0.79513155761340460
               0.46509258174600643 0.84990218852181580 0.60891841773224540 0.17762546899785780
               0.91392437387649970 0.15121278917171620 0.83577024792282850 0.25057409400489294
           ]
           m, n = size(A)
           model = Model(opt)
           @variable(model, x[1:m], Bin)
           @variable(model, t)
           @constraint(model, sum(x) == 6)
           @objective(model, Max, t)
           @expression(
               model, 
               Y[i in 1:n, j in 1:n],
               sum(A[k, i] * x[k] * A[k, j] for k in 1:m)
           )
           @constraint(model, con, Symmetric(Y - t * I(n)) in PSDCone())
           optimize!(model)
           display(solution_summary(model))
           return model
       end
main (generic function with 1 method)

julia> model_6 = main(1e-6);
solving continuous relaxation
continuous relaxation status is ALMOST_OPTIMAL
separated 0 rays before imposing integrality
 iter     cuts          obj        bound          gap
continuous subproblem status is OPTIMAL
    1       22  -2.5718e-01  -4.0723e-01   5.8342e-01
continuous subproblem status is OPTIMAL
    2       25  -3.7325e-01  -4.0723e-01   9.1020e-02
continuous subproblem status is OPTIMAL
    3       28  -3.7325e-01  -3.8649e-01   3.5457e-02
continuous subproblem status is OPTIMAL
    4       31  -3.7325e-01  -3.8326e-01   2.6818e-02
continuous subproblem status is OPTIMAL
    5       34  -3.7325e-01  -3.8157e-01   2.2271e-02
┌ Warning: integral solution repeated
└ @ Pajarito ~/.julia/packages/Pajarito/gSNvz/src/algorithms.jl:103
separation cuts could not be added
    6       34  -3.7325e-01  -3.7325e-01   0.0000e+00
objective relative gap 0.0 reached; terminating
OA solver finished with status OPTIMAL, after 0.928272008895874 seconds and 34 cuts
iterative method used 6 iterations

solution_summary(; result = 1, verbose = false)
├ solver_name          : Pajarito
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : OPTIMAL
│ └ objective_bound    : 3.73252e-01
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 3.73252e-01
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 9.28272e-01
  └ node_count         : 21

julia> model_8 = main(1e-8);
solving continuous relaxation
continuous relaxation status is ALMOST_OPTIMAL
separated 0 rays before imposing integrality
 iter     cuts          obj        bound          gap
continuous subproblem status is OPTIMAL
    1       22  -3.3826e-01  -3.7920e-01   1.2103e-01
┌ Warning: integral solution repeated
└ @ Pajarito ~/.julia/packages/Pajarito/gSNvz/src/algorithms.jl:103
separation cuts could not be added
    2       22  -3.3826e-01  -3.3826e-01   0.0000e+00
objective relative gap 0.0 reached; terminating
OA solver finished with status OPTIMAL, after 0.049635887145996094 seconds and 22 cuts
iterative method used 2 iterations

solution_summary(; result = 1, verbose = false)
├ solver_name          : Pajarito
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : OPTIMAL
│ └ objective_bound    : 3.38256e-01
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 3.38256e-01
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 4.96359e-02
  └ node_count         : 1

Perhaps @chriscoey wants to chime in here. I don't know if there are any good rules for choosing the tolerances of the sub solvers when the original model is driving to the edge of a PSD feasible region.

[odow]

If the MIP solver supports callbacks (like Gurobi), then Pajarito will use a single tree method:

julia> using JuMP, LinearAlgebra

julia> import Clarabel, Gurobi, Pajarito

julia> function main()
           opt = optimizer_with_attributes(
               Pajarito.Optimizer,
               "oa_solver" =>
                   optimizer_with_attributes(Gurobi.Optimizer, MOI.Silent() => true), 
               "conic_solver" =>
                   optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true),
           )
           A = [
               0.64056012568833400 0.89889246517054630 0.60263152937506580 0.09942414283811873
               0.32444206489382490 0.21960977858670605 0.87524270261036450 0.73307212672747010
               0.42841511839317460 0.06037385523635341 0.44794371234184440 0.40800090600302463
               0.35399810064424840 0.11431932418843094 0.67974180291789910 0.36760109702917910
               0.62406872766825460 0.82767604149988340 0.89514600386325570 0.94086225685952470
               0.60934050871599350 0.54753972332526070 0.37538046657276670 0.62833230776638370
               0.81489456708779870 0.81593348771181280 0.42139733688904235 0.45963637018224910
               0.85817617610292050 0.67914894341581310 0.81627844357023740 0.38054865157982964
               0.40373546123371820 0.92340734699131370 0.49251051874475880 0.42021960804067870
               0.04656322429042170 0.58860802682728900 0.21073940556559967 0.93571956441385520
               0.46079375823300184 0.96554351244268410 0.33129739228282684 0.87490973906326070
               0.23653413851585414 0.96897178904418530 0.42856878727029510 0.98871165896343280
               0.56168235706407570 0.37005114102775850 0.77053411974137210 0.88504423129369580
               0.82953475672257020 0.92527316728048140 0.11200245309187484 0.65435454953348100
               0.68772308416368060 0.36562522511027085 0.53806541738276430 0.33035378144723415
               0.48457722652079827 0.15673193660803542 0.20485906499945394 0.63629429897193930
               0.85741836035625680 0.57611713922122330 0.65996164661992610 0.58338425278464910
               0.50196882181441650 0.48157924272816110 0.50514990661216540 0.79513155761340460
               0.46509258174600643 0.84990218852181580 0.60891841773224540 0.17762546899785780
               0.91392437387649970 0.15121278917171620 0.83577024792282850 0.25057409400489294
           ]
           m, n = size(A)
           model = Model(opt)
           set_silent(model)
           @variable(model, x[1:m], Bin, start = 1)
           @variable(model, t, start = 0)
           @constraint(model, sum(x) == 6)
           @objective(model, Max, t)
           @constraint(model, con, Symmetric(A'*diagm(x)*A - t*I(n)) in PSDCone())
           optimize!(model)
           display(solution_summary(model))
           return model
       end
main (generic function with 2 methods)

julia> model = main();
Set parameter LicenseID to value 890341
solution_summary(; result = 1, verbose = false)
├ solver_name          : Pajarito
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : OPTIMAL
│ └ objective_bound    : 3.73252e-01
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 3.73252e-01
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 4.19950e-02
  └ node_count         : 53

I think the lesson is to try multiple solvers.

[dhendryc]

I see, thank you! I'll test some more with tightening tolerances.

[odow]

Note that tightening the mip_feasibility_tolerance made the solve worse in this case. I don't know that I understand why.

[dhendryc]

So I tested around a bit more, tightening tolerances and using different MIP solvers. Here is an example with SCIP. Even then changing the tolerance, it converges to the same suboptimal solution and the difference between the solution from Pajarito and the other one is quite substantial.

using JuMP
using Pajarito
using SCIP
using Hypatia
using LinearAlgebra

A = [0.20397479276512875 0.5600940343125469 0.45472524981501117 0.35278726155098583 0.2881693713038024; 0.8906161322634301 0.3335383513927278 0.8141953306761535 0.37554125644589864 0.9965546254596083; 0.8700506073165521 0.40558671075956976 0.5915539781362535 0.9602320903821391 0.002528412650643075; 0.42335436127898696 0.1470437521939535 0.409583252155334 0.9955322862363234 0.3561965507509096; 0.45204419323139544 0.17584310786405166 0.4698068119250204 0.7978224527861346 0.061869171925145405; 0.9938435069453978 0.7607506861442169 0.13546783713561372 0.8436400602066185 0.412351140464992; 0.7888326559680227 0.21660655041080557 0.32358423430652017 0.7694832580788346 0.8668872564364225; 0.027135689021179576 0.3349548147347662 0.7115694885069926 0.6418700549520155 0.08100612763461557; 0.3269857804459574 0.10777278945058633 0.9004654179603605 0.38571749521335286 0.395157495883918; 0.9792698588382195 0.9484728035434385 0.4021090528623319 0.18110220183559056 0.33822476076104013; 0.47225482425801946 0.6999383110475934 0.7434579046476046 0.06515798892285385 0.6441088270000772; 0.7701665490135011 0.4881769438341823 0.5767565597443073 0.573869290723939 0.20832506898287184; 0.727854587866495 0.7755879481131438 0.08760362801288879 0.20088484180907973 0.49137252537724185; 0.1581492940428847 0.5137375049340204 0.11347016740145532 0.5060836042999693 0.16475316149130337; 0.4260373399702665 0.048105183487624936 0.6764412459496885 0.8485116638850635 0.4318280767036963; 0.1657225990804625 0.8157009907353074 0.8717719921326075 0.5241401406144064 0.85025157632051; 0.2352649154510822 0.29237813893064035 0.7828386987122108 0.5543790660828773 0.7538299898689885; 0.921716669609433 0.6939900120670813 0.5966424657601748 0.9150786006579544 0.9558569076471641; 0.5311477811962938 0.48132070682224304 0.06322226331186254 0.17943478772892596 0.004182201317213363; 0.9899313604017254 0.5749617319544268 0.280951123629432 0.6442924765916151 0.4736902161087042; 0.3477776317605237 0.7807396957193539 0.9166393540193443 0.8512578867844222 0.8825343308831509; 0.33864098343590165 0.5504498104779956 0.35727443644448875 0.5401396318763271 0.6447223256297139; 0.7145221805541991 0.1255772751758224 0.2878764812168828 0.7100474288791633 0.643308145452051; 0.29973268231824923 0.8510437636247872 0.6976712939628285 0.9223738043230265 0.6833324422939229; 0.9717592154668062 0.5268980438617441 0.39933851486002003 0.26267184953206313 0.6478661326489296; 0.7369783524265642 0.7463998580598168 0.4885594555708441 0.4849074758070735 0.4840946680594752; 0.16822212020239136 0.7340438357074418 0.5443124901171333 0.3808734950664604 0.2195372744044719; 0.6592078860366614 0.624868249518463 0.6910101221053022 0.3686018090810178 0.6129289148987004; 0.6170457394198341 0.4423415090190661 0.0033886418945704433 0.22528927224136563 0.1165356815132812; 0.6697706034599619 0.778284374323255 0.451937939936528 0.38470268304966204 0.5684069338632726]
N = 7

# Build E-optimal Pajarito model (based on build_E_pajarito_model)
function build_E_pajarito_model_minimal(m, n, time_limit; verbose=true)
    # Generate data
    #A, N, ub = build_data(seed, m, n, zero_one=zero_one)
    ub = fill(1.0, m)
    
    println("Problem size: m=$m, n=$n, N=$N, sum(ub)=$(sum(ub))")

    # Setup solvers
    # MIP solver 
    oa_solver = optimizer_with_attributes(SCIP.Optimizer,
        MOI.Silent() => !verbose,
        "limits/absgap" => 1e-8,
        "limits/gap" => 1e-8,
    )
    
    # SDP solver
    conic_solver = optimizer_with_attributes(Hypatia.Optimizer, 
        MOI.Silent() => !verbose,
    )
    
    # Pajarito optimizer
    opt = optimizer_with_attributes(Pajarito.Optimizer,
        "time_limit" => time_limit, 
        "iteration_limit" => 100000,
        "oa_solver" => oa_solver, 
        "conic_solver" => conic_solver,
        "tol_rel_gap" => 5e-2,
        "tol_abs_gap" => 1e-6,
        MOI.Silent() => !verbose,
    )

    # Create model
    model = Model(opt)
    
    # Add variables
    @variable(model, x[1:m])
    set_integer.(x)
    @variable(model, t)
    
    # Constraints
    @constraint(model, sum(x) == N)
    @constraint(model, x in MOI.Nonnegatives(m))
    @constraint(model, x .<= ub)
    
    # Objective: maximize t (minimize negative largest eigenvalue)
    @objective(model, Max, t)

    # PSD constraint: A' * diag(x) * A - t*I ⪰ 0
    # This ensures that the smallest eigenvalue of A' * diag(x) * A is >= t
    # Information matrix: A' * diag(x) * A - t*I ⪰ 0
    info_matrix = [
        @expression(model, 
            (i == j ? -t : 0.0) + sum(A[k, i] * x[k] * A[k, j] for k in 1:m)
        ) for i in 1:n, j in 1:n
    ]
    # Add PSD constraint
    @constraint(model, info_matrix in PSDCone())


    return model, x, t, A, N, ub
end

# Main solving function
function solve_e_optimal_pajarito_minimal(;m=50, n=7, time_limit=300, verbose=true)
    println("=== E-Optimal Design with Pajarito ===")
    println("Parameters: m=$m, n=$n")
    
    # Build model
    model, x, t, A, N, ub = build_E_pajarito_model_minimal(m, n, time_limit, verbose=verbose)
    
    # Solve the model
    println("\nSolving...")
    optimize!(model)
    
    # Extract results
    status = termination_status(model)
    solution_x = value.(x)
    solution_t = value(t)
    solve_time_sec = solve_time(model)
    
    # Get additional solver info
    paja_opt = JuMP.unsafe_backend(model)
    num_iterations = paja_opt.num_iters
    num_cuts = paja_opt.num_cuts
    
    println("\n=== Results ===")
    println("Status: $status")
    println("Optimal t value: $solution_t")
    println("Solve time: $(solve_time_sec) seconds")
    println("Number of iterations: $num_iterations")
    println("Number of cuts: $num_cuts")
    println("Solution x: $(solution_x[1:min(20, length(solution_x))])")
    println("Sum of x: $(sum(solution_x))")
    
    return solution_x, solution_t, status, solve_time_sec
end

# Example usage
if abspath(PROGRAM_FILE) == @__FILE__
    # Run a small example
    println("Running minimal E-optimal design example...")
    solution_x, solution_t, status, solve_time_sec = solve_e_optimal_pajarito_minimal(
        m=size(A, 1),     # m (number of experiments)
        n=size(A, 2),      # n (number of parameters)  
        time_limit=300,    # time_limit
        verbose=true,   # verbose
    )
    
    println("\n=== Compare with a different solver ===")
    println("Solution found with a different solver: 0.456240835254216")
    println("x = [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]")
end