When adding a SumOfSquares constraint, is it possible to get the JuMP variable corresponding to the Gram matrix?

Specific Domains Optimization (Mathematical)
1

Hello!

I have an optimization model with an SOS constraint and I would like to add extra constraints on the corresponding Gram matrix. Unfortunately, I could not find a way to obtain the JuMP variable corresponding to the Gram matrix (the function ‘gram_matrix’ can only be called after the model has been optimized).

Behind the hood, a constrain p in SOSCone() must be converted in p = Phi^T * G * Phi, with G≥0 (in the psd sense) and Phi being the polynomial basis. I would like to add extra constraints on the Gram matrix G. More precisely, I would like to constraint the condition number of G to be not greater than a user-defined κ. To this end, I want to add the constraint λI ≤ G ≤ λκ*I, where λ≥0 is an optimization variable.

The code would be something like this:

using JuMP
using SumOfSquares
using CSDP
using LinearAlgebra

model = SOSModel(CSDP.Optimizer)

@polyvar x
@variable(model, a)

p = x^2 + a

@objective(model, Min, a)
cons = @constraint(model, p in SOSCone())

κ = 1e4 # upper bound on the condition number of the Gram matrix

G = gram_matrix(cons) # this is throwing an error: OptimizeNotCalled()

@variable(model, λ≥0)
@constraint(model, G-λ*I in PSDCone())
@constraint(model, λ*κ*I - G in PSDCone())

optimize!(model)

Ideally, it would be possible to do the same thing when a domain is used in the SOS constraint, i.e.,

cons = @constraint(model, p in SOSCone(), domain = @set(x>=1))

In that case, there are several Gram matrices associated with that constraint.

Is there a way to do that?

PS: It is my first post here, I am happy to edit it if I am not following the guidelines properly :slightly_smiling_face:

1 Like
2

Hi @aaspeel, welcome to the forum :smile: you’ve followed the guidelines perfectly.

I don’t know the answer, but @blegat can help.

1 Like
3

At the moment, getting this matrix of variables is unfortunately not supported. For future reference, I have opened an issue to track this feature.
As a workaround, you can do

@variable(model, q, SOSPoly(monomials(x, 0:1)))
@constraint(model, p == q)
G = value_matrix(q)
1 Like