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.PolyVar{true}[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
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
Iter: 7 Ap: 9.00e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.1796787e-13 Success: SDP solved Primal objective value: 0.0000000e+00 Dual objective value: 0.0000000e+00 Relative primal infeasibility: 6.80e-09 Relative dual infeasibility: 2.83e-50 Real Relative Gap: 0.00e+00 XZ Relative Gap: 8.00e-50 DIMACS error measures: 1.02e-08 0.00e+00 2.83e-50 0.00e+00 0.00e+00 8.00e-50 CSDP 6.2.0 Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 9.71e-01 Pobj: -2.4909930e+02 Ad: 8.28e-01 Dobj: 7.2430886e+02 Iter: 2 Ap: 9.99e-01 Pobj: -1.6816440e+02 Ad: 9.32e-01 Dobj: 1.6018610e+02 Iter: 3 Ap: 9.96e-01 Pobj: -6.7498290e+00 Ad: 9.24e-01 Dobj: 6.1399335e+01 Iter: 4 Ap: 1.00e+00 Pobj: -1.4722921e+00 Ad: 9.01e-01 Dobj: 8.9329470e+00 Iter: 5 Ap: 1.00e+00 Pobj: -5.0478524e-01 Ad: 7.90e-01 Dobj: 2.7911688e+00 Iter: 6 Ap: 1.00e+00 Pobj: -1.5409216e-01 Ad: 7.76e-01 Dobj: 8.7431136e-01 Iter: 7 Ap: 1.00e+00 Pobj: -4.8559418e-02 Ad: 8.12e-01 Dobj: 2.4102522e-01 Iter: 8 Ap: 1.00e+00 Pobj: -1.2896403e-02 Ad: 8.39e-01 Dobj: 6.6196298e-02 Iter: 9 Ap: 1.00e+00 Pobj: -1.7033049e-03 Ad: 1.00e+00 Dobj: 1.0847092e-02 Iter: 10 Ap: 1.00e+00 Pobj: -2.4388656e-04 Ad: 1.00e+00 Dobj: 7.6100068e-04 Iter: 11 Ap: 9.96e-01 Pobj: -7.3346824e-06 Ad: 1.00e+00 Dobj: -4.5795945e-05 Iter: 12 Ap: 9.97e-01 Pobj: -2.2313479e-07 Ad: 1.00e+00 Dobj: -7.5985091e-05 Iter: 13 Ap: 1.00e+00 Pobj: -9.6177081e-09 Ad: 1.00e+00 Dobj: -4.9851670e-06
-3.045195207107554e-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.
extractatoms(ν, 1e-6)
Atomic measure on the variables x[1], x[2], x[3] with 1 atoms: at [2.0703595321702734e-6, 1.2367907094787687e-6, 2.2362397702348983e-8] with weight 1.0000000014903554
We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables.
ν.basis
MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₂x₃, x₁, x₂, x₃, 1])
Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.
bound, ν = sos_min(Sparsity.Monomial())
bound
Iter: 14 Ap: 9.68e-01 Pobj: -3.0451952e-10 Ad: 9.70e-01 Dobj: -2.0275242e-07 Success: SDP solved Primal objective value: -3.0451952e-10 Dual objective value: -2.0275242e-07 Relative primal infeasibility: 1.62e-13 Relative dual infeasibility: 1.55e-09 Real Relative Gap: -2.02e-07 XZ Relative Gap: 2.86e-09 DIMACS error measures: 1.75e-13 0.00e+00 2.56e-09 0.00e+00 -2.02e-07 2.86e-09 CSDP 6.2.0 Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 9.71e-01 Pobj: -2.4909930e+02 Ad: 8.28e-01 Dobj: 7.2430886e+02 Iter: 2 Ap: 9.99e-01 Pobj: -1.6816440e+02 Ad: 9.32e-01 Dobj: 1.6018610e+02 Iter: 3 Ap: 9.96e-01 Pobj: -6.7498290e+00 Ad: 9.24e-01 Dobj: 6.1399335e+01 Iter: 4 Ap: 1.00e+00 Pobj: -1.4722921e+00 Ad: 9.01e-01 Dobj: 8.9329470e+00 Iter: 5 Ap: 1.00e+00 Pobj: -5.0478524e-01 Ad: 7.90e-01 Dobj: 2.7911688e+00 Iter: 6 Ap: 1.00e+00 Pobj: -1.5409216e-01 Ad: 7.76e-01 Dobj: 8.7431136e-01 Iter: 7 Ap: 1.00e+00 Pobj: -4.8559418e-02 Ad: 8.12e-01 Dobj: 2.4102522e-01 Iter: 8 Ap: 1.00e+00 Pobj: -1.2896403e-02 Ad: 8.39e-01 Dobj: 6.6196298e-02 Iter: 9 Ap: 1.00e+00 Pobj: -1.7033049e-03 Ad: 1.00e+00 Dobj: 1.0847092e-02 Iter: 10 Ap: 1.00e+00 Pobj: -2.4388656e-04 Ad: 1.00e+00 Dobj: 7.6100068e-04 Iter: 11 Ap: 9.96e-01 Pobj: -7.3346824e-06 Ad: 1.00e+00 Dobj: -4.5795945e-05 Iter: 12 Ap: 9.97e-01 Pobj: -2.2313479e-07 Ad: 1.00e+00 Dobj: -7.5985091e-05 Iter: 13 Ap: 1.00e+00 Pobj: -9.6177081e-09 Ad: 1.00e+00 Dobj: -4.9851670e-06
-3.045195207107554e-10
Which is not suprising as no sparsity reduction could be performed.
[sub.basis for sub in ν.sub_moment_matrices]
1-element Vector{MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}}: MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₂x₃, x₁, x₂, x₃, 1])
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
Iter: 14 Ap: 9.68e-01 Pobj: -3.0451952e-10 Ad: 9.70e-01 Dobj: -2.0275242e-07 Success: SDP solved Primal objective value: -3.0451952e-10 Dual objective value: -2.0275242e-07 Relative primal infeasibility: 1.62e-13 Relative dual infeasibility: 1.55e-09 Real Relative Gap: -2.02e-07 XZ Relative Gap: 2.86e-09 DIMACS error measures: 1.75e-13 0.00e+00 2.56e-09 0.00e+00 -2.02e-07 2.86e-09 CSDP 6.2.0 Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 9.17e-01 Pobj: -3.5364812e+02 Ad: 6.71e-01 Dobj: 5.8532180e+02 Iter: 2 Ap: 8.26e-01 Pobj: -3.6133495e+02 Ad: 9.36e-01 Dobj: 1.9736716e+02 Iter: 3 Ap: 9.70e-01 Pobj: -1.1800069e+02 Ad: 8.23e-01 Dobj: 1.4747690e+02 Iter: 4 Ap: 1.00e+00 Pobj: -1.2898462e+01 Ad: 8.32e-01 Dobj: 5.5286769e+01 Iter: 5 Ap: 1.00e+00 Pobj: -2.4900158e+00 Ad: 9.07e-01 Dobj: 8.7767866e+00 Iter: 6 Ap: 1.00e+00 Pobj: -1.2196157e+00 Ad: 8.79e-01 Dobj: 2.3717906e+00 Iter: 7 Ap: 1.00e+00 Pobj: -4.3509371e-01 Ad: 6.82e-01 Dobj: 1.1617211e+00 Iter: 8 Ap: 1.00e+00 Pobj: -1.7554721e-01 Ad: 7.67e-01 Dobj: 3.9978327e-01 Iter: 9 Ap: 1.00e+00 Pobj: -6.1256598e-02 Ad: 8.61e-01 Dobj: 9.4968927e-02 Iter: 10 Ap: 1.00e+00 Pobj: -2.0602943e-02 Ad: 1.00e+00 Dobj: 8.3003823e-03 Iter: 11 Ap: 1.00e+00 Pobj: -8.1086444e-03 Ad: 1.00e+00 Dobj: 2.1149126e-04 Iter: 12 Ap: 1.00e+00 Pobj: -5.0556265e-03 Ad: 1.00e+00 Dobj: -2.8106237e-03 Iter: 13 Ap: 1.00e+00 Pobj: -3.7777679e-03 Ad: 1.00e+00 Dobj: -3.4791027e-03 Iter: 14 Ap: 1.00e+00 Pobj: -3.5703395e-03 Ad: 1.00e+00 Dobj: -3.5995200e-03 Iter: 15 Ap: 1.00e+00 Pobj: -3.5520509e-03 Ad: 1.00e+00 Dobj: -3.6089795e-03 Iter: 16 Ap: 1.00e+00 Pobj: -3.5512394e-03 Ad: 1.00e+00 Dobj: -3.5545003e-03
-0.0035512009904666852
However, this bound was obtained with an SDP with 4 matrices of size 3x3.
[sub.basis for sub in ν.sub_moment_matrices]
4-element Vector{MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}}: MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₂x₃, x₂, x₃]) MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₂x₃, x₂, 1]) MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁, x₂, 1]) MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₁, 1])
This notebook was generated using Literate.jl.