Rewriting a polynomial product into sums and differences of squares (a type of Lagrange reduction)

Question

Say I have a product of two polynomials and I need to rewrite it into as many squares as possible using similar sums and differences. Is there a way to achieve this in Mathematica ?

E.g.

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

Note that each subtracted square is of a similar form to the square it is subtracted from, while there may be extra added squares.

I am trying to do such a rewrite for a more complicated expression, such as:

(3a - 4b + c)(9a - 24b + 22c - 8d + e)

To me it seems similar to Writing an expression as sum of squares of expressions (I am no longer sure it is similar)

Any hints please? I need this for a smaller part of a Numerical Analysis / PDEs problem.

UPDATE

Thanks to reading up on @Daniel Lichtblau's suggestion about diagonalization, I just realized this is a type of Lagrange reduction, but with some constraints (on the coefficients used in the linear combinations), such that the terms will sum up from 1 to N in the way described below. I hope this helps? https://www.encyclopediaofmath.org/index.php/Lagrange_method

UPDATE (BRUTE FORCE)

I think maybe this might be brute-forced as follows: We need a linear combination of squares of a linear combination, where the inner linear combination might have zero coefficients.

e.g. a(a-b) = 0.5*[(1a+0b)^2 - (0a+1b)^2 + (1a-1b)^2]

and

a(3a-4b+c) = 0.5[(1a+0b+0c)^2 - (0a+1b+0c)^2 + (2a-1b+0c)^2 - (0a+2b-1c)^2 + (1a-2b+1c)^2] = 0.5(a^2 - b^2 + (2a-b)^2 - (2b-c)^2 + (a-2b+c)^2)

Here is what I am thinking:

