Positve semidefinite constraint on SDP program in JuMP

Specific Domains Optimization (Mathematical)
1

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?

2

Did you check termination_status(model) and primal_status(model)?

3

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))
4

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.

5

I just tried SCS.Optimizer, it also gives a result that is not positive semidefinite. what else would you suggest me to try?

6

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

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?

8

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.

9

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.

10

I just tried Pkg.add("Hypatia"), but it wouldn’t work. Do you know how I can import the package?

12

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)

13

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: