Positve semidefinite constraint on SDP program in JuMP

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')
Σ_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:


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)
        @variable(model, Σ[1:n, 1:n], PSD)
        @objective(model, Min, tr(Σ*B)) # fixed paren
        @constraint(model, Σ[1:d,1:d] .>= A*Σ*A')
        Σ_sol = JuMP.value.(Σ)
        if termination_status(model) == JuMP.MOI.OPTIMAL
            @show N
            @show (solver, eigvals(Σ_sol), termination_status(model), primal_status(model))

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)
1 Like

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)

1 Like

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 :slight_smile:

1 Like