Adapted from: Example 4 of [L09]
[L09] Lofberg, Johan. Pre-and post-processing sum-of-squares programs in practice. IEEE transactions on automatic control 54, no. 5 (2009): 1007-1011.
using DynamicPolynomials
@polyvar x[1:3]
(DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[x₁, x₂, x₃],)
We would like to determine whether the following polynomial is a sum-of-squares.
poly = 1 + x[1]^4 + x[1] * x[2] + x[2]^4 + x[3]^2
In order to do this, we can solve the following Sum-of-Squares program.
import CSDP
solver = CSDP.Optimizer
using SumOfSquares
function sos_check(sparsity)
model = Model(solver)
con_ref = @constraint(model, poly in SOSCone(), sparsity = sparsity)
optimize!(model)
println(solution_summary(model))
return gram_matrix(con_ref)
end
g = sos_check(Sparsity.NoPattern())
g.basis.monomials
Success: SDP solved Primal objective value: 0.0000000e+00 Dual objective value: 0.0000000e+00 Relative primal infeasibility: 5.23e-17 Relative dual infeasibility: 5.00e-11 Real Relative Gap: 0.00e+00 XZ Relative Gap: 1.00e-10 DIMACS error measures: 7.85e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.00e-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: 8.10e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.0878274e+01 Iter: 2 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 2.7225538e+01 * Solver : CSDP * Status Result count : 1 Termination status : OPTIMAL Message from the solver: "Problem solved to optimality." * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : 0.00000e+00 Dual objective value : 0.00000e+00 * Work counters Solve time (sec) : 4.61102e-04
7-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
1
x₃
x₂
x₁
x₂²
x₁x₂
x₁²
As detailed in the Example 4 of [L09], we can exploit the sign symmetry of the polynomial to decompose the large positive semidefinite matrix into smaller ones.
g = sos_check(Sparsity.SignSymmetry())
monos = [sub.basis.monomials for sub in g.blocks]
Success: SDP solved Primal objective value: 0.0000000e+00 Dual objective value: 0.0000000e+00 Relative primal infeasibility: 7.46e-17 Relative dual infeasibility: 5.00e-11 Real Relative Gap: 0.00e+00 XZ Relative Gap: 1.21e-10 DIMACS error measures: 1.21e-16 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.21e-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: 7.15e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.7634216e+01 Iter: 2 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.5035802e+01 * Solver : CSDP * Status Result count : 1 Termination status : OPTIMAL Message from the solver: "Problem solved to optimality." * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : 0.00000e+00 Dual objective value : 0.00000e+00 * Work counters Solve time (sec) : 6.58989e-04
3-element Vector{DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}:
[x₂, x₁]
[x₃]
[1, x₂², x₁x₂, x₁²]
This notebook was generated using Literate.jl.