# 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. – Andrew Jun 2 '13 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? – J. M. Jun 4 '13 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. – Daniel Lichtblau Jun 4 '13 at 14:05