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