A trivial SOS decomposition example

Contributed by: votroto

using DynamicPolynomials
using SumOfSquares
import CSDP

The polynomial p = x^2 - x*y^2 + y^4 + 1 is SOS. We can, for example, decompose it as p = 3/4*(x - y^2)^2 + 1/4*(x + y)^2 + 1, which clearly proves that p is SOS, and there are infinitely many other ways to decompose p into sums of squares.

We can use SumOfSquares.jl to find such decompositions.

First, setup the polynomial of interest.

@polyvar x y
p = x^2 - x*y^2 + y^4 + 1

Secondly, constrain the polynomial to be nonnegative. SumOfSquares.jl transparently reinterprets polyonmial nonnegativity as the appropriate SOS certificate for polynomials nonnegative on semialgebraic sets.

model = SOSModel(CSDP.Optimizer)
@constraint(model, cref, p >= 0)

cref : $ (1)y^{4} + (-1)xy^{2} + (1)x^{2} + (1) \text{ is SOS} $

Thirdly, optimize the feasibility problem!

optimize!(model)

Iter: 17 Ap: 9.60e-01 Pobj: -3.5512010e-03 Ad: 9.60e-01 Dobj: -3.5513625e-03
Success: SDP solved
Primal objective value: -3.5512010e-03
Dual objective value: -3.5513625e-03
Relative primal infeasibility: 7.16e-14
Relative dual infeasibility: 1.65e-09
Real Relative Gap: -1.60e-07
XZ Relative Gap: 2.91e-09
DIMACS error measures: 7.72e-14 0.00e+00 4.44e-09 0.00e+00 -1.60e-07 2.91e-09
CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.7659722e+01
Iter:  2 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  1.2562719e+01

Lastly, recover a SOS decomposition. In general, SOS decompositions are not unique!

sos_dec = sos_decomposition(cref, 1e-4)

(-0.949275737222999*y^2 + 0.5916342145954685*x + 0.7423566403018638)^2 + (1.1201589623631871*y)^2 + (1.9762107540267777e-16*y^2 - 0.7820242435285937*x + 0.6232480104529262)^2 + (-0.31444486753600004*y^2 - 0.19597713808894598*x - 0.24590350966628918)^2

Converting, rounding, and simplifying - Huzza, Back where we began!

polynomial(sos_dec, Float32)

A deeper explanation and the unexplained 1e-4 parameter

p = x^2 - x*y^2 + y^4 + 1 can be represented in terms of its Gram matrix as

gram = gram_matrix(cref)

GramMatrix{Float64,MonomialBasis{DynamicPolynomials.Monomial{true},DynamicPolynomials.MonomialVector{true}},Float64,SymMatrix{Float64}}([0.9999999999999999 -0.5 0.0 -0.627378050481286; -0.5 1.0 0.0 -2.7061686225238353e-18; 0.0 0.0 1.254756100962572 0.0; -0.627378050481286 -2.7061686225238353e-18 0.0 1.0], MonomialBasis{DynamicPolynomials.Monomial{true},DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[y², x, y, 1]))

gram.basis.monomials' * gram.Q * gram.basis.monomials

where the matrix gram.Q is positive semidefinite, because p is SOS. If we could only get the decomposition gram.Q = V' * V, the SOS decomposition would simply be ||V * monomials||^2.

Unfortunately, we can not use Cholesky decomposition, since gram Q is only semidefinite, not definite. Hence, SumOfSquares.jl uses SVD decomposition instead and discards small singular values (in our case 1e-4).

This page was generated using Literate.jl.