# Express polynomials as sum of squares

Well, while SOStools of Matlab can express a polynomials as SOS (Sum Of Squares) by the following theorem:

A multivariate polynomial p in n variables and of degree 2d is a sum of squares if and only if there exists a positive semidefinite matrix Q (often called the Gram matrix) such that $$p\left( x \right) = {z^T}Qz,$$ where $$z$$ is the vector of monomials of degree up to $$d:$$ $$z = \left[ {1,{x_1},{x_2}, \cdots ,{x_n},{x_1}{x_2}, \cdots ,x_n^d} \right].$$

I wondered if there is a similar package in Mathematica.

I found two answer in Expressing a polynomial as a sum of squares. But unfortunately, for one of these answers (the other I can't get SOS in my Mathematica), it seems not to work with all example.

Example 1. The following is Vasile Cîrtoaje's problem: $$(a^2+b^2+c^2)^2 \geqslant 3(a^3b+b^3c+c^3a)\forall a,b,c \in \mathbb{R}.$$ f = a^4 - 3 a^3 b + 2 a^2 b^2 + b^4 - 3 b^3 c + 2 a^2 c^2 + 2 b^2 c^2 - 3 a c^3 + c^4

Example 2. (one could see an example SOS here, but this is only the result)

(a^2-b*c+4*(a+b+c)^2/9)*(b^2-c*a+4*(a+b+c)^2/9)*(c^2-a*b+4*(a+b+c)^2/9) -64*(a+b+c)^6/3^6 -4*(a-b)^2*(b-c)^2*(c-a)^2 where $$a,b,c$$ are non-negative.

Thanks a lot.

There might be several viable approaches. One I showed in a prior thread "almost" works. Except I missed that there can be certain cancellations and thus oversimplified in the processing. There was another issue in handling of coefficients where I sort of assumed they would be integers in a case where rationals also work (maybe algebraics should be considered too, I'm not sure). Here I'll try to get it right, using the current example that was not handled by the prior post.

To recap a bit, we have a homogeneous polynomial of degree 4. So we'll need to work with quadratic monomials as these are all that can appear in any sum of squares rewrite.

I'll show some of the shorter intermediate results along the way since that might help for understanding the method.

poly = a^4 - 3 a^3 b + 2 a^2 b^2 + b^4 - 3 b^3 c + 2 a^2 c^2 +
2 b^2 c^2 - 3 a c^3 + c^4;
allvars = Variables[poly]

(* Out[256]= {a, b, c}

Out[257]= {a^2, a b, b^2, a c, b c, c^2}

Out[258]= 6 *)


Get all monomials in poly and create a list of surrogate replacements for each (we'll need it later).

terms = Apply[List, poly] /. n_Integer*ab_ :> ab
tvars = Array[t, Length[terms]];

(* Out[134]= {a^4, a^3 b, a^2 b^2, b^4, b^3 c, a^2 c^2, b^2 c^2,
a c^3, c^4} *)


A sum of squares representation is equivalent to a having a matrix m such that quads.m.m^t.quads = poly (here m^t is meant as the transpose of m). What this means is that m.m^t is positive semidefinite. So we form a symmetric matrix mat that we hope will have a Cholesky or similar decomposition into the desired m.m^t.

mat = Array[x, {n, n}] /. x[i_, j_] /; i < j :> x[j, i];
mvars = Flatten[mat];


Expand out the quadratic form representation.

qpoly = Expand[quads . mat . quads]

(* Out[263]=
a^4 x[1, 1] + 2 a^3 b x[2, 1] + a^2 b^2 x[2, 2] + 2 a^2 b^2 x[3, 1] +
2 a b^3 x[3, 2] + b^4 x[3, 3] + 2 a^3 c x[4, 1] +
2 a^2 b c x[4, 2] + 2 a b^2 c x[4, 3] + a^2 c^2 x[4, 4] +
2 a^2 b c x[5, 1] + 2 a b^2 c x[5, 2] + 2 b^3 c x[5, 3] +
2 a b c^2 x[5, 4] + b^2 c^2 x[5, 5] + 2 a^2 c^2 x[6, 1] +
2 a b c^2 x[6, 2] + 2 b^2 c^2 x[6, 3] + 2 a c^3 x[6, 4] +
2 b c^3 x[6, 5] + c^4 x[6, 6] *)


Next we make a polynomial of the differences between this and the given input.

diffpoly = qpoly - poly /. reps

(* Out[270]= -t[1] + 3 t[2] - 2 t[3] - t[4] + 3 t[5] - 2 t[6] - 2 t[7] +
3 t[8] - t[9] + t[1] x[1, 1] + 2 t[2] x[2, 1] + t[3] x[2, 2] +
2 t[3] x[3, 1] + 2 a b^3 x[3, 2] + t[4] x[3, 3] + 2 a^3 c x[4, 1] +
2 a^2 b c x[4, 2] + 2 a b^2 c x[4, 3] + t[6] x[4, 4] +
2 a^2 b c x[5, 1] + 2 a b^2 c x[5, 2] + 2 t[5] x[5, 3] +
2 a b c^2 x[5, 4] + t[7] x[5, 5] + 2 t[6] x[6, 1] +
2 a b c^2 x[6, 2] + 2 t[7] x[6, 3] + 2 t[8] x[6, 4] +
2 b c^3 x[6, 5] + t[9] x[6, 6] *)


Get all possible quartics in the variables. Make a replacement rule that turns each into a distinct new variable.

quartics = Union[Flatten[Outer[Times, quads, quads]]];
qrule = Thread[quartics -> Array[qq, Length[quartics]]];


Here is one place where the prior post oversimplified matters. What we want to do now is obtain the "constant" term of the linear map in tvars defined by diffpoly. We need to have it removed, that is, solve for matrix elements so as to force this to be zero. In particular, we need to force each quartic monomial in this term to vanish. This will give linear equations in the matrix elements, that we then solve.

linpolys = Normal[CoefficientArrays[diffpoly, tvars]]
morelinpolys =
Normal[CoefficientArrays[linpolys[[1]] /. qrule,
quartics /. qrule]][[2]] /. 0 -> Nothing;
qsolns = Solve[morelinpolys == 0]

(* Out[272]= {2 a b^3 x[3, 2] + 2 a^3 c x[4, 1] + 2 a^2 b c x[4, 2] +
2 a b^2 c x[4, 3] + 2 a^2 b c x[5, 1] + 2 a b^2 c x[5, 2] +
2 a b c^2 x[5, 4] + 2 a b c^2 x[6, 2] +
2 b c^3 x[6, 5], {-1 + x[1, 1],
3 + 2 x[2, 1], -2 + x[2, 2] + 2 x[3, 1], -1 + x[3, 3],
3 + 2 x[5, 3], -2 + x[4, 4] + 2 x[6, 1], -2 + x[5, 5] + 2 x[6, 3],
3 + 2 x[6, 4], -1 + x[6, 6]}}

Out[274]= {{x[3, 2] -> 0, x[4, 1] -> 0, x[5, 1] -> -x[4, 2],
x[5, 2] -> -x[4, 3], x[6, 2] -> -x[5, 4], x[6, 5] -> 0}} *)


Now reexpand the difference polynomial, using these solutions to simplify.

qpoly2 = Expand[quads . mat . quads] /. qsolns;
diffpoly2 = qpoly2 - poly /. reps

(* Out[282]= {-t[1] + 3 t[2] - 2 t[3] - t[4] + 3 t[5] - 2 t[6] -
2 t[7] + 3 t[8] - t[9] + t[1] x[1, 1] + 2 t[2] x[2, 1] +
t[3] x[2, 2] + 2 t[3] x[3, 1] + t[4] x[3, 3] + t[6] x[4, 4] +
2 t[5] x[5, 3] + t[7] x[5, 5] + 2 t[6] x[6, 1] + 2 t[7] x[6, 3] +
2 t[8] x[6, 4] + t[9] x[6, 6]} *)


Again extract coefficient arrays. Note that the first term is now zero.

linpolys2 = Normal[CoefficientArrays[diffpoly2, tvars]]

(* ut[283]= {{0}, {{-1 + x[1, 1],
3 + 2 x[2, 1], -2 + x[2, 2] + 2 x[3, 1], -1 + x[3, 3],
3 + 2 x[5, 3], -2 + x[4, 4] + 2 x[6, 1], -2 + x[5, 5] + 2 x[6, 3],
3 + 2 x[6, 4], -1 + x[6, 6]}}} *)


We reform our matrix using the solutions we have found. This of course reduces our variable pool.

solns = Solve[linpolys2[[2]] == 0];
bmat = mat /. solns[[1]] /. qsolns[[1]]
bvars = Variables[bmat]

(* Out[285]= {{1, -(3/2), 1 - 1/2 x[2, 2], 0, -x[4, 2],
1 - 1/2 x[4, 4]}, {-(3/2), x[2, 2], 0,
x[4, 2], -x[4, 3], -x[5, 4]}, {1 - 1/2 x[2, 2], 0, 1,
x[4, 3], -(3/2), 1 - 1/2 x[5, 5]}, {0, x[4, 2], x[4, 3], x[4, 4],
x[5, 4], -(3/2)}, {-x[4, 2], -x[4, 3], -(3/2), x[5, 4], x[5, 5],
0}, {1 - 1/2 x[4, 4], -x[5, 4], 1 - 1/2 x[5, 5], -(3/2), 0, 1}}

Out[286]= {x[2, 2], x[4, 2], x[4, 3], x[4, 4], x[5, 4], x[5, 5]} *)


We now want to solve for the remaining variables in such a way as to force the matrix to be positive semidefinite. One equivalent way to state this is that we want the smallest eigenvalue to be nonnegative. So we use an objective function that maximizes it.

obj[vals : {_?NumberQ ..}] :=
{max, vals} = FindMaximum[obj[bvars],

(* During evaluation of In[294]:= FindMaximum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

Out[295]= {-3.16655*10^-7, {x[2, 2] -> 2.99999, x[4, 2] -> -1.49999,
x[4, 3] -> 1.49999, x[4, 4] -> 3., x[5, 4] -> -1.49999,
x[5, 5] -> 2.99999}} *)


From this we can surmise that we have a "nice" exact result with rational values. I will remark that not all randomized runs give this good of a result.

vals2 = vals /. Rule[a_, b_] :> Rule[a, Round[b, 1/100]]

(* Out[313]= {x[2, 2] -> 3, x[4, 2] -> -(3/2), x[4, 3] -> 3/2,
x[4, 4] -> 3, x[5, 4] -> -(3/2), x[5, 5] -> 3} *)


Now we can show the numeric matrix mat. We check that eigenvalues are nonnegative and that it corredtly reconstructs poly as a bilinear form.

matsolved = (bmat /. vals2)
Eigenvalues[matsolved]

(* Out[319]= {{1, -(3/2), -(1/2), 0, 3/2, -(1/2)}, {-(3/2), 3,
0, -(3/2), -(3/2), 3/2}, {-(1/2), 0, 1, 3/
2, -(3/2), -(1/2)}, {0, -(3/2), 3/2, 3, -(3/2), -(3/2)}, {3/
2, -(3/2), -(3/2), -(3/2), 3, 0}, {-(1/2), 3/2, -(1/2), -(3/2), 0,
1}}

Out[320]= {6, 6, 0, 0, 0, 0}

Out[321]= 0 *)


Use the eigensystem to construct the equivalent of the Cholesky decomposition for mat; this gives the desired m matrix. Note that (I missed this in the first iteration) since there are repeated eigenvalues, the corresponding eigenspaces need not be orthogonalized. So we do this explicitly.

{evals, evecs} = Eigensystem[matsolved];
evecs2 = Orthogonalize[evecs];
cholmat = Sqrt[DiagonalMatrix[evals]] . evecs2

(* Out[335]= {{-(1/2), 3/2, -(1/2), -(3/2), 0, 1}, {Sqrt[3]/
2, -(Sqrt[3]/2), -(Sqrt[3]/2), -(Sqrt[3]/2), Sqrt[3], 0}, {0, 0, 0,
0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0,
0}} *)


Last we form the explicit sum of squares and check that it recovers the input polynomial.

sos = (quads . Transpose[cholmat]) . (cholmat . quads)
Expand[sos - poly]

(* Out[341]= ((Sqrt[3] a^2)/2 - 1/2 Sqrt[3] a b - (Sqrt[3] b^2)/2 -
1/2 Sqrt[3] a c + Sqrt[3] b c)^2 + (-(a^2/2) + (3 a b)/2 - b^2/
2 - (3 a c)/2 + c^2)^2

Out[342]= 0 *)


Okay, that wasn't quite last. We can prettify it too.

Map[Factor[FactorTerms[#]] &, sos]

(* Out[345]=
3/4 (a^2 - a b - b^2 - a c + 2 b c)^2 +
1/4 (a^2 - 3 a b + b^2 + 3 a c - 2 c^2)^2 *)


I hope some of this made sense.

• Thank you very much sir, +1. I will try it with some others. Jun 27 at 22:58
• Sir, maybe every polynomial has a little difference in SOS code? Because when I try poly = a^4+b^4+c^4-a^2*b^2-b^2*c^2-c^2*a^2, I get i.imgur.com/G581OzP.png Jun 27 at 23:19
• @tthnew Thanks for catching that. In addition to the typos (ff from old post instead of poly) I had a glitch with the final matrix factorization. We need to orthogonalize the eigenvectors. The edited code should work better. Jun 28 at 14:51
• Thank you for editing. It works much better now. But for some expressions like poly=a^4+b^4+c^4-a^3*b-b^3*c-c^3*a, the cholmat matrix is still... like: i.imgur.com/YvaSPpF.png . Maybe when we round some value to use instead, for this example it make large coefficients. Jun 28 at 14:58
• Actually this is a more difficult situation. There is not a unique solution to the optimization, hence the one selected is not going to rationalize nicely. I'll give it some thought but I'm not optimistic I'll get anything better than just an approximate numeric solution. Jun 28 at 15:08