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—the free variables. Because of this, I have to add a presolve loop (
).
The full code is as follows (this is a very mature (≈definitive) 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.