# Factoring a quadratic form associated with a $4 \times 4$ matrix into a sum of four squares

The matrix W = $$\left( \begin{array}{cccc} 5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 & 10 \\ \end{array} \right)$$

has a quadratic form associated with it given by:

$$(x3 + x4)^2 + (x3 + 2 x4)^2 + (x1 + x2 + 2 x3 + x4)^2 + (2 x1 + 3 x2 + 2 x3 + 2 x4)^2$$

Now if I execute the following:

w = ( {
{5, 7, 6, 5},
{7, 10, 8, 7},
{6, 8, 10, 9},
{5, 7, 9, 10}
} );

x = {x1, x2, x3, x4};

Transpose[x].w.x // Expand

(* 5 x1^2 + 14 x1 x2 + 10 x2^2 + 12 x1 x3 + 16 x2 x3 + 10 x3^2 + 10 x1 x4 + 14 x2 x4 + 18 x3 x4 + 10 x4^2 *)

Transpose[x].w.x == (x3 + x4)^2 + (x3 + 2 x4)^2 + (x1 + x2 + 2 x3 + x4)^2 + (2 x1 + 3 x2 + 2 x3 + 2 x4)^2 // Simplify

(* True *)


I am wondering is there a way of factorising the initial answer to get the sum of four squares form directly without the need to take the answer and do a comparison with the sum of squares?

I've played around with Expand, Factor, Together, Apart, Simplify etc. but don't seem to get anywhere.

• This will get a sum of squares, but not the one that you want: In[27]:= ch = CholeskyDecomposition[w]; sos = Total[(Array[x, 4].Transpose[ch])^2] Out[28]= (x[2]/Sqrt[5] - (2 x[3])/Sqrt[5])^2 + x[4]^2/2 + (Sqrt[2] x[3] + (3 x[4])/Sqrt[2])^2 + (Sqrt[5] x[1] + ( 7 x[2])/Sqrt[5] + (6 x[3])/Sqrt[5] + Sqrt[5] x[4])^2 Feb 4, 2022 at 15:13
• Feb 5, 2022 at 5:49
• @MichaelE2 - I'm very familiar with the Higham et al paper. I have been using some of the techniques in there for other things I am doing. The matrix I used in this question is the quite well known Wilson Matrix which crops up a lot in this sort of research. Feb 5, 2022 at 10:08
• The context may help others and thereby help you, especially since they seem to explain how they solved this very problem (using Mathematica, apparently). I'm not familiar with the paper, btw, and just happened upon it & glanced through it. Feb 5, 2022 at 16:15

Here's a Euclid-algorithm approach, which works on the OP's example. The basic idea is to reduce the matrix with integer-scalar symmetric row and column operations. Probably something about this is already known, but I don't know it. (Maybe ask on MathOverflow?) I don't have time to make it into a proper function right now. The steps below reflect an assumption that the input matrix w is symmetric. Other than that, I don't know how to tell when a problem is unsolvable, if unsolvable problems exist. (I did this one by hand, one row-column op at a time, trying to reduce to diagonal entries somewhat like row-reduction. The trick here is that we cannot simply swap rows, but have to swap the corresponding columns, too. The effect is to swap the diagonal entries. So you want to swap the absolutely least diagonal entry to the top-left of the current row-column one is working on. So I have some notion of an algorithm, but I don't have time right now to debug it and test its robustness. And if it is an algorithm, as I said, it's probably already known.)

(* elementary row ops *)
rop // ClearAll;
rop[{r1_Integer, r2_Integer}, m_, n_Integer] := (* add m x r1 to r2 *)
SparseArray[{{i_, i_} :> 1, {r2, r1} :> m}, {n, n}];
rop[{r1_Integer, r2_Integer}, n_Integer] :=        (* swap r1 <> r2 *)
SparseArray[{{i_, i_} /; i != r1 && i != r2 :> 1, {r1, r2} :>
1, {r2, r1} :> 1}, {n, n}];
rop[r1_Integer, m_, n_Integer] :=                  (* scale r1 by m *)
SparseArray[{{i_, i_} /; i != r1 :> 1, {r1, r1} :> m}, {n, n}];

{w // MatrixForm,             (* display new mat *)
a = rop[{1, 2}, -1, 4];     a . w . a\[Transpose], (* row/col = 1 *)
a = rop[{1, 2}, 4] . a;     a . w . a\[Transpose], (* swap 1 <> 2 *)
a = rop[{1, 2}, -2, 4] . a; a . w . a\[Transpose],
a = rop[{1, 3}, -2, 4] . a; a . w . a\[Transpose],
a = rop[{1, 4}, -2, 4] . a; a . w . a\[Transpose],
a = rop[{2, 3}, -2, 4] . a; a . w . a\[Transpose], (* row/col = 2 *)
a = rop[{2, 4}, -1, 4] . a; a . w . a\[Transpose],
a = rop[{3, 4}, -1, 4] . a; a . w . a\[Transpose], (* row/col = 3 *)
a = rop[{3, 4}, 4] . a;     a . w . a\[Transpose], (* swap 3 <> 4 *)
a = rop[{3, 4}, -1, 4] . a; a . w . a\[Transpose]  (* row/col = 4 *)
} // Map@MatrixForm

a // MatrixForm

