# Sum-Of-Squares SOC constraint in JuMP

I’m trying to add a constraint of type `(x'Ax <= b)` in JuMP to be interpreted as SOC by the solver.
A is Positive definite.
Seems to be the same error as presented in this issue, but `A+A'` has positive eigenvalues:

``````eigvals(cov(d).data + cov(d).data')
10-element Array{Float64,1}:
2.0014222400461272
2.0605866115342453
2.2799706986926673
2.5317489320096396
2.7815912421610607
3.459123499910021
4.181503965703311
4.970808766218942
6.5318231430488245
53.903463261426204
``````

`A` also has positive eignenvalues:

``````eigvals(cov(d).data)
10-element Array{Float64,1}:
1.0007111200230636
1.0302933057671226
1.1399853493463337
1.2658744660048198
1.3907956210805303
1.7295617499550104
2.0907519828516556
2.485404383109471
3.2659115715244122
26.951731630713102
``````

From stacktrace, I can see that this happens at the call of `cholesky` (ref) and seems that the “error” was already reported: https://github.com/JuliaLang/julia/issues/30831

Could the solution be to call `cholesky(Hermitian(matrix))` or the same with `Symmetric ` in MOI?

Minimal working case:

``````using Random
using JuMP
using ECOS

n = 2
rng = MersenneTwister(1)
means = randn(rng, n)
sqrt_cov = rand(rng, n, n) ./ 2.0
Σ = sqrt_cov' * sqrt_cov

model = Model(ECOS.Optimizer)
@variable(model, x[i=1:n] >= 0);
@variable(model, y[i=1:n] <= 0);
v = x + y
@objective(model, Min, 0);
@constraint(model, v'Σ*v / 321.0^2 <= 2);

optimize!(model)
``````

Curious enough, the following doesn’t error:

``````using Random
using JuMP
using ECOS

n = 2
rng = MersenneTwister(1)
means = randn(rng, n)
sqrt_cov = rand(rng, n, n) ./ 2.0
Σ = sqrt_cov' * sqrt_cov

model = Model(ECOS.Optimizer)
@variable(model, x[i=1:n] >= 0);
@objective(model, Min, 0);
@constraint(model, x'Σ*x / 321.0^2 <= 2);

optimize!(model)
``````

and even more curious, the following doesn’t error:

``````using Random
using JuMP
using ECOS

n = 2
rng = MersenneTwister(1)
means = randn(rng, n)
sqrt_cov = rand(rng, n, n) ./ 2.0
Σ = sqrt_cov' * sqrt_cov

model = Model(ECOS.Optimizer)
@variable(model, x[i=1:n] >= 0);
@variable(model, y[i=1:n] <= 0);
@variable(model, v[i=1:n])
@constraint(model, v .== x + y)
@objective(model, Min, 0);
@constraint(model, v'Σ*v / 321.0^2 <= 2);

optimize!(model)

``````

Just use the `SecondOrderCone` syntax:

``````@constraint(model, [t; x] in SecondOrderCone())
``````
1 Like

yeah I did that, I was just curious why this didn’t work. However, Joaquim just showed me that indeed when `v` is not a variable, the constraint uses a matrix that is Positive-Semidefinite instead of Positive-Definite.

``````v'Σ*v = (x + y)'Σ*(x + y) = x'Σ*(x + y) + y'Σ*(x + y) = x'Σ*x + x'Σ*y + y'Σ*x + y'Σ*y
Since x'Σ*y == y'Σ*x, then:
v'Σ*v = (x + y)'Σ*(x + y) = x'Σ*(x + y) + y'Σ*(x + y) = x'Σ*x + x'Σ*y + y'Σ*x + y'Σ*y = x'Σ*x + 2x'Σ*y + y'Σ*y
And:
x'Σ*x + 2x'Σ*y + y'Σ*y == [x; y]' [Σ Σ;Σ Σ] [x; y]
At last:
[Σ Σ;Σ Σ] >= 0 and not > 0
``````