SumOfSquares Code only works for single constraint

Dear All,
I am trying out SOS decomposition of a polynomial matrix on a bounded (rectangular) domain. I generate an SOS matrix as the product P’*P and then multiply that by (1 - x[i]^2) for each of the variables that I am using. If I multiply by a single term, the code works correctly, (returns feasible). However if I multiply P’P by (1 - x[1]^2)(1-x[2]^2), the code returns infeasible! I do not understand why this happens. In the code changing NVarTst (line 10) from 1 to 2 causes the transition from feasible to infeasible.

using DynamicPolynomials
using SumOfSquares
using LinearAlgebra
using Plots
import CSDP
"""
   SOSPSDTst -- script to test SOSPSD algorithm
"""
NVar = 2 #number variables
NVarTst = 2 #number used in test ing for SOS matrix
NDeg = 1 # degree/2 of polynomials in the test matrix
NSize = 2 # size of the test matrix
xL = 1.0 #domain of test -xL..xL

solver = optimizer_with_attributes(CSDP.Optimizer)

@polyvar x[1:NVar]

#There should be a more elegant method to do this!
if NVar == 1
   S = @set xL^2 - x[1]^1 >= 0.0
end
if NVar == 2
   S = @set xL^2 - x[1]^2 >= 0 && xL^2 - x[2]^2 >= 0.0
end
if NVar == 3
   S = @set xL^2 - x[1]^2 >= 0.0 && xL^2 - x[2]^2 >= 0.0 && xL^2 - x[3]^2 >= 0.0
end
model = SOSModel(solver)

#fill in the polynomials in the test matrix
xmons = monomials(x, 0:NDeg)
P = Matrix{DynamicPolynomials.Polynomial}(undef, NSize, NSize)
for i = 1:NSize
   for j = 1:NSize
       local coefs = randn(Float64, length(xmons))
       P[i, j] = coefs'*xmons
   end
end