p1 = w . vars . vars
p2 = (a . w . a\[Transpose]) . (Inverse[a\[Transpose]] .
vars) . (Inverse[a\[Transpose]] . vars)
(*
x[1] (5 x[1] + 7 x[2] + 6 x[3] + 5 x[4]) +
x[2] (7 x[1] + 10 x[2] + 8 x[3] + 7 x[4]) +
x[3] (6 x[1] + 8 x[2] + 10 x[3] + 9 x[4]) +
x[4] (5 x[1] + 7 x[2] + 9 x[3] + 10 x[4])

(x[3] + x[4])^2 +                     (* same as OP's answer :) *)
(x[1] + x[2] + 2 x[3] + x[4])^2 +
(x[3] + 2 x[4])^2 +
(2 x[1] + 3 x[2] + 2 x[3] + 2 x[4])^2
*)

p2 - p1 // Expand

(*  0  *)


First answer: Before it was clear we were to work over the integers

I was wondering how to put....into the nice sum of squares form

How about, how to put....into the nice a sum of squares form?

w = ({{5, 7, 6, 5}, {7, 10, 8, 7}, {6, 8, 10, 9}, {5, 7, 9, 10}});
vars = {x1, x2, x3, x4};
p1 = w . vars . vars;

{eval, evec} = Eigensystem[w];
evec = Normalize /@ evec;
p2 = (Sqrt[eval] evec . vars // ExpandAll)^2 // RootReduce // Total;
Collect[p1 - p2, vars] // RootReduce
(*  0  <-- p2 is the answer *)

• I started thinking I sort of asked the impossible. You need to know the answer ahead of time to get the succinct sum of four squares where the coefficients in those squares are integers. Feb 4, 2022 at 19:38
• @1729taxi Maybe, but maybe there's a way...Someone else might come up with it. Feb 4, 2022 at 19:44
• @1729taxi I added a way that works for your example. It should for other matrices, but I'm not sure it works in general. Might need to prove diagonal keeps getting smaller after every $k$ steps for some $k$ (like the dimension of the matrix). Feb 5, 2022 at 0:14
• That is very interesting. I need to think about this - I thank you a lot for this. I'll accept this as an answer . Feb 5, 2022 at 10:27
• Could you elaborate further on the sentence "So you want to swap the absolutely least diagonal entry to the top-left of the current row-column one is working on." ----- Also, what are the comments (* row/col = 1 ) and ( row/col = 3 *) trying to convey? Feb 5, 2022 at 13:34

As you may know, in a Eigenvector base the quadratic form is purely quadratic. Therefore we may get the eigenvectors, transform x1,..x4 to the eigenvectors and the quadratic form, written in the new vectors, will be purely quadratic (although different from your example):

The (normalized) eigenvectors and values:

w = ({{5, 7, 6, 5}, {7, 10, 8, 7}, {6, 8, 10, 9}, {5, 7, 9, 10}});
x = {x1, x2, x3, x4};
{eval, evec} = Eigensystem[w] // N;
evec = Normalize /@ evec;


x1..x4 written in the eigenvector base:

nvec = evec . x;


We can now check if the new vectors really diagonalize the quadratic form (I use Chop to get ride of small numeric errors):

nvec . DiagonalMatrix[eval] . nvec == x . w . x // Simplify // Chop
(* True *)


If we call the components of n1,..n4 we see that the above is purely quadratic:

Clear[n1, n2, n3, n4]
{n1, n2, n3, n4} . DiagonalMatrix[eval] . {n1, n2, n3, n4} //
Expand // Chop


Or written in the original x1,..x4:

{n1, n2, n3, n4} . DiagonalMatrix[eval] . {n1, n2, n3, n4} /.
Thread[{n1, n2, n3, n4} -> nvec] // Simplify // Chop


• I understand that but I was wondering how to put the 5 x1^2 + 14 x1 x2 + 10 x2^2 + 12 x1 x3 + 16 x2 x3 + 10 x3^2 + 10 x1 x4 + 14 x2 x4 + 18 x3 x4 + 10 x4^2 into the nice sum of squares form. Feb 4, 2022 at 16:33
• The Choleskey approach can give a sum of the form a1*p1^2+...+a4*p4^2 where the aj are rationals and the pj are linear with rational coefficients. Now whether it is possible to do better in general is something I don't know offhand. Feb 4, 2022 at 16:47
• You could also use the singular value decomposition SVD. Feb 4, 2022 at 16:56

The answer posted by @MichaelE2 reminded me that I'd encountered this problem in past. It was in the context of what is called a "cogredient transformation" of a matrix.

One approach is to use an LU decomposition. But this needs to maintain ordering, that is, avoid pivots. So the method I show below will work provided no zero pivots are encountered, which in turn means it will be fine whenever the symmetric matrix has full rank.

LUDecomposition of an integer matrix will pivot based on the smallest candidate. To enforce no pivoting we multiply rows and columns by progessively increasing powers of a sufficiently large integer. We later use the inverse of this transformation.

Further details are noted in at the link provided above. Here I show the code relevant for this example.

mat = ({{5, 7, 6, 5}, {7, 10, 8, 7}, {6, 8, 10, 9}, {5, 7, 9, 10}});
max = 2*Max[Abs[mat]];
diag = max^Range[0, 3];
mult = DiagonalMatrix[diag];
invmult = DiagonalMatrix[1/diag];
newmat = mult . mat . mult;
lud = LUDecomposition[newmat];
l = LowerTriangularize[lud[[1]], -1] + IdentityMatrix[4];
diag = Diagonal[lud[[1]]];
x = {x1, x2, x3, x4};
sos = diag . Together[(x . invmult . l)^2]

(* Out[138]=
1/5 (x2 - 2 x3)^2 + x4^2/2 + 1/2 (2 x3 + 3 x4)^2 +
1/5 (5 x1 + 7 x2 + 6 x3 + 5 x4)^2 *)


This is not quite the desired result, but certainly is a sum of squares with rationals only appearing as "weights", that is, multipliers of the squares.