I have to solve an SDP program. For example:
min tr(Σ * A )
s.t. Σ[1:2,1:2] >= AΣA’
Σ >= 0
Here is what I have:
n = 4
d = 2
A = rand(d,n)
B = rand(n,n)
model = Model(ProxSDP.Optimizer)
@variable(model, Σ[1:n, 1:n], PSD)
@objective(model, Min, tr(Σ*B)
@constraint(model, Σ[1:d,1:d] .>= A*Σ*A')
JuMP.optimize!(model)
Σ_sol = JuMP.value.(Σ)
After optimizing this, I found that the optimal value returned is not positive semidefinite. What could be the problem here?
Did you check termination_status(model) and primal_status(model)?
It says that the termination status is optimal and the primal status is feasible_point
PS: I printed this after the optimize method, i.e:
JuMP.optimize!(model)
println(primal_status(model))
println(termination_status(model))
Do you get the same results with a different solver? I’ve noticed ProxSDP is not always precise, at least with the default settings.
By the way, Home · ConvexTests.jl shows the results of various solvers with various test problems, including SDPs, and what solver settings were used.
I just tried SCS.Optimizer, it also gives a result that is not positive semidefinite. what else would you suggest me to try?
Hm, I just tried it myself and am getting mostly infeasible statuses. I suppose it depends on the random choices. With
using JuMP, StableRNGs, ProxSDP, SCS, LinearAlgebra
while true
N = rand(1:1000)
for solver in (SCS, ProxSDP)
n = 4
d = 2
rng = StableRNG(N)
A = rand(rng, d,n)
B = rand(rng, n,n)
model = Model(solver.Optimizer)
set_silent(model)
@variable(model, Σ[1:n, 1:n], PSD)
@objective(model, Min, tr(Σ*B)) # fixed paren
@constraint(model, Σ[1:d,1:d] .>= A*Σ*A')
JuMP.optimize!(model)
Σ_sol = JuMP.value.(Σ)
if termination_status(model) == JuMP.MOI.OPTIMAL
@show N
@show (solver, eigvals(Σ_sol), termination_status(model), primal_status(model))
end
end
end
I could find some feasible points, but the solutions were often zero, e.g.
N = 202
(solver, eigvals(Σ_sol), termination_status(model), primal_status(model)) = (SCS, [-7.122216469209867e-14, -3.308630615603461e-14, 2.6572571167405607e-14, 3.266987164674012e-14], MathOptInterface.OPTIMAL, MathOptInterface.FEASIBLE_POINT)
N = 202
(solver, eigvals(Σ_sol), termination_status(model), primal_status(model)) = (ProxSDP, [0.0, 0.0, 0.0, 0.0], MathOptInterface.OPTIMAL, MathOptInterface.FEASIBLE_POINT)
I was actually doing it in another SDP program, the one I posted is just a minimal example. But when I was using ProxSDP, do you know why the solution is not positive semidefinite even if the terminal condition is optimal?
how non-PSD is it? e.g. what is eigmin(Σ_sol)? If it’s only a little bit negative (like -1e-8) then that’s pretty much to be expected as numerical error. If it’s very negative, then that could be a problem with the solver or some kind numerical instability (though it should not give you OPTIMAL in that case).
I would try other solvers like COSMO and Hypatia as well.
It means a lot to know how much non-PSD the solution is.
ProxSDP can be less efficient in finding a solution, but it should tell the user if the solution is feasible (within tolerances) correctly, if it does not, then its a bug.
On the non-PSD thing, ProxSDP tends to be really good at enforcing that, so ProxSDP solutions tend to be more SDP feasible and linear feasible, as opposed to SCS/COSMO/SDPNAL.
Case studies in the paper highlighted this feature: https://www.tandfonline.com/doi/abs/10.1080/02331934.2020.1823387
I would recommend testing Mosek, which is very stable.
Hypatia also have a good chance of returning good quality solutions since it is second order.
I just tried Pkg.add("Hypatia"), but it wouldn’t work. Do you know how I can import the package?
You might be having a dependency conflict, you should probably try adding it in a new pkg environment. (see 4. Working with Environments · Pkg.jl)
also, Hypatia requires Julia 1.5, so if you’re on an older version you would need to upgrade. But it’s easier to help diagnose the issue if you include the error message 