Term sparsity

Adapted from: Example 3.5 of [WML20b]

[WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre. TSSOS: A Moment-SOS hierarchy that exploits term sparsity. arXiv preprint arXiv:1912.08899 (2020).

[WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre. Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension. arXiv preprint arXiv:2003.03210 (2020).

using DynamicPolynomials
@polyvar x[1:3]

(DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[x₁, x₂, x₃],)

We would like to find the minimum value of the polynomial

poly = x[1]^2 - 2x[1]*x[2] + 3x[2]^2 - 2x[1]^2*x[2] + 2x[1]^2*x[2]^2 - 2x[2]*x[3] + 6x[3]^2 + 18x[2]^2*x[3] - 54x[2]*x[3]^2 + 142x[2]^2*x[3]^2

\[ 6x_{3}^{2} - 2x_{2}x_{3} + 3x_{2}^{2} - 2x_{1}x_{2} + x_{1}^{2} - 54x_{2}x_{3}^{2} + 18x_{2}^{2}x_{3} - 2x_{1}^{2}x_{2} + 142x_{2}^{2}x_{3}^{2} + 2x_{1}^{2}x_{2}^{2} \]

The minimum value of the polynomial can be found to be zero.

import CSDP
solver = CSDP.Optimizer
using SumOfSquares
function sos_min(sparsity)
    model = Model(solver)
    @variable(model, t)
    @objective(model, Max, t)
    con_ref = @constraint(model, poly - t in SOSCone(), sparsity = sparsity)
    optimize!(model)
    return value(t), moment_matrix(con_ref)
end

bound, ν = sos_min(Sparsity.NoPattern())
bound

-3.0596397637694395e-10

We find the corresponding minimizer (0, 0, 0) by matching the moments of the moment matrix with a dirac measure centered at this minimizer.

atomic_measure(ν, 1e-6)

Atomic measure on the variables x[1], x[2], x[3] with 1 atoms:
 at [2.0746548419820583e-6, 1.2388378355883834e-6, 2.3096229027153988e-8] with weight 1.000000001500628

We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables.

ν.basis

MonomialBasis([1, x₃, x₂, x₁, x₂x₃, x₁x₂])

Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.

bound, ν = sos_min(Sparsity.Monomial())
bound

-3.0596397637694395e-10

Which is not suprising as no sparsity reduction could be performed.

[sub.basis for sub in ν.blocks]

1-element Vector{MonomialBasis{DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}:
 MonomialBasis([1, x₃, x₂, x₁, x₂x₃, x₁x₂])

Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0.

bound, ν = sos_min(Sparsity.Monomial(ChordalCompletion()))
bound

-0.0035511999768404467

However, this bound was obtained with an SDP with 4 matrices of size 3x3.

[sub.basis for sub in ν.blocks]

4-element Vector{MonomialBasis{DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}:
 MonomialBasis([x₂, x₁, x₁x₂])
 MonomialBasis([x₂, x₂x₃, x₁x₂])
 MonomialBasis([x₃, x₂, x₂x₃])
 MonomialBasis([1, x₂x₃, x₁x₂])

This page was generated using Literate.jl.