In [1]:
using Pkg
Pkg.activate(".")
# Pkg.instantiate()

  Activating project at `~/Mathematics/Teaching/2022/Computational Group Theory/2022 WS Karlsruhe/Symbolic SOS`
In [2]:
# Pkg.build("PyCall");
In [ ]:
using SymPy
using LinearAlgebra

Example¶

When is a quadratic polynomial in one variable decomposable as sum of squares?

Solution Let us assume that $p = ax^2 + bx + c$ is an arbitrary polynomial in $\mathbb{R}[x]$. We may assume that $a > 0$.

In [ ]:
a = symbols("a", positive=true)
x, b, c = symbols("x b c", real=true);
In [ ]:
p = a*x^2 + b*x + c

Let us define the monomial basis: $m_1 = [1, x]$ (as a row vector)

In [ ]:
m₁ = transpose(Sym[1, x])

and define $P$ as

In [ ]:
P = Sym[c b/2; b/2 a]

Let us check, that the $m·P·m^T$ does what it is supposed to do:

In [ ]:
ex = m₁*P*m₁'
In [ ]:
display(simplify(ex))
simplify(ex) == p

Now we take $LDL^T$-decomposition of $P$, i.e. we find: a lower-triangular matrix $L$ and a diagonal matrix $D$ such that $$L·D·L^T = P$$

$$L·\sqrt D · (\sqrt D )^T·L^T = L·\sqrt D · \left(L·\sqrt D \right)^T = P$$

Then $\sqrt P = L\cdot\sqrt D$

In [ ]:
L, D = sympy.Matrix(P).LDLdecomposition()
L*D*L' == P

Using $L$ and $D$ we may define $$Q = L·\sqrt{D},$$ which is a matrix square root of $P$, i.e. $$Q·Q^T = \left(L\sqrt{D}\right)·\left(L\sqrt{D}\right)^T = L·\left(\sqrt{D}\sqrt{D}\right)·L^T = P.$$

Thus the columns of $Q$ provide candidates for summands in SOS-decomposition: $$ p = [1,x]·P·\begin{bmatrix}1\\x\end{bmatrix} = \left([1,x]·(L\sqrt{D})\right)·\left([1,x]·(L\sqrt{D})\right)^T = \left([1,x]Q\right)·\left([1,x]Q\right)^T $$

In [ ]:
Q = L*sqrt.(D)
Q

$p = f₁^2 + f₂^2$

In [ ]:
f = m₁*Q
In [ ]:
simplify(f[1]^2 + f[2]^2)

For this to be a sum of squares we just need that $c>0$ and $$a - \frac{b^2}{4c} ≥ 0,$$ i.e. $$b^2 - 4ac ≤ 0.$$

In [ ]:

Example¶

Formulate the sum of squares decomposition problem as the feasibility of a semi-definite program for a polynomial in one variable of degree no larger than $4$.

Solution: Let $f\in \mathbb{R}[x]$ be an arbitrary polynomial in $x$ of degree not greater than $4$, i.e. $$f = C_0 + C_1x + C_2x^2 + C_3x^3 + C_4x^4.$$

This time our monomial basis $m_2$ will consist of three elements:

In [ ]:
m₂ = [1 x x^2]
In [ ]:
P = let p = collect(symbols("p(0:3)(0:3)", real=true))
    P = Symmetric(reshape(p, 3,3))
end

We need to evaluate

In [ ]:
ex = (m₂*P*m₂')[1]

How does the resulting polynomial $m·P·m^T$ depends on $P$?

In [ ]:
collect(expand(ex), x)

Comparing this to the given (arbitrary) polynomial $f$:

In [ ]:
C = symbols("C0:5", real=true)
f = sum(C[i]*x^(i-1) for i in 1:5)
In [ ]:
pex = ex.as_poly(x)
pf = f.as_poly(x)
for (a,b) in zip(pex.coeffs(), pf.coeffs())
    display(Eq(a, b))
end

So the final semi-definite problem is as follows: $$\begin{align*} \quad & \text{feasibility}\\ \text{Subject to:}\quad & P = (p_{i,j}) \succcurlyeq 0\\ & p_{0,0} = C_0\\ & 2p_{0,1} = C_1\\ & 2p_{0,2} + p_{1,1} = C_2\\ & 2p_{1,2}= C_3\\ & p_{2,2} = C_4\\ \end{align*}$$

$\max t$ s.t. $f(x) - t \in \Sigma^2\mathbb{R}[x]$

In [ ]:
C = [1, 2, 3, 4, 5];
f = dot([1 x x^2 x^3 x^4], C)
In [ ]:
using JuMP
M = Model()
@variable M P[1:3,1:3]
@constraints(M, begin # shifting indices by one!
    P[1,1] == C[1]
    2P[1,2] == C[2]
    2P[1,3] + P[1,1] == C[3]
    2P[2,3] == C[4]
    P[3,3] == C[5]
    P in PSDCone() # P >= 0
    end
)
print(M)

Motzkin example¶

In [ ]:
x,y = symbols("x y", real=true);
Motzkin = x^4*y^2 + x^2*y^4 - 3x^2*y^2 + 1
In [ ]:
monomials(d::Int, x, y) = [x^(d-i)*y^(i) for i in 0:d]
function monomials(iter, x, y)
    m = Sym[]
    for d in iter
        append!(m, monomials(d, x, y))
    end
    return m
end

# display(monomials(3, x, y))
# monomials(2:3, x, y)
In [ ]:
m₃ = transpose(monomials(0:3, x, y))
In [ ]:
P = let l = length(m₃), p = collect(symbols("p(0:$l)(0:$l)", real=true))
    Symmetric(reshape(p, l,l))
end
In [ ]:
sympy.MatMul(Matrix(m₃), P, Matrix(m₃'))
In [ ]:
f = (m₃*P*m₃')[1]
In [ ]:
pf = sympy.Poly(f, x,y)
pm = sympy.Poly(Motzkin, x,y)
In [ ]:
for m in monomials(0:6, x, y)
#     display(m)
    display(Eq(pf.coeff_monomial(m), pm.coeff_monomial(m)))
end
In [ ]:
Motzkin

Proving that Motzkin is not a SOS¶

Solution: Suppose that Motzkin polynomial IS a sum of squares of polynomials and take one summand $f^2$ of the sum.

In [ ]:
Motzkin
In [ ]:
a₀ = symbols("a0", real=true)
a = collect(symbols("a1:10", real=true))
f = sum(c*m for (c,m) in zip([a₀;a], m₃))
In [ ]:
sympy.Poly(expand(f^2), x, y)(
    a[6] => 0, # it's a coefficient in front of x^6
    a[9] => 0, # it's a coefficient in front of y^6
    a[3] => 0, # it's a coefficient in front of x^4
    a[5] => 0, # it's a coefficient in front of y^4
    a[1] => 0, # it's a coefficient in front of x^2
    a[2] => 0, # it's a coefficient in front of y^2
)

So if Motzikn polynomial $m(x,y)$ is a sum of squares, $$ m(x,y) = \sum_i f_i^2(x,y)\enspace, $$ its coefficient in front of $x^2y^2$ is the sum of the corresponding coefficients in $f_i$: $$ -3 = \sum_i (f_i^2)^{x^2y^2} = \sum_i \left((a_4)_{(i)}\right)^2\enspace . $$

In [ ]:
In [ ]: