# How to get non-trivial solution with SDP in SCS?

I am trying to solve the following, but I want to find the nonzero solution (right now it returns only trivial solutions). Is there anyway to do that?

Basically I want to find matrix G (symmetric and PSD) and matrix D (symmetric and PSD) such that the block matrix M is PSD.

using JuMP
using SCS
using LinearAlgebra

model = Model(SCS.Optimizer)

s = 2
r = 2

A = [1/2  0;
4/3  1/2];
U = [1 -1/2;
1 -5/6];
B = [19/16 19/16;
1/4  3/4];
V = [1 -3/4;
0  0];

@variable(model, G[1:r, 1:r],PSD)
@variable(model, D[1:s, 1:s],PSD)

M = [D*A + A'*D - B'*G*B D*U-B'*G*V; U'*D-V'*G*B G-V'*G*V];
@constraint(model, Symmetric(M)>=0 ,PSDCone())

optimize!(model)
# solution_summary(model)
print(value(G))
print(value(D))


Your equations appear to be invariant under scaling D,G \to \alpha D, \alpha G for any \alpha>0.

So, if you are looking for a solution where D is PD and not just PSD, you could eliminate the latter solution by adding a PSD constraint D - I \succeq 0, for example, instead of D \succeq 0. Or, equivalently, change variables from D to D = X + I where X is PSD.

@stevengj Thanks. I tried your approach, but when I examine the eigenvalues of G I found that one of them is negative (so G is not PSD). Am I misunderstanding something? Here is the modification

using JuMP
using SCS
using LinearAlgebra

model = Model(SCS.Optimizer)
set_optimizer_attribute(model, "max_iters", 1000000)
set_optimizer_attribute(model, "eps_abs", 1e-16)

s = 2
r = 2

A = [1/2  0;
4/3  1/2];
U = [1 -1/2;
1 -5/6];
B = [19/16 19/16;
1/4  3/4];
V = [1 -3/4;
0  0];

@variable(model, G[1:r, 1:r])
@variable(model, D[1:s, 1:s])

M = [D*A + A'*D - B'*G*B D*U-B'*G*V; U'*D-V'*G*B G-V'*G*V];

@expression(model, GG, G - I(r))
@expression(model, DD, D - I(s))

@constraint(model, Symmetric(GG)>=0 ,PSDCone())
@constraint(model, Symmetric(DD)>=0 ,PSDCone())
@constraint(model, Symmetric(M)>=0 ,PSDCone())

optimize!(model)
G_opt = value.(G)
D_opt = value.(D)
M_opt = [D_opt*A + A'*D_opt - B'*G_opt*B D_opt*U-B'*G_opt*V; U'*D_opt-V'*G_opt*B G_opt-V'*G_opt*V]

println(G_opt)
println(minimum(eigvals(G_opt)))
println(D_opt)
println(minimum(eigvals(D_opt)))
println(M_opt)
println(minimum(eigvals(Symmetric(M_opt))))


The output is

[6.3955477328463965 -4.796660799634771; -7.0815796227759416 5.299305393113521]
-0.006492824312413603
[4.34668485201598 1.477633175045142; 0.5842394825405693 1.6992887331325797]
1.405745872320629
[1.2722728261505036 0.10538009221944522 1.1546319456101628e-14 0.963535115083558; 0.5686569509549435 0.278765718069089 9.769962616701378e-15 0.004452476207420064; -1.464623398289853 -0.8202954248512961 0.0 2.6645352591003757e-14; 2.136452138176338 0.8430224679720428 -2.284918823141144 1.7018097933874223]
1.0731483923287292e-16
------------------------------------------------------------------
SCS v3.2.0 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 8, constraints m: 16
cones: 	  s: psd vars: 16, ssize: 3
settings: eps_abs: 1.0e-16, eps_rel: 1.0e-04, eps_infeas: 1.0e-16
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 1000000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
nnz(A): 44, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
0| 1.00e+00  4.15e-02  9.29e-02  4.65e-02  1.00e-01  2.38e-03
175| 3.82e-14  5.30e-18  7.08e-20 -3.54e-20  2.46e-04  3.83e-03
------------------------------------------------------------------
status:  solved
timings: total: 5.99e-03s = setup: 2.14e-03s + solve: 3.84e-03s
lin-sys: 5.65e-05s, cones: 3.13e-03s, accel: 5.13e-05s
------------------------------------------------------------------
objective = -0.000000
------------------------------------------------------------------


Did it find a feasible solution? Otherwise it seems like your constraint G - I \succeq 0 is violated if G has any eigenvalue ≤ 1, so the optimizer must be failing.

(Note that it should be sufficient to constrain D - I \succeq 0 and G \succeq 0 to exclude the trivial solution. But if a solution exists where both D and G are PD, not just PSD, then you can use your formulation)

(I’m don’t really use JuMP myself, so maybe I’m mis-reading your syntax.)

Thanks a lot. I have included the output that I got which shows that the problem was solved to high precision. Perhaps someone else could correct anything that I did not do correctly.

I haven’t tried your problem, but I noticed once obvious thing:

Using Symmetric is a hint to tell JuMP that the matrix GG is symmetric. If the matrix is not symmetric, then you can get incorrect results because JuMP will not add constraints to ensure the matrix is symmetric.

Try removing Symmetric: @constraint(model, GG >= 0, PSDCone()).

From a conceptual point of view though, I think you need to reconsider what you’re trying to do. The solver will find you a solution that satisfies the constraints as you have written them. If you want a different solution, it means that you are missing some constraints. That might, for example, mean adding a constraint like @constraint(model, G[1, 1] == 1) to fix one of the variables to a known value, around which the rest of the solution can be built.

Thanks. But I think we should have Symmetric() because I want G to be symmetric right? I tried it without it, and the solver did not converge anyway.

I think you are right with regards to your second point. The problem is that I don’t know what the solution would look like so adding a constraint like this may make the problem infeasible.

If you want symmetric or PSD variables, use:

@variable(model, G[1:r, 1:r], Symmetric)
@variable(model, D[1:s, 1:s], PSD)


Documentation: Variables · JuMP

If I want G to be BOTH Symmetric and PSD, how do I do that? (It is not clear from the documentation how to do that)

Use PSD, which ensures that it is symmetric and positive semidefinite.

(It is not clear from the documentation how to do that)

I’ll update the documentation to clarify the point