# Writing an expression as sum of squares of expressions

Suppose we have a symmetric homogeneous polynomial expression $P$ in $X=(x_1,\cdots, x_n)$. I want to check whether there are functions $g(X)$ so that $P$ is of the form $\sum _{1\le i<j\le n} g_{ij}(X) (x_i-x_j)^2$. Is there a function to find $g$?

As an example:

$2 a^4 +2 b^4 + 2 c^4+ 4 a^3 b + 4 ab^3 + 4 a^3 c + 4 a c^3 + 4 b^3 c + 4 b c^3 - 4 a^2 b c- 4 ab^2 c - 4 a b c^2 - 6 a^2 b^2 - 6 b^2 c^2 - 6 a^2 c^2= (a^2+b^2+2c^2+6ab)(a-b)^2+(b^2+c^2+2a^2+6bc)(b-c)^2+(c^2+a^2+2b^2+6ca)(c-a)^2$

In Mathematica:

2 a^4 + 2 b^4 + 2 c^4 + 4 a^3 b + 4 a b^3 + 4 a^3 c + 4 a c^3 +  4 b^3 c + 4 b c^3 -
4 a^2 b c - 4 a b^2 c - 4 a b c^2 - 6 a^2 b^2 - 6 b^2 c^2 - 6 a^2 c^2 ==
(a^2 + b^2 + 2 c^2 + 6 a b) (a - b)^2 + (b^2 + c^2 + 2 a^2 + 6 b c) (b - c)^2 +
(c^2 + a^2 + 2 b^2 + 6 c a) (c - a)^2

• А necessary condition which can be checked easily: fix a variable, say $x_1$, and put in the expression $x_2=x_3=\ldots=x_n$, then the result should be divisible by $(x_1-x_2)^2$. And so on for all variables. Commented Jun 2, 2013 at 17:52

To get the g's one can proceed in steps. The first is not entirely trivial, but we have code for it.

(1) Find a Groebner basis for the difference-of-square polynomials. Also find a conversion matrix to write the basis in terms of these square differences.

(2) Reduce the input by the Groebner basis. This shows how to rewrite the input in terms of the basis polynomials.

(3) Use the conversion matrix to now rewrite in terms of the difference-of-square polynomials.

I'll give code to find the basis and conversion matrix, then illustrate usage with your example. The code is cribbed immediately from

http://forums.wolfram.com/mathgroup/archive/2011/Mar/msg00362.html

It goes back around 15 years though.

moduleGroebnerBasis[polys_, vars_, cvars_, opts___] :=
Module[{newpols, rels, len = Length[cvars], gb, j, k, rul},
rels = Flatten[
Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
newpols = Join[polys, rels];
gb = GroebnerBasis[newpols, Join[cvars, vars], opts];
rul = Map[(# :> {}) &, rels];
gb = Flatten[gb /. rul];
Collect[gb, cvars]]

conversionMatrix[polys_, vars_] :=
Module[{aa, coords, pmat, len = Length[polys], newpolys, mgb, gb,
convmat, fvar, rvars}, coords = Array[aa, len + 1];
fvar = First[coords];
rvars = Rest[coords];
pmat = Transpose[Join[{polys}, IdentityMatrix[len]]];
newpolys = pmat.coords;
mgb = moduleGroebnerBasis[newpolys, vars, coords];
gb = mgb /. Join[{fvar -> 1}, Thread[rvars -> 0]] /. 0 :> Sequence[];
convmat = Select[mgb, ! FreeQ[#, fvar] &] /. fvar -> 0;
{gb, convmat /.
Thread[rvars -> Table[UnitVector[len, j], {j, len}]]}]


Now set up the example.

ee = 2 a^4 + 2 b^4 + 2 c^4 + 4 a^3 b + 4 a b^3 + 4 a^3 c + 4 a c^3 +
4 b^3 c + 4 b c^3 - 4 a^2 b c - 4 a b^2 c - 4 a b c^2 -
6 a^2 b^2 - 6 b^2 c^2 - 6 a^2 c^2;
vars = Variables[ee];
diffsq = Flatten[
Table[(vars[[j]] - vars[[i]])^2, {j, 2, Length[vars]}, {i, 1,
j - 1}]];


Step 1:

{gb, cmat} = conversionMatrix[diffsq, vars]

(* Out[81]= {{b^2 - 2 b c + c^2, 2 a b - 2 a c - 2 b c + 2 c^2,
a^2 - 2 a c + c^2}, {{0, 0, 1}, {-1, 1, 1}, {0, 1, 0}}} *)


Step 2:

{vec, rem} = PolynomialReduce[ee, gb, vars]

(* Out[82]= {{-6 a^2 + 4 a b + 2 b^2 + 4 a c + 8 b c + 8 c^2,
2 a^2 - 6 a c - 6 c^2, 2 a^2 + 12 a c + 6 c^2}, 0} *)


Step 3:

gvec = vec.cmat

(* Out[88]= {-2 a^2 + 6 a c + 6 c^2,
4 a^2 + 6 a c, -4 a^2 + 4 a b + 2 b^2 - 2 a c + 8 b c + 2 c^2} *)


We check the result.

Expand[gvec.diffsq - ee]

(* Out[89]= 0 *)

• Table[UnitVector[len, j], {j, len}] - it's shorter to use IdentityMatrix[len], no? Commented Jun 4, 2013 at 0:26
• @J. M. Yeah. Not sure why I had it that way. Maybe just because I was thinking of a module basis as comprised of unit vectors. Commented Jun 4, 2013 at 14:05
• If not the form $x_j-x_i,$ is there any way to get SOS? Like SOStool of matlab? Commented Feb 22, 2021 at 7:29
• @tthnew could check this MSE thread Commented Feb 22, 2021 at 20:27
• I saw it, but code by Daniel Lichtblau maybe is not generate. And the code of Cesareo I can't use to get SOS. Commented Feb 23, 2021 at 1:17