Writing an expression as sum of squares of expressions

Question

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

Answer 1

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 *)