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

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.

• It seems from the examples that you want to do this exclusively for polynomials that are quadratic. Is that correct? Commented Sep 22, 2017 at 16:18
• @DanielLichtblau Yes, this is true. Multiplying it out gives a quadratic in at least one term. Commented Sep 22, 2017 at 16:28
• You can rewrite as a difference of sums of squares for the quadratic case by diagonalizing the quadratic form (can look up that phrase if it is not familiar). This, however, will typically give rise to square roots in the coefficients. On the plus side, it is minimal in the sense of the number of positive and negative terms it requires. Commented Sep 22, 2017 at 16:45
• To give an example of what I mean, but using approximate numbers to avoid the huge algebraic expressions, (3a - 4b + c)(9a - 24b + 22c - 8d + e) can be rewritten as a positive plus a negative square: -4.79720198444 (1. a - 0.283636811685 b - 1.32868615928 c + 0.699797681099 d - 0.0874747101373 e)^2 + 31.7972019844 (1. a - 1.74105454634 b + 0.978891920593 c - 0.271814142004 d + 0.0339767677505 e)^2 Commented Sep 22, 2017 at 16:50
• @DanielLichtblau Thanks I will check out "diagonalization". However, the result in the example above won't sum up telescopically to give a final form, an initial form and a sum that is certainly positive (please see my example in the OP). The summing is the very reason I need this rewrite. In order to be correct, every item Y^2 that is subtracted, must be subtracted from a similar item X^2, such that they will cascade in their sum from 1 to N. e.g. a result that has (2a+3b)^2 - (2b+3c)^2 + c^2 would be fine because the only subtracted term (2b+3c) is from a similar term (2a+3b) Commented Sep 22, 2017 at 18:34

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.

• Thank you! I honestly am a total noob with Mathematica, but it looks like this approach should work! I will patiently wait for your update. Thanks again. Commented Sep 20, 2017 at 1:27
• @Cogicero. I think the updated version will work nicely. Check it out. Again, note that I'm assuming that you have a polynomial consisting only of second-degree monomials. Commented Sep 20, 2017 at 3:12
• Thanks a lot, March. I am accepting this answer because you have put in a whole lot of work and I am thankful, but sadly I'm not there yet. As I said, each subtracted square is of a similar form to the square it is subtracted from. Using this code, then a (a + b) = (3 a^2)/2 + b^2/2 - 1/2 (a + b)^2 isn't correct because the 0.5(a+b)^2 is subtracted, and not from a similar form. Whereas, something like a(a-b) = 0.5(a^2 - b^2 + (a-b)^2) is correct cos the only subtraction is b^2, which is subtracted from similar form i.e. a^2. I am adding a note to the original post for more context. Commented Sep 21, 2017 at 17:37
• @Cogicero. I'm not sure how to do this manually, but you can change the basis that you want to work with. In other words, for your examples in your comment, use the basis {a^2, b^2, (a-b)^2} instead. I'm not sure how to automate it, but I'll give it a shot. In the meantime, if I didn't completely answer your question (since you did make the comment "each subtracted square is of a similar form to the square it is subtracted from" in your OP, although I skipped it because I didn't understand it), maybe hold off on accepting. Commented Sep 21, 2017 at 18:09
• @Cogicero. Your second update makes me think that this is a hard problem in general, because your squared terms can be complicated. Perhaps, if you have some idea of what terms would appear, you can use those to form the basis, and then use my code to get the expansion coefficients. Otherwise, I definitely have no idea how to proceed in automating this procedure, i.e. finding the correct basis. Commented Sep 21, 2017 at 18:12