# Positve semidefinite constraint on SDP program in JuMP
**URL:** https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972
**Category:** Optimization (Mathematical)
**Created:** [January 26, 2021, 11:56am UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972 "2021-01-26T11:56:58Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@cheng\_chen](https://discourse.julialang.org/u/cheng_chen)
#### Post date: [January 26, 2021, 11:56am UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/1 "2021-01-26T11:56:58Z")
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:
```julia
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?
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [January 26, 2021, 12:18pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/2 "2021-01-26T12:18:57Z")
Did you check `termination_status(model)` and `primal_status(model)`?
---
### Author:  [@cheng\_chen](https://discourse.julialang.org/u/cheng_chen)
#### Post date: [January 26, 2021, 12:24pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/3 "2021-01-26T12:24:57Z")
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:
```julia
JuMP.optimize!(model)
println(primal_status(model))
println(termination_status(model))
```
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [January 26, 2021, 12:32pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/4 "2021-01-26T12:32:29Z")
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](https://ericphanson.github.io/ConvexTests.jl/dev/) shows the results of various solvers with various test problems, including SDPs, and what solver settings were used.
---
### Author:  [@cheng\_chen](https://discourse.julialang.org/u/cheng_chen)
#### Post date: [January 26, 2021, 12:35pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/5 "2021-01-26T12:35:33Z")
I just tried SCS.Optimizer, it also gives a result that is not positive semidefinite. what else would you suggest me to try?
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [January 26, 2021, 12:54pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/6 "2021-01-26T12:54:24Z")
Hm, I just tried it myself and am getting mostly infeasible statuses. I suppose it depends on the random choices. With
```julia
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.
```julia
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)
```
---
### Author:  [@cheng\_chen](https://discourse.julialang.org/u/cheng_chen)
#### Post date: [January 26, 2021, 1:09pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/7 "2021-01-26T13:09:11Z")
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?
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [January 26, 2021, 1:29pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/8 "2021-01-26T13:29:39Z")
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.
---
### Author:  [@joaquimg](https://discourse.julialang.org/u/joaquimg)
#### Post date: [January 26, 2021, 1:48pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/9 "2021-01-26T13:48:57Z")
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](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.
---
### Author:  [@cheng\_chen](https://discourse.julialang.org/u/cheng_chen)
#### Post date: [January 27, 2021, 3:37am UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/10 "2021-01-27T03:37:28Z")
I just tried `Pkg.add("Hypatia")`, but it wouldn’t work. Do you know how I can import the package?
---
### Author:  [@joaquimg](https://discourse.julialang.org/u/joaquimg)
#### Post date: [January 27, 2021, 1:46pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/12 "2021-01-27T13:46:59Z")
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](https://julialang.github.io/Pkg.jl/v1/environments/))
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [January 27, 2021, 2:51pm UTC](https://discourse.julialang.org/t/positve-semidefinite-constraint-on-sdp-program-in-jump/53972/13 "2021-01-27T14:51:56Z")
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 🙂