Putinar certificates for SOS matrices on domains

[giofantuzzi]

Summary

Putinar cerificates to verify that a polynomial matrix $M(x)$ is nonnegative on a semialgebraic set $K$ are more expensive than necessary due to scalarization.

Example

Problem

Test if a 2x2 univariate polynomial matrix $M(x)$ is positive semidefinite when $x$ satisfies $1-x^2\geq 0$.

Matrix-valued certificates

At relaxation order $\omega$, the matrix-valued Putinar SOS strengthening of the problem is

$$M(x) = \Sigma_0(x) + (1-x^2) \Sigma_1(x)$$

where $\Sigma_0$ and $\Sigma_1$ are SOS matrices with $\deg(\Sigma_0) = 2\omega$ and $\deg(\Sigma_1)=2\omega-1$.

Scalarized version of the certificate

An equivalent formulation is to produce a nonnegativity certificate for the scalarized polynomial $m(x,y) = y^\top M(x) y$ with $y \in \mathbb{R}^2$. The desired certificate is

$$m(x,y) = \sigma_0(x,y) + (1-x^2) \sigma_1(x,y)$$

where $\sigma_0(x,y)= y^\top \Sigma_0(x) y^\top$ and $\sigma_1(x,y)= y^\top \Sigma_1(x) y^\top$ are SOS polynomials. Crucially, these polynomials are homogeneous quadratic in $y$.

The issue

SumOfSquares.jl works with the scalarized certificate but ignores the homogeneous quadratic structure of the SOS polynomials. This leads to more expensive SOS certificates than expected.

A minimum working example follows.

# half degree of polynomial matrix inequality
ω = 2 
# Random 2x2 SOS polynomial matrix of degree 2ω
@polyvar x
MON = monomials(x, 0:ω)
V = rand(2, 2*length(MON)) * kron( I(2), MON )
M = Symmetric(V'*V)
# Domain
K = @set 1 - x^2 ≥ 0
# SOS program
model = SOSModel(Mosek.Optimizer)
@constraint(model, PMI, M ∈ PSDCone(), domain=K)
optimize!(model)
certificate_monomials(PMI)

The output of certificate_monomials(PMI) is:

10-element MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}:
 ##418
 ##416
 ##416##418
 x##418
 x##416
 x²
 x##416##418
 x²##418
 x²##416
 x³

Monomials in this list are redundant and could (should?) be dropped if they are

1. independent of the $y$ variables (variables Remove block constructor #416 and Simpler mapreduce #418), or
2. quadratic in the $y$ variables.

Desired behaviour

For a polynomial matrix $M(x)$ of even degree, the constraint

@constraint(model, M ∈ PSDCone(), domain=K)

should construct a Putinar certificate for $y^\top M(x) y$ that

1. has the same degree in $x$ as $M(x)$ and
2. is homogeneous quadratic in $y$

[blegat]

Good point, I never thought of this. The polynomials σ1 indeed won't have the right form and hence when the monomials of the polynomial σ0 are computed using newton polytope, since we look at m - σ1, it's not multipartite anymore in x, y. When the user writes @constraint(model, M ∈ PSDCone(), domain=K), we see that the multiplier for the domain should be matrices so we could do better indeed.