Q1= P'*P
Q = (Q1 + Q1')/2.0 # so that julia recognizes Q as symmetric
for i = 1:NVarTst
   global Q = Q .* (xL^2 - x[i]^2)
end

constr = @constraint(model, Q in PSDCone(), domain = S)

println(model)
optimize!(model)
println("\nResult is: \n")
print(solution_summary(model))
println("\n---------------\n")

@polyvar y[1:NSize]
model2 = SOSModel(solver)
z = y'*Q*y
constr2 = @constraint(model2, z >= 0.0, domain = S)
println(model2)
optimize!(model2)
println("\nResult is: \n")
println(solution_summary(model2))
println("\n-----------------------\n")

This is a question for @blegat.

The SumOfSquares condition is only sufficient for nonnegativity. When creating a model with SOSModel instead of Model, you are asking that the nonnegativity constraint z >= 0 is interpreted as z in SOSCone().
Now, you can reduce the conservativeness by increasing the degree of the certificate, at the expense of increased computational cost.
To do this, use the maxdegree keyword of @constraint. Here is an example: Nonconvex QP · SumOfSquares

1 Like

Thank you for your prompt answer. I tried your suggestion, to maxdegree = 16, but it did not help. Is there something wrong with my example?

What’s the termination_status of the solver, does it finish with OPTIMAL ?

1 Like

Thank you again for your reply. I think that the problem in JuMP is solved on the basis of Putinar’s Positivestellsatz. (Is this correct?) The theorem states that there exist SOS matrices…, but without specifying the degrees of these matrices. Is it possible that degree 16 is insufficient? In any case here is the output, with maxdegree = 16
Feasibility
Subject to
(1.3995440949783284)##243² + (0.5869250663560481)##241##243 + (0.10471918050754031)##241² + (0.6567392441399948)x₂##243² + (-0.3848851434790659)x₂##241##243 + (-0.12170245815222104)x₂##241² + (-3.20460361818109)x₁##243² + (-3.7699288635098194)x₁##241##243 + (-0.9189556297418807)x₁##241² + (-1.283839095005359)x₂²##243² + (-0.7239699758839327)x₂²##241##243 + (-0.0393391038321398)x₂²##241² + (-0.8520445217709228)x₁x₂##243² + (-0.2921718700158021)x₁x₂##241##243 + (0.6290662093637827)x₁x₂##241² + (0.4997626588825679)x₁²##243² + (3.196015550725679)x₁²##241##243 + (1.9866079395131657)x₁²##241² + (-0.6567392441399948)x₂³##243² + (0.3848851434790659)x₂³##241##243 + (0.12170245815222104)x₂³##241² + (3.20460361818109)x₁x₂²##243² + (3.7699288635098194)x₁x₂²##241##243 + (0.9189556297418807)x₁x₂²##241² + (-0.6567392441399948)x₁²x₂##243² + (0.3848851434790659)x₁²x₂##241##243 + (0.12170245815222104)x₁²x₂##241² + (3.20460361818109)x₁³##243² + (3.7699288635098194)x₁³##241##243 + (0.9189556297418807)x₁³##241² + (-0.11570499997296935)x₂⁴##243² + (0.1370449095278846)x₂⁴##241##243 + (-0.06538007667540051)x₂⁴##241² + (0.8520445217709228)x₁x₂³##243² + (0.2921718700158021)x₁x₂³##241##243 + (-0.6290662093637827)x₁x₂³##241² + (-0.6154676588555372)x₁²x₂²##243² + (-3.0589706411977944)x₁²x₂²##241##243 + (-2.051988016188566)x₁²x₂²##241² + (0.8520445217709228)x₁³x₂##243² + (0.2921718700158021)x₁³x₂##241##243 + (-0.6290662093637827)x₁³x₂##241² + (-1.8993067538608963)x₁⁴##243² + (-3.782940617081727)x₁⁴##241##243 + (-2.091327120020706)x₁⁴##241² + (0.6567392441399948)x₁²x₂³##243² + (-0.3848851434790659)x₁²x₂³##241##243 + (-0.12170245815222104)x₁²x₂³##241² + (-3.20460361818109)x₁³x₂²##243² + (-3.7699288635098194)x₁³x₂²##241##243 + (-0.9189556297418807)x₁³x₂²##241² + (0.11570499997296935)x₁²x₂⁴##243² + (-0.1370449095278846)x₁²x₂⁴##241##243 + (0.06538007667540051)x₁²x₂⁴##241² + (-0.8520445217709228)x₁³x₂³##243² + (-0.2921718700158021)x₁³x₂³##241##243 + (0.6290662093637827)x₁³x₂³##241² + (1.8993067538608963)x₁⁴x₂²##243² + (3.782940617081727)x₁⁴x₂²##241##243 + (2.091327120020706)x₁⁴x₂²##241² is SOS

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: 4.58e-01 Pobj: 0.0000000e+00 Ad: 8.42e-01 Dobj: 1.7570511e+01
Iter: 2 Ap: 4.71e-01 Pobj: 0.0000000e+00 Ad: 8.00e-01 Dobj: 3.1316099e+01
Iter: 3 Ap: 5.77e-01 Pobj: 0.0000000e+00 Ad: 6.81e-01 Dobj: 3.7915219e+01
Iter: 4 Ap: 6.29e-01 Pobj: 0.0000000e+00 Ad: 8.78e-01 Dobj: 4.4841314e+01
Iter: 5 Ap: 7.05e-01 Pobj: 0.0000000e+00 Ad: 8.19e-01 Dobj: 4.0867220e+01
Iter: 6 Ap: 4.49e-01 Pobj: 0.0000000e+00 Ad: 8.96e-01 Dobj: 3.3956470e+01
Iter: 7 Ap: 7.06e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 2.4975713e+01
Iter: 8 Ap: 7.70e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 1.1624175e+01
Iter: 9 Ap: 8.19e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.4118701e+00
Iter: 10 Ap: 6.94e-01 Pobj: 0.0000000e+00 Ad: 8.81e-01 Dobj: 1.5336094e+00
Iter: 11 Ap: 3.86e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.3982636e-01
Iter: 12 Ap: 7.22e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.1824011e-01
Iter: 13 Ap: 4.34e-01 Pobj: 0.0000000e+00 Ad: 5.75e-01 Dobj: 2.6649839e-01
Iter: 14 Ap: 1.93e-01 Pobj: 0.0000000e+00 Ad: 4.14e-01 Dobj: 8.8697271e-02
Iter: 15 Ap: 9.37e-02 Pobj: 0.0000000e+00 Ad: 1.25e-01 Dobj: -1.0731679e-02
Declaring primal infeasibility.
Success: SDP is primal infeasible
Certificate of primal infeasibility: a’*y=-1.00000e+00, ||A’(y)-Z||=1.65423e-10

Result is:

  • Solver : CSDP

  • Status
    Result count : 1
    Termination status : INFEASIBLE
    Message from the solver:
    “Problem is primal infeasible.”

  • Candidate solution (result #1)
    Primal status : INFEASIBLE_POINT
    Dual status : INFEASIBILITY_CERTIFICATE
    Objective value : 0.00000e+00
    Dual objective value : 1.00000e+00

  • Work counters
    Solve time (sec) : 2.33713e+02


Feasibility
Subject to
(1.3995440949783284)y₂² + (0.5869250663560481)y₁y₂ + (0.10471918050754031)y₁² + (0.6567392441399948)x₂y₂² + (-0.3848851434790659)x₂y₁y₂ + (-0.12170245815222104)x₂y₁² + (-3.20460361818109)x₁y₂² + (-3.7699288635098194)x₁y₁y₂ + (-0.9189556297418807)x₁y₁² + (-1.283839095005359)x₂²y₂² + (-0.7239699758839327)x₂²y₁y₂ + (-0.0393391038321398)x₂²y₁² + (-0.8520445217709228)x₁x₂y₂² + (-0.2921718700158021)x₁x₂y₁y₂ + (0.6290662093637827)x₁x₂y₁² + (0.4997626588825679)x₁²y₂² + (3.196015550725679)x₁²y₁y₂ + (1.9866079395131657)x₁²y₁² + (-0.6567392441399948)x₂³y₂² + (0.3848851434790659)x₂³y₁y₂ + (0.12170245815222104)x₂³y₁² + (3.20460361818109)x₁x₂²y₂² + (3.7699288635098194)x₁x₂²y₁y₂ + (0.9189556297418807)x₁x₂²y₁² + (-0.6567392441399948)x₁²x₂y₂² + (0.3848851434790659)x₁²x₂y₁y₂ + (0.12170245815222104)x₁²x₂y₁² + (3.20460361818109)x₁³y₂² + (3.7699288635098194)x₁³y₁y₂ + (0.9189556297418807)x₁³y₁² + (-0.11570499997296935)x₂⁴y₂² + (0.1370449095278846)x₂⁴y₁y₂ + (-0.06538007667540051)x₂⁴y₁² + (0.8520445217709228)x₁x₂³y₂² + (0.2921718700158021)x₁x₂³y₁y₂ + (-0.6290662093637827)x₁x₂³y₁² + (-0.6154676588555372)x₁²x₂²y₂² + (-3.0589706411977944)x₁²x₂²y₁y₂ + (-2.051988016188566)x₁²x₂²y₁² + (0.8520445217709228)x₁³x₂y₂² + (0.2921718700158021)x₁³x₂y₁y₂ + (-0.6290662093637827)x₁³x₂y₁² + (-1.8993067538608963)x₁⁴y₂² + (-3.782940617081727)x₁⁴y₁y₂ + (-2.091327120020706)x₁⁴y₁² + (0.6567392441399948)x₁²x₂³y₂² + (-0.3848851434790659)x₁²x₂³y₁y₂ + (-0.12170245815222104)x₁²x₂³y₁² + (-3.20460361818109)x₁³x₂²y₂² + (-3.7699288635098194)x₁³x₂²y₁y₂ + (-0.9189556297418807)x₁³x₂²y₁² + (0.11570499997296935)x₁²x₂⁴y₂² + (-0.1370449095278846)x₁²x₂⁴y₁y₂ + (0.06538007667540051)x₁²x₂⁴y₁² + (-0.8520445217709228)x₁³x₂³y₂² + (-0.2921718700158021)x₁³x₂³y₁y₂ + (0.6290662093637827)x₁³x₂³y₁² + (1.8993067538608963)x₁⁴x₂²y₂² + (3.782940617081727)x₁⁴x₂²y₁y₂ + (2.091327120020706)x₁⁴x₂²y₁² is SOS

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: 4.58e-01 Pobj: 0.0000000e+00 Ad: 8.42e-01 Dobj: 1.7570511e+01
Iter: 2 Ap: 4.71e-01 Pobj: 0.0000000e+00 Ad: 8.00e-01 Dobj: 3.1316099e+01
Iter: 3 Ap: 5.77e-01 Pobj: 0.0000000e+00 Ad: 6.81e-01 Dobj: 3.7915219e+01
Iter: 4 Ap: 6.29e-01 Pobj: 0.0000000e+00 Ad: 8.78e-01 Dobj: 4.4841314e+01
Iter: 5 Ap: 7.05e-01 Pobj: 0.0000000e+00 Ad: 8.19e-01 Dobj: 4.0867222e+01
Iter: 6 Ap: 4.49e-01 Pobj: 0.0000000e+00 Ad: 8.96e-01 Dobj: 3.3956443e+01
Iter: 7 Ap: 7.06e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 2.4975717e+01
Iter: 8 Ap: 7.70e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 1.1624169e+01
Iter: 9 Ap: 8.19e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.4118886e+00
Iter: 10 Ap: 6.94e-01 Pobj: 0.0000000e+00 Ad: 8.81e-01 Dobj: 1.5336118e+00
Iter: 11 Ap: 3.86e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.3983330e-01
Iter: 12 Ap: 7.22e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.1824039e-01
Iter: 13 Ap: 4.34e-01 Pobj: 0.0000000e+00 Ad: 5.75e-01 Dobj: 2.6649112e-01
Iter: 14 Ap: 1.93e-01 Pobj: 0.0000000e+00 Ad: 4.14e-01 Dobj: 8.8738723e-02
Iter: 15 Ap: 9.37e-02 Pobj: 0.0000000e+00 Ad: 1.25e-01 Dobj: -1.0676615e-02
Declaring primal infeasibility.
Success: SDP is primal infeasible
Certificate of primal infeasibility: a’*y=-1.00000e+00, ||A’(y)-Z||=1.64882e-10

Result is:

  • Solver : CSDP

  • Status
    Result count : 1
    Termination status : INFEASIBLE
    Message from the solver:
    “Problem is primal infeasible.”

  • Candidate solution (result #1)
    Primal status : INFEASIBLE_POINT
    Dual status : INFEASIBILITY_CERTIFICATE
    Objective value : 0.00000e+00
    Dual objective value : 1.00000e+00

  • Work counters
    Solve time (sec) : 2.38339e+02


You get the SOS certificate by construction here, it is the sum of the squares of the entries of the vector P * y multiplied by (xL^2 - x[1]^2) * (xL^2 - x[2]^2) * (xL^2 - x[3]^2)
The certificate used by default which is Putinar’s certificate only considers SOS polynomials multiplied by xL^2 - x[i]^2 alone.
Here, you need to multiply by the product of the three inequalities so you’d need Schumdgen’s certificate, see Certificate · SumOfSquares

To answer to the comment There should be a more elegant method to do this!, you can use SemidalgebraicSets.basis_semialgebraic_sets(SemidalgebraicSets.FullSpace, [xL^2 - x[i]^2 for i in 1:3]) I believe

I am sorry, please be patient with me, but I did not understand your answer. As I understand Putinar’s Positvstellensatz, under certain compactness conditions, a positive function, f can be expressed as
f= sig_0(x) + \sum_i sig_i(x)*g_i(x)
with the sig_i(x) being SOS polynomials, and the g_i(x) are polynomials that define the basic semi-algebraic set. Is it the way SumOfSquares translate the problem, that requires the use of Schmudgen’s Positivstellensatz? As I understand Schmudgen’s certificate requires ALL different products of the g_i(x). Is this automatically generated in SumOfSquares?