How can I forcefully bridge sum-of-squares constraints into PSD constraints?

Hi everyone,

I have model written using SumOfSquares.jl, it contains “F in S” pairs such as
Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Float64, Polynomial{true, Float64}, FullSpace}, Monomial{true}, MonomialVector{true}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, Tuple{}}, SOSCone, MonomialBasis}}

I’m guessing that before calling the solver (when calling optimize!() ?), the constraints get transformed into some regular PSD constraints like MOI.PositiveSemidefiniteConeTriangle

How can I do/force the transformation manually?

The reason I’m asking is that I want to work with something closer to the PSD representation of the model, to be able to use packages like
https://github.com/matbesancon/MathOptSetDistances.jl (which work only with the basic constraint types defined by MOI)

1 Like

It’s getting a little into the depths of JuMP and MOI, but you could do something like this:

using DynamicPolynomials
using SumOfSquares

@polyvar x1 x2
f1 = x1 + x2 + 1
f2 = 19 - 14x1 + 3x1^2 - 14x2 + 6x1*x2 + 3x2^2
f3 = 2x1 - 3x2
f4 = 18 - 32x1 + 12x1^2 + 48x2 - 36x1*x2 + 27x2^2
f = (1 + f1^2 * f2) * (30 + f3^2 * f4)

MOI.Utilities.@model(
    MyNewModel,
    (),
    (),
    (MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle),
    (),
    (),
    (),
    (),
    (MOI.VectorAffineFunction,),
    true,
)

model = SOSModel(MyNewModel{Float64})
@variable(model, γ)
@objective(model, Max, γ)
@constraint(model, f >= γ)
MOI.Utilities.attach_optimizer(model)
moi_model = backend(model).optimizer.model

julia> MOI.get(moi_model, MOI.ListOfConstraints())
2-element Vector{Tuple{DataType, DataType}}:
 (MathOptInterface.VectorAffineFunction{Float64}, MathOptInterface.Nonnegatives)
 (MathOptInterface.VectorAffineFunction{Float64}, MathOptInterface.PositiveSemidefiniteConeTriangle)

(p.s. The next version of JuMP will make this much easier by exposing an unsafe_backend function to avoid the .optimizer.model mess.)

1 Like

@blegat pointed out you could also use use

SOSModel(MOI.Utilities.Model{Float64})

to avoid having to create MyNewModel.

1 Like

Thanks for your reply,
your code works on my problem, but I am struggling with solving the model now.

Normally I would create a model by calling model = SOSModel(Mosek.Optimizer),
add some constraints, and call optimize!(model).
If I understand correctly from MOI API reference, MyNewModel now acts instead of the Mosek.Optimizer, but doesn’t implement the optimize!() method, obviously.

I’d like to get primal/dual variables for the problem specified by moi_model,
for which I need to solve it in that particular form (unless there is some other way).
How can I do that?

I solved it by copying to a new model.

using MosekTools
m = Model(Mosek.Optimizer)
MOI.copy_to(m.moi_backend, moi_model)
optimize!(m)

# get the value of γ by using its `MOI.VariableIndex`
x = MOI.get(m.moi_backend, MOI.VariablePrimal(), index(γ))