A trivial SOS decomposition example

Contributed by: votroto Adapted from: Examples 3.25 of [BPT12]

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

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

\[ 1 + x^{2} - xy^{2} + y^{4} \]

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)

\[ (1) + (1)x^{2} + (-1)xy^{2} + (1)y^{4} \text{ is SOS} \]

Thirdly, optimize the feasibility problem!

optimize!(model)

Success: SDP solved
Primal objective value: 0.0000000e+00
Dual objective value: 0.0000000e+00
Relative primal infeasibility: 4.15e-17
Relative dual infeasibility: 5.00e-11
Real Relative Gap: 0.00e+00
XZ Relative Gap: 3.63e-10
DIMACS error measures: 6.34e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 3.63e-10
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.7423566403018639 - 0.5916342145954684*x + 0.9492757372229985*y^2)^2 + (-6.60921131771784e-17 + 1.1201589623631873*y - 2.4872525425116695e-16*x - 1.0534655002111817e-16*y^2)^2 + (0.6232480104529261 - 0.782024243528594*x + 1.1102230246251565e-16*y^2)^2 + (-0.24590350966628946 - 1.2342687879752308e-17*y - 0.19597713808894593*x - 0.31444486753600015*y^2)^2

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

polynomial(sos_dec, Float32)

\[ 1.0 - 1.4199712 \cdot 10^{-16}y - 1.9428903 \cdot 10^{-16}x + 2.220446 \cdot 10^{-16}y^{2} - 5.5238587 \cdot 10^{-16}xy + x^{2} - 2.2824759 \cdot 10^{-16}y^{3} - xy^{2} + y^{4} \]

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 with row/column basis:
 MonomialBasis([1, y, x, y^2])
And entries in a 4×4 SymMatrix{Float64}:
  1.0                    0.0                 …  -0.6273780504812863
  0.0                    1.2547561009625725      0.0
  4.933553565678038e-17  0.0                    -0.5
 -0.6273780504812863     0.0                     1.0

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

\[ 1.0 + 9.867107131356075 \cdot 10^{-17}x + x^{2} - xy^{2} + y^{4} \]

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.