• Given product P = (a + b + ...)(a + b + c + d +...)
• List all the terms $a,b,c,d,$ etc that are available in the product P
• Create squares of random linear combinations of the terms e.g. $(2a + 0b + 1c - d)^2$ where the coefficients can be any integer.
• Create random linear combinations of the above using coefficients $1$ and $-1$ (the kicker here is that I don't know how many terms from step 3 above, to linearly combine)
• Compare result to the product P. If it is equivalent upon expansion of P, end. If not, continue.

Old Update:

Here is a simplified explanation of what I am trying to do:

The terms a, b, c, etc are some variables $x,x_{n-1},x_{n-2}$ etc. This rearrangement is to help me obtain a telescoping sum when the product is summed from n=1 to a final value N. The final result must be positive sum, an initial expression at zero and a final expression at N.

e.g. $\displaystyle\sum_{n=1}^N(x_n - x_{n-1})x_n = 0.5\displaystyle\sum_{n=1}^N(x_n^2 -x_{n-1}^2 + (x_n-x_{n-1})^2)$ $=0.5x_N^2 - 0.5x_0^2 + 0.5\displaystyle\sum_{n=1}^N(x_n-x_{n-1})^2$

And we see that, thanks to the fact that the subtractions are deducted from an expression of similar form, the final answer will have an initial form, a final form and a sum that is definitely positive.

Similarly, using the rewriting above, I would be able to have

$\displaystyle\sum_{n=1}^N(3x_n - 4x_{n-1}+x_{n-2})x_n$

$=0.5x_N^2 - 0.5x_0^2 + 0.5(2x_N-x_{N-1})^2 - 0.5(2x_1-x_0)^2 + 0.5\displaystyle\sum_{n=1}^N(x_n-2x_{n-1}+x_{n-2})^2$

I think @march's graceful attempt in the answer below might be modified to include this constraint (?) on the subtracted forms, but I am not sure how to do this.

Answer 1

Update

I will proceed assuming that the polynomials are of the specific form that it seems you are assuming---i.e., a sum of second-degree monomials. In that case, a basis for this subspace consists of all polynomials of the form (x + y)^2, where x and y range over all of the variables in the polynomial and can be the same. Therefore, we create the following function:

polynomialToSquares[poly_, vars_] := Module[
  {basis = Plus[##]^2 & @@@ Tuples[vars, {2}] // Union, params}
  , params = Block[{x}, Array[Unique[x] &, Length@basis]]
  ; basis.params /. First@SolveAlways[poly == basis.params, vars]
 ]
polynomialToSquares[poly_] := polynomialToSquares[poly, Variables@poly]

Usage:

poly = a (a + b);
polynomialToSquares[poly]
poly - % // Expand
(* (3 a^2)/2 + b^2/2 - 1/2 (a + b)^2 *)
(* 0 *)

poly = a (3 a - 4 b + 2 c);
polynomialToSquares[poly]
poly - % // Expand
(* 4 a^2 + 2 b^2 - 2 (a + b)^2 - c^2 + (a + c)^2 *)
(* 0 *)

SeedRandom[10];
poly = Union[Times @@@ Tuples[vars, {2}]].RandomInteger[{-20, 20}, 10]
polynomialToSquares[poly]
poly - % // Expand
(* 3 a^2 + 15 a b + 18 b^2 + 2 a c + 12 b c + c^2 + 10 a d - 9 b d + 20 c d - 13 d^2 *)
(* -((21 a^2)/2) + 9 b^2 + 15/2 (a + b)^2 - 16 c^2 + (a + c)^2 + 6 (b + c)^2 - (47 d^2)/2 + 5 (a + d)^2 - 9/2 (b + d)^2 + 10 (c + d)^2 *)
(* 0 *)

---

Update 2

We can also use a direct linear algebra approach. This finds a matrix basis transformation between the basis of second-degree monomials and the basis of "squared monomials", inverts the transformation, and applies it to the original polynomial written in the monomial basis. As written, it seems a little inelegant and can probably be written nicer, but here's what I've got.

ClearAll@polynomialToSquares
polynomialToSquares[poly_, vars_] := Module[{
   monomialBasis = Times @@@ Tuples[vars, {2}] // Union
   , squaredBasis = Plus[##]^2 & @@@ Tuples[vars, {2}] // Union
   , tempVars = Block[{x}, Array[Unique[x] &, 2 Length@vars]]
   , basisTransformation
  }
  , basisTransformation = Normal@Last@CoefficientArrays[
      Expand@squaredBasis /. Thread[monomialBasis -> tempVars],  tempVars
     ]
  ; Last@CoefficientArrays[poly /. Thread[monomialBasis -> tempVars], tempVars].Inverse[basisTransformation].squaredBasis
 ]
polynomialToSquares[poly_] := polynomialToSquares[poly, Variables@poly]

Usage is the same.

---

Original Post

I don't have the math at my fingertips, so I can't guarantee that this will work for all polynomials; perhaps we have to expand the "basis". In any case, consider the following.

We consider a (possibly incomplete) "basis" for polynomials in two variables:

basis = {a^2, b^2, (a + b)^2};

Then, notice that

First@SolveAlways[a (a - b) == Array[$x, 3].basis, {a, b}]
(* {$x[1] -> 3/2, $x[3] -> -(1/2), $x[2] -> 1/2} *)

yields the correct expansion:

Array[$x, 3].basis - a (a - b) /. {$x[1] -> 3/2, $x[3] -> -(1/2), $x[2] -> 1/2} // Expand
(* 0 *)

This can be generalized. Taking

basis = {a^2, b^2, c^2, (a + b)^2, (a + c)^2, (b + c)^2};

we can get a random polynomial of the specific form that it seems you are assuming (i.e. a sum of second-degree monomials) via

RandomSeed[10];
poly = RandomInteger[{-10, 10}, 9].(Times @@@ Tuples[{a, b, c}, {2}])
(* -2 a^2 - 2 a b - 7 b^2 - 2 a c - 16 b c - 3 c^2 *)

Then, using

soln = First@SolveAlways[poly == Array[$x, 6].basis, {a, b, c}];
Array[$x, 6].basis - poly /. soln // Expand
(* 0 *)

I will generalize this later when I have the time. We'll be able to package this up into a nice little routine that programmatically generates the list of variables, the basis, etc.