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")
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
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
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?