# Expressing a polynomial as a sum of squares

I encountered a degree 4 polynomial in 8 variables f(a1,a2,a3,a4,b1,b2,b3,b4) that I suspect can be written as a sum of squares. While sostools in Matlab would find such a sum of squares decomposition, I am wondering whether a similar package exists for Mathematica.

f= 2 a1^2 b1^2 + a2^2 b1^2 + a3^2 b1^2 + 2 a1 a2 b1 b2 +
2 a3 a4 b1 b2 + a1^2 b2^2 + 2 a2^2 b2^2 + a4^2 b2^2 +
2 a1 a3 b1 b3 + 2 a2 a4 b1 b3 - 4 a2 a3 b2 b3 + 4 a1 a4 b2 b3 +
a1^2 b3^2 + 2 a3^2 b3^2 + a4^2 b3^2 + 4 a2 a3 b1 b4 -
4 a1 a4 b1 b4 + 2 a1 a3 b2 b4 + 2 a2 a4 b2 b4 + 2 a1 a2 b3 b4 +
2 a3 a4 b3 b4 + a2^2 b4^2 + a3^2 b4^2 + 2 a4^2 b4^2;

• What exactly do you mean with a sum of squares? Are you trying to simplify your expression or do you mean squares of polynomials? If so, you might have to use Grobner Basis. – Mirko Aveta Dec 9 '16 at 14:55
• I mean sum of squares of polynomials, so indeed one would use Grobner Basis to do this. SOSTOOLS does this in Matlab, i want to know if there is an analogous code in Mathematica. – pizzazz Dec 9 '16 at 15:04
• Take a look at this. I had a similar problem last year and back then Daniel Lichtblau illuminated me with the use of Grobner bases. mathematica.stackexchange.com/questions/103007/… – Mirko Aveta Dec 9 '16 at 15:07
• Nothing built-in, IIRC; someone will have to implement this or a more modern variant of it. – J. M. will be back soon Dec 9 '16 at 16:48

Here is a somewhat heuristic approach, that relies in this case on the SOS having integer coefficients. I will indicate an alteration that gives a numerical approximation in the general case.

ff = 2 a1^2 b1^2 + a2^2 b1^2 + a3^2 b1^2 + 2 a1 a2 b1 b2 +
2 a3 a4 b1 b2 + a1^2 b2^2 + 2 a2^2 b2^2 + a4^2 b2^2 +
2 a1 a3 b1 b3 + 2 a2 a4 b1 b3 - 4 a2 a3 b2 b3 + 4 a1 a4 b2 b3 +
a1^2 b3^2 + 2 a3^2 b3^2 + a4^2 b3^2 + 4 a2 a3 b1 b4 -
4 a1 a4 b1 b4 + 2 a1 a3 b2 b4 + 2 a2 a4 b2 b4 + 2 a1 a2 b3 b4 +
2 a3 a4 b3 b4 + a2^2 b4^2 + a3^2 b4^2 + 2 a4^2 b4^2;


Extract variables. It is homogeneous of degree 4 so we know we seek squares involving products of two (not necessarily distinct) variables.

allvars = Variables[ff]


Create a symbolic matrix for the bilinear form. Use it to make a symbolic such form. Extract the terms, create a new variable for each one.

mat = Array[x, {n, n}] /. x[i_, j_] /; i < j :> x[j, i];
mvars = Flatten[mat];
terms = Apply[List, ff] /. n_Integer*ab_ :> ab
tvars = Array[t, Length[terms]];


Now form the difference between this symbolic polynomial and the one of interest. Replace the terms by their new variables. We'll set all of those to zero to solve for the matrix variables.

reps = Thread[terms -> tvars];
diffpoly = qpoly - ff /. reps;


Now we separate out terms to get a linear system. The "constant" part is independent of terms of interest and we'll just set variables therein to zero. In general this might be a problem but it works fine for this example. We have that constant part, and the matrix variables within it, only because we included more quadratic monomials than we really needed.

linpolys = Normal[CoefficientArrays[diffpoly, tvars]];
unneeded = Cases[Variables[linpolys[[1]]], _x];


Now use the zeroed variables to form a smaller set of equations.

diffpoly2 = qpoly2 - ff /. reps;
linpolys2 = Normal[CoefficientArrays[diffpoly2, tvars]];

(* {0, {-2 + x[11, 11] + 2 x[15, 1], -1 + x[12, 12] + 2 x[15, 3], -1 +
x[13, 13] + 2 x[15, 6], -2 + 2 x[16, 12] + 2 x[17, 11] +
2 x[20, 2], -2 + 2 x[18, 14] + 2 x[19, 13] + 2 x[20, 9], -1 +
x[16, 16] + 2 x[21, 1], -2 + x[17, 17] + 2 x[21, 3], -1 +
x[19, 19] + 2 x[21, 10], -2 + 2 x[22, 13] + 2 x[24, 11] +
2 x[26, 4], -2 + 2 x[23, 14] + 2 x[25, 12] + 2 x[26, 8],
4 + 2 x[23, 18] + 2 x[24, 17] + 2 x[27, 5], -4 + 2 x[22, 19] +
2 x[25, 16] + 2 x[27, 7], -1 + x[22, 22] + 2 x[28, 1], -2 +
x[24, 24] + 2 x[28, 6], -1 + x[25, 25] + 2 x[28, 10], -4 +
2 x[30, 13] + 2 x[31, 12] + 2 x[33, 5],
4 + 2 x[29, 14] + 2 x[32, 11] + 2 x[33, 7], -2 + 2 x[29, 18] +
2 x[31, 16] + 2 x[34, 4], -2 + 2 x[30, 19] + 2 x[32, 17] +
2 x[34, 8], -2 + 2 x[29, 23] + 2 x[30, 22] + 2 x[35, 2], -2 +
2 x[31, 25] + 2 x[32, 24] + 2 x[35, 9], -1 + x[30, 30] +
2 x[36, 3], -1 + x[31, 31] + 2 x[36, 6], -2 + x[32, 32] +
2 x[36, 10]}} *)


Solve the system that makes bmat into a correct bilinear form for this polynomial. Plug in the solutions. We do not get values for all variables so some work remains.

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


We'll adjust the "free parameters" so as to make the matrix positive definite. To that end, create a numeric function that finds the smallest eigenvalue, given numeric values for these variables.

obj[vals : {_?NumberQ ..}] :=


Now we want to maximize that smallest eigenvalue. This is fairly slow with NMaximize, faster with FindMaximum. The latter gives a weaker result and in fact does not quite deliver a viable matrix.

AbsoluteTiming[{max, vals} = NMaximize[obj[bvars], bvars]]

(* Out[255]= {111.376033, {-0.00907241340679, {x[11, 11] ->
2.01682304768, x[12, 12] -> 0.998652993047,
x[13, 13] -> 1.00192521005, x[16, 12] -> 0.995046130758,
x[16, 16] -> 0.99609862964, x[17, 11] -> 0.00336991638564,
x[17, 17] -> 2.00227052115, x[18, 14] -> 0.00906737456166,
x[19, 13] -> 0.986910720933, x[19, 19] -> 1.01376718577,
x[22, 13] -> 0.890037586566, x[22, 19] -> 0.985669306087,
x[22, 22] -> 1.00413101831, x[23, 14] -> -0.00906106269306,
x[23, 18] -> -0.00904096742103, x[24, 11] -> 0.110559251075,
x[24, 17] -> -1.98448194582, x[24, 24] -> 2.00494690946,
x[25, 12] -> 1.00921451048, x[25, 16] -> 1.01108280433,
x[25, 25] -> 1.01123628949, x[29, 14] -> -0.00906007908607,
x[29, 18] -> -0.00903941520969, x[29, 23] -> 0.00907237266086,
x[30, 13] -> 0.987345323683, x[30, 19] -> 1.02376193674,
x[30, 22] -> 0.987031961614, x[30, 30] -> 1.01561654825,
x[31, 12] -> 1.00883984845, x[31, 16] -> 1.00918739894,
x[31, 25] -> 1.01897692633, x[31, 31] -> 1.00867831993,
x[32, 11] -> -1.98669409665, x[32, 17] -> -0.032831980588,
x[32, 24] -> -0.0259559271915, x[32, 32] -> 1.99690556219}}} *)


One will notice that the smallest eigenvalue is somewhat below zero, which is bad since it means the matrix will not be quite positive semidefinite (so no positive sum-of-squares). But we also see that every value is "near" an integer. We'll round them off and use that instead.

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


Check that this gives our bilinear form.

matsolved = (bmat /. vals2);

(* 0 *)


Check that the matrix is positive semidefinite (all eigenvalues must be >= 0).

Eigenvalues[matsolved]

(* Out[446]= {4, 4, 4, 4, 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, 0, 0, 0, 0, 0, 0, 0} *)


Form the eigensystem. Normalize the eigenvectors to that the transpose and inverse matrices are the same.

{evals, evecs} = Eigensystem[matsolved];
evecs2 = Map[#/Norm[#] &, evecs];


Now use square roots of the eigenvalues to form the Cholesky square root matrix.

cholmat = Chop[Sqrt[DiagonalMatrix[evals]].evecs2];


Check the result.

Expand[quads.Transpose[cholmat].cholmat.quads - ff]

(* Out[499]= 0 *)


Here is the explicit SoS.

(quads.Transpose[cholmat]).(cholmat.quads)

(* Out[501]= (-Sqrt[2] a2 b2 + Sqrt[2] a3 b3)^2 + (a3 b1 + a4 b2 +
a1 b3 + a2 b4)^2 + (a2 b1 + a1 b2 + a4 b3 +
a3 b4)^2 + (-Sqrt[2] a1 b1 + Sqrt[2] a4 b4)^2 *)


Had this rounding trick not worked we could still have found an SoS that approximates the input polynomial. Simply work with the approximate values for that matrix, alter negative eigenvalues to be zero, and use that to obtain the Cholesky matrix. For this purpose the lines below could be used.

matsolved = (bmat /. vals);
{evals, evecs} = Eigensystem[matsolved];
evals2 = evals /. aa_ /; aa < 0 :> 0;


When the dust settles one gets an SoS correct to two or so decimal places.

An alternative I have not tried is to assume that NMinimize is failing to deliver a high precision result, and take what it gives as an initialization for FindMinimum. It is possible that a local optimization step with a good starting value might improve on matters. In principle NMinimize ought to be doing this already, but it might be using heuristics for the method or other settings that are not good for this particular problem.

• It's annoying that CholeskyDecomposition[] chokes in the positive semidefinite case; otherwise, we easily get a test for semidefiniteness for free, and is definitely cheaper than an eigendecomposition. – J. M. will be back soon Dec 11 '16 at 16:55
• @J.M. There is a reason Cholesky decompositions are not used for solving the "normal equations" for a least-squares solution, in the rank-deficient case. That's exactly when the matrix in question is only psd and not pd. – Daniel Lichtblau Dec 11 '16 at 17:02
• One could do symmetric pivoting in the semidefinite case (see e.g. Golub/van Loan) if needed, tho. Of course, that adds a bit of complication... – J. M. will be back soon Dec 11 '16 at 17:06
• That's neat, thanks! – pizzazz Dec 11 '16 at 18:35