# SLOW\_PROGRESS or NUMERICAL\_ERROR: a simple but nontrivial SDP example, solved with LP relaxation

**URL:** https://discourse.julialang.org/t/slow-progress-or-numerical-error-a-simple-but-nontrivial-sdp-example-solved-with-lp-relaxation/127399
**Category:** Optimization (Mathematical)
**Tags:** question, examples, mosek, sdp
**Created:** [March 27, 2025, 4:27am UTC](https://discourse.julialang.org/t/slow-progress-or-numerical-error-a-simple-but-nontrivial-sdp-example-solved-with-lp-relaxation/127399 "2025-03-27T04:27:49Z")
**Posts on this page:** 1
**Showing post:** 6

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 28, 2025, 9:05am UTC](https://discourse.julialang.org/t/slow-progress-or-numerical-error-a-simple-but-nontrivial-sdp-example-solved-with-lp-relaxation/127399/6 "2025-03-28T09:05:31Z")

</div>

Due to their poor performance on this simple test case, I propose to solve it via LP relaxation.  
The crux, is to introduce a proper artificial bound.  
After adding about 23 linear cuts, we converge to the global optimal value `0.0`.

```julia
import JuMP, Gurobi, LinearAlgebra
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
ϵ = 1e-13 # tolerance on the smallest eigenvalue of a Symmetric matrix
ib = 9. # 🍏 pick an artificial initial bound
model = JuMP.Model(Gurobi.Optimizer) 
JuMP.@variable(model, -ib <= y <= ib); JuMP.@variable(model, -ib <= x <= ib)
JuMP.@expression(model, X, LinearAlgebra.Symmetric([0 y; y x]))
# JuMP.@constraint(model, X in PSDCone()) # We'll substitute this PSD cut by linear cuts
JuMP.@objective(model, Min, 2 * y)
for i in 0:typemax(Int) # iteratively tighten the LP relaxation
    optimise(model)
    lb = JuMP.objective_bound(model); @info "i = $i, lb = $lb"
    Xt = LinearAlgebra.Symmetric(JuMP.value.(X))
    valmin, t = my_eigmin(Xt)
    (valmin > -ϵ ? break : JuMP.@constraint(model, t' * X * t >= 0))
end
# during the above for-loop, there is only one (mild) Warning reported by Gurobi
# julia> [ Info: i = 22, lb = -4.365443091908161e-6
# julia> Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
optimise(model)
S = [
    -JuMP.dual(JuMP.UpperBoundRef(x)) -JuMP.dual(JuMP.UpperBoundRef(y)); 
    JuMP.dual(JuMP.LowerBoundRef(x)) JuMP.dual(JuMP.LowerBoundRef(y))
] # We observe that sensitivity == 0
# Hence we can delete our artificial bounds
JuMP.delete_upper_bound(x); JuMP.delete_upper_bound(y) 
JuMP.delete_lower_bound(x); JuMP.delete_lower_bound(y) 
optimise(model) # then we reoptimize
lb = JuMP.objective_bound(model); @info "lb = $lb"
Xt = LinearAlgebra.Symmetric(JuMP.value.(X))
valmin, t = my_eigmin(Xt)
@assert valmin >= -ϵ # double check the primal feasibility
ub = 2 * Xt[1, 2]
# ✅ here lb = ub = 0.0
# ✅ The LP relaxation method converges to global optimality

```

Remark: tuning `ib` is nontrivial, e.g. `ib = 7.8` would be an improper choice, where we cannot leave `for` quickly.

There are 2 to-be-improved issues: [smallest](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/7) eigenval and [Symmetric](https://discourse.julialang.org/t/retaining-the-property-of-symmetric/127437).

---

_[View the full topic](https://discourse.julialang.org/t/slow-progress-or-numerical-error-a-simple-but-nontrivial-sdp-example-solved-with-lp-relaxation/127399)._