SLOW_PROGRESS or NUMERICAL_ERROR: a simple but nontrivial SDP example, solved with LP relaxation

Specific Domains Optimization (Mathematical)
, , ,
8

This should be a final perspective on the troublesome SDP problem proposed by @odow

(primal problem) \quad \min_{x, y} \{2y | M \in \mathrm{PSD} \}, where M = [0 y; y x].

The difficulty of this problem lies in the unboundedness of the variable x and y.
When I employ the LP relaxation to tackle the primal problem, it appears that it is a feat to decide the initial artificial bound on x and y鈥攖he free variables. Because of this, I have to add a presolve loop (:blue_circle:).

The full code is as follows (this is a very mature (鈮坉efinitive) version)

import JuMP, LinearAlgebra
import Gurobi; GRB_ENV = Gurobi.Env()
function my_eigmin(A)
    valmin, vecmin = LinearAlgebra.eigen(A, 1:1)
    return first(valmin), vecmin[:, 1]
end
function optimise(model)
    JuMP.optimize!(model)
    JuMP.assert_is_solved_and_feasible(model; allow_local = false, dual = true)
end
CoT = 1e-7 # linear cut-off tolerance, set according to my experience of using Gurobi
function get_result_with(ib) # `ib` is an indispensable artificial initial bound
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(model)
    JuMP.@variable(model, -ib <= y <= ib); JuMP.@variable(model, -ib <= x <= ib)
    JuMP.@expression(model, X, LinearAlgebra.Symmetric([0 y; y x]))
    JuMP.@objective(model, Min, 2 * y)
    while true
        optimise(model)
        Xt = LinearAlgebra.Symmetric(JuMP.value.(X)) # a trial solution that may be infeasible to the primal problem
        _, t = my_eigmin(Xt) # generate a linear cut coefficient at the current trial solution
        if t' * Xt * t <= -CoT # the current trial solution `Xt` can be strongly cut off
            JuMP.@constraint(model, t' * X * t >= 0) # add a linear cut, as a valid substituting cut for the PSDCut
        else # no more linear cuts can be added
            break
        end
    end
    yu = -JuMP.dual(JuMP.UpperBoundRef(y))
    xu = -JuMP.dual(JuMP.UpperBoundRef(x))
    yl =  JuMP.dual(JuMP.LowerBoundRef(y))
    xl =  JuMP.dual(JuMP.LowerBoundRef(x))
    aggregated_sensitivity = +(yu, xu, yl, xl)
    lb = JuMP.objective_bound(model)
    return aggregated_sensitivity, lb
end

# 馃數 The following loop is deemed a "presolve" procedure
for i in 1:17 # The aim is to find a good artificial initial bound
    as, lb = get_result_with(exp(i))
    @info "i = $i, sensitivity = $as, lb = $lb"
end # we observe that `exp(15)` might be a good artificial bound

ib = exp(15) # an artificial initial bound
model = JuMP.Model(Gurobi.Optimizer); JuMP.set_silent(model)
JuMP.@variable(model, -ib <= y <= ib); JuMP.@variable(model, -ib <= x <= ib)
JuMP.@expression(model, X, LinearAlgebra.Symmetric([0 y; y x]))
JuMP.@objective(model, Min, 2 * y)
while true
    optimise(model)
    Xt = LinearAlgebra.Symmetric(JuMP.value.(X))
    _, t = my_eigmin(Xt)
    (t' * Xt * t <= -CoT ? JuMP.@constraint(model, t' * X * t >= 0) : break)
end
optimise(model)
yu = -JuMP.dual(JuMP.UpperBoundRef(y))
xu = -JuMP.dual(JuMP.UpperBoundRef(x))
yl =  JuMP.dual(JuMP.LowerBoundRef(y))
xl =  JuMP.dual(JuMP.LowerBoundRef(x))
aggregated_sensitivity = +(yu, xu, yl, xl)
@assert aggregated_sensitivity == 0
JuMP.delete_lower_bound(x);JuMP.delete_lower_bound(y);
JuMP.delete_upper_bound(x);JuMP.delete_upper_bound(y);
optimise(model) # At this time, the `model` is indeed a valid LP relaxation of the primal SDP
@assert LinearAlgebra.norm(JuMP.value.(X)) == 0 # we succeed in finding the optimal solution of the primal SDP
@assert JuMP.objective_bound(model) == JuMP.objective_value(model) == 0

It is nontrivial to decide a good initial artificial bound, no wonder Mosek and Clarabel both fails on this problem.