# Write $(x^2 + y^2 + z^2)^2 - 3 ( x^3 y + y^3 z + z^3 x)$ as a sum of (three) squares of quadratic forms

The quartic form

$$(x^2 + y^2 + z^2)^2 - 3 ( x^3 y + y^3 z + z^3 x)$$

is non-negative for all real $x$, $y$, $z$, as one can check (with some effort). A theorem of Hilbert implies that there exist quadratic forms $Q_1$, $Q_2$, $Q_3$ so that

$$(x^2 + y^2 + z^2)^2 - 3( x^3 y + y^3 z + z^3 x) = Q_1^2 + Q_2^2 + Q_3^2$$

I would like to find an explicit writing of the quartic forms, with rational quadratic forms $Q_i$. Maybe more than $3$ terms are necessary.

• Why $3$ quadratic forms? What is special about $3$? – Rodrigo de Azevedo Sep 5 '17 at 15:55
• In this case ( positive quartic forms) $3$ squares will do. I don't know how HIlbert proved his result. It is not valid for more variables or higher degree I suppose, there are counter examples. So $3$ is just the number that will do. – Orest Bucicovschi Sep 5 '17 at 15:59
• I believe $2$ will do in this case. I will update my answer soon. – Rodrigo de Azevedo Sep 5 '17 at 16:14
• @Rodrigo de Azevedo: That would be pretty nice then. So you can get a positive semidefinite matrix of rank $2$ in that subspace. Interesting. If you could just produce the matrix that would do it. – Orest Bucicovschi Sep 5 '17 at 16:18

$$(x^2 + y^2 + z^2)^2 - 3 ( x^3 y + y^3 z + z^3 x)=\frac{1}{2}\sum_{cyc}(x^2-y^2-xy-xz+2yz)^2=$$ $$=\frac{1}{6}\sum_{cyc}(x^2-2y^2+z^2-3xz+3yz)^2.$$

• The first sum, is there some mistake since $-xy$ appears twice ? – Orest Bucicovschi Aug 30 '17 at 12:11
• @orangeskid It was typo. See now please. By the way, if $Q_i$ are polynomials then your statement is wrong. – Michael Rozenberg Aug 30 '17 at 12:14
• Excellent! Great answer! – Orest Bucicovschi Aug 30 '17 at 12:24
• @orangeskid Thank you! Good luck! – Michael Rozenberg Aug 30 '17 at 12:25
• @Jaideep Khare For the first solution we can write $\sum\limits_{cyc}(x^4+x^2y^2-3x^3y)=\sum\limits_{cyc}(ax^2+by^2+cz^2+dxy+exz+fyz)^2.$ Now, it's better to assume $a=-b$, $c=0$ and from here $d+e+f=0$ because the equality in $\sum\limits_{cyc}(x^4+x^2y^2-3x^3y)=0$ occurs for $x=y=z$. Now, we can write the system on $a$, $b$, $d$, $e$ and $f$ and playing with this system gives the answer. The second solution we can get by the similar way. – Michael Rozenberg Aug 30 '17 at 13:14

We are given a trivariate quartic form

$$f (x, y, z) := (x^2 + y^2 + z^2)^2 - 3 ( x^3 y + y^3 z + z^3 x)$$

which can be rewritten as the sum of $9$ monomials

$$f (x, y, z) = x^4 + y^4 + z^4 + 2 x^2 y^2 + 2 x^2 z^2 + 2 y^2 z^2 - 3 x^3 y - 3 y^3 z - 3 z^3 x$$

and we would like to write it as the sum of squares (SOS) of quadratic forms. The fewer the better.

## SOS decomposition

A general trivariate quartic form is parametrized as follows

$$\begin{array}{c} \begin{bmatrix} x^2\\ y^2\\ z^2\\ x y\\ x z\\ y z\end{bmatrix}^\top \begin{bmatrix} q_{11} & q_{12} & q_{13} & q_{14} & q_{15} & q_{16}\\ q_{12} & q_{22} & q_{23} & q_{24} & q_{25} & q_{26}\\ q_{13} & q_{23} & q_{33} & q_{34} & q_{35} & q_{36}\\ q_{14} & q_{24} & q_{34} & q_{44} & q_{45} & q_{46}\\ q_{15} & q_{25} & q_{35} & q_{45} & q_{55} & q_{56}\\ q_{16} & q_{26} & q_{36} & q_{46} & q_{56} & q_{66}\end{bmatrix} \begin{bmatrix} x^2\\ y^2\\ z^2\\ x y\\ x z\\ y z\end{bmatrix} = \\\\ = q_{11} x^4 + q_{22} y^4 + q_{33} z^4 + 2 q_{14} x^3 y + 2 q_{15} x^3 z + 2 q_{24} x y^3 + 2 q_{26} y^3 z + 2 q_{35} x z^3 + 2 q_{36} y z^3 + (2 q_{12} + q_{44}) x^2 y^2 + (2 q_{13} + q_{55}) x^2 z^2 + (2 q_{23} + q_{66}) y^2 z^2 + (2 q_{16} + 2 q_{45}) x^2 y z + (2 q_{46} + 2 q_{25}) x y^2 z + (2 q_{34} + 2 q_{56}) x y z^2\end{array}$$

Note that the $6 \times 6$ matrix above, which we denote by $\mathrm Q$, is symmetric by construction. Hence, we have $1 + 2 + \cdots + 6 = 21$ unknowns (instead of $6^2 = 36$). We recover the particular quartic $f$ if the following $15$ affine equality constraints are satisfied

$$\begin{array}{rl} q_{11} &= 1\\ q_{22} &= 1\\ q_{33} &= 1\\ 2 q_{14} &= -3\\ q_{15} &= 0\\ q_{24} &= 0\\ 2 q_{26} &= -3\\ 2 q_{35} &= -3\\ q_{36} &= 0\\ 2 q_{12} + q_{44} &= 2\\ 2 q_{13} + q_{55} &= 2\\ 2 q_{23} + q_{66} &= 2\\ 2 q_{16} + 2 q_{45} & = 0\\ 2 q_{46} + 2 q_{25} &= 0\\ 2 q_{34} + 2 q_{56} &= 0\end{array}$$

Note that we have $21 - 15 = 6$ degrees of freedom. These $15$ equality constraints can be written much more economically in the following form

$$\mathcal A (\mathrm Q) = \mathrm b$$

where $\mathcal A : \mbox{Sym}_6 (\mathbb R) \to \mathbb R^{15}$ is a linear function and $\mathrm b \in \mathbb R^{15}$.

Our goal is to find a (symmetric) positive semidefinite $6 \times 6$ matrix $\rm Q$ such that the $15$ equality constraints above are satisfied. If such a matrix $\rm Q$ exists, then there exists a ("square root") matrix $\rm R$ such that $\rm Q = R R^T$ and, thus, $f$ can be written as a SOS of quadratic forms.

## Approximate solution

We would like matrix $\rm Q$ to be low-rank, in order to have as few terms in the SOS decomposition as possible. Hence, we have the following rank-minimization problem

$$\begin{array}{ll} \text{minimize} & \mbox{rank} (\mathrm Q)\\ \text{subject to} & \mathcal A (\mathrm Q) = \mathrm b\\ & \mathrm Q \succeq \mathrm O_6\end{array}$$

Unfortunately, the objective function, $\mbox{rank} (\mathrm Q)$, is non-convex. One popular heuristic is to minimize the nuclear norm of $\rm Q$, which is a convex proxy for $\mbox{rank} (\mathrm Q)$. Fortunately, the nuclear norm of a symmetric, positive semidefinite matrix is its trace, i.e., a linear function of the matrix's (diagonal) entries. Thus, we obtain the following semidefinite program (SDP)

$$\begin{array}{ll} \text{minimize} & \mbox{tr} (\mathrm Q)\\ \text{subject to} & \mathcal A (\mathrm Q) = \mathrm b\\ & \mathrm Q \succeq \mathrm O_6\end{array}$$

which is convex and, thus, computationally tractable. The following Python + NumPy + CVXPY script solves the SDP above:

from cvxpy import *
import numpy as np

np.set_printoptions(linewidth=125)

# matrix variable
Q = Semidef(6)

# objective function
objective = Minimize( trace(Q) )

# constraints
constraints = [   Q[0,0]            ==  1,
Q[1,1]            ==  1,
Q[2,2]            ==  1,
2*Q[0,3]            == -3,
Q[0,4]            ==  0,
Q[1,3]            ==  0,
2*Q[1,5]            == -3,
2*Q[2,4]            == -3,
Q[2,5]            ==  0,
2*Q[0,1] +   Q[3,3] ==  2,
2*Q[0,2] +   Q[4,4] ==  2,
2*Q[1,2] +   Q[5,5] ==  2,
2*Q[0,5] + 2*Q[3,4] ==  0,
2*Q[3,5] + 2*Q[1,4] ==  0,
2*Q[2,3] + 2*Q[4,5] ==  0 ]

# create optimization problem
prob = Problem(objective, constraints)

# solve optimization problem
prob.solve()

print "Solution = \n", Q.value
print "Solution after rounding = \n", np.round(Q.value, 1)

# compute eigendecomposition of the solution
Lambda, V = np.linalg.eigh( Q.value )
print "Eigenvalues = \n", Lambda


This script produces the following output:

Solution =
[[  1.00052427e+00  -4.93986993e-01  -4.93986993e-01  -1.50131075e+00   8.16263094e-04   1.48128400e+00]
[ -4.93986993e-01   1.00052427e+00  -4.93986993e-01   8.16263094e-04   1.48128400e+00  -1.50131075e+00]
[ -4.93986993e-01  -4.93986993e-01   1.00052427e+00   1.48128400e+00  -1.50131075e+00   8.16263095e-04]
[ -1.50131075e+00   8.16263094e-04   1.48128400e+00   2.98994780e+00  -1.48087234e+00  -1.48087234e+00]
[  8.16263094e-04   1.48128400e+00  -1.50131075e+00  -1.48087234e+00   2.98994780e+00  -1.48087234e+00]
[  1.48128400e+00  -1.50131075e+00   8.16263095e-04  -1.48087234e+00  -1.48087234e+00   2.98994780e+00]]
Solution after rounding =
[[ 1.  -0.5 -0.5 -1.5  0.   1.5]
[-0.5  1.  -0.5  0.   1.5 -1.5]
[-0.5 -0.5  1.   1.5 -1.5  0. ]
[-1.5  0.   1.5  3.  -1.5 -1.5]
[ 0.   1.5 -1.5 -1.5  3.  -1.5]
[ 1.5 -1.5  0.  -1.5 -1.5  3. ]]
Eigenvalues =
[ -3.66863824e-04   1.62150777e-03   1.62150779e-03   4.11202771e-02   5.96370990e+00   5.96370990e+00]


The SDP solver produces a matrix $\rm Q$ whose (numerical) rank is $2$. Note that one of the eigenvalues of $\rm Q$ is slightly negative due to numerical noise.

## Exact solution

The approximate solution produced by the SDP solver suggests that the exact solution is

$$\mathrm Q = \begin{bmatrix} 1 & - \frac{1}{2} & - \frac{1}{2} & - \frac{3}{2} & 0 & \frac{3}{2}\\- \frac{1}{2} & 1 & - \frac{1}{2} & 0 & \frac{3}{2} & - \frac{3}{2}\\- \frac{1}{2} & - \frac{1}{2} & 1 & \frac{3}{2} & - \frac{3}{2} & 0\\- \frac{3}{2} & 0 & \frac{3}{2} & 3 & - \frac{3}{2} & - \frac{3}{2}\\0 & \frac{3}{2} & - \frac{3}{2} & - \frac{3}{2} & 3 & - \frac{3}{2}\\ \frac{3}{2} & - \frac{3}{2} & 0 & - \frac{3}{2} & - \frac{3}{2} & 3\end{bmatrix}$$

It is easy to verify that, indeed, this is the exact solution. Using Python + SymPy:

from sympy import *

# build Q matrix
U = Matrix([[    1,   -1,   -1,   -3,    0,    3],
[    0,    1,   -1,    0,    3,   -3],
[    0,    0,    1,    3,   -3,    0],
[    0,    0,    0,    3,   -3,   -3],
[    0,    0,    0,    0,    3,   -3],
[    0,    0,    0,    0,    0,    3]])
Q = (U + U.T) / 2

print "Determinant = ", Q.det()
print "Rank        = ", Q.rank()
print "Eigenvalues = ", Q.eigenvals()

# compute Cholesky decomposition
L = Q.cholesky()

print "6x2 square root = ", L[:,[0,1]]

#####################################
# compute the exact SOS decomposition
#####################################

x, y, z = symbols('x y z')

# quartic form
f = (x**2 + y**2 + z**2)**2 - 3 * ( x**3 * y + y**3 * z + z**3 * x)

q1 = (L[:,0].T * Matrix([x**2, y**2, z**2, x*y, x*z, y*z]))[0,0]
q2 = (L[:,1].T * Matrix([x**2, y**2, z**2, x*y, x*z, y*z]))[0,0]

print "Error polynomial = ", simplify(q1**2 + q2**2 - f)


This script produces the following output:

Determinant =  0
Rank        =  2
Eigenvalues =  {0: 4, 6: 2}
6x2 square root =  Matrix([[   1,          0],
[-1/2,  sqrt(3)/2],
[-1/2, -sqrt(3)/2],
[-3/2, -sqrt(3)/2],
[   0,    sqrt(3)],
[ 3/2, -sqrt(3)/2]])
Error polynomial =  0


Thus, a SOS decomposition of $f$ is

$$\boxed{\begin{array}{rl} f (x,y,z) &= \left( x^{2} - \frac{y^{2}}{2} - \frac{z^{2}}{2} - \frac{3}{2} x y + \frac{3}{2} y z \right)^2 + \left( \frac{\sqrt{3}}{2} \left( y^{2} - z^{2} - x y + 2 x z - y z \right) \right)^2\\ &= \frac{1}{4} \left( 2 x^{2} - y^{2} - z^{2} - 3 x y + 3 y z \right)^2 + \frac{3}{4} \left( y^{2} - z^{2} - x y + 2 x z - y z \right)^2\end{array}}$$

## References

• Oh, great, so this is a SOS method, right? – Orest Bucicovschi Sep 3 '17 at 0:37
• Yes, it is Parrilo's method. – Rodrigo de Azevedo Sep 3 '17 at 0:39
• That is indeed impressive! Did not expect to have a solution in $2$ terms. Do you think that there are just finitely many solutions of rank $2$? – Orest Bucicovschi Sep 8 '17 at 16:41
• The column space of $\rm Q$ is $2$-dimensional. Every orthonormal basis for this $2$-dimensional column space leads to a SOS decomposition with 2 terms. – Rodrigo de Azevedo Sep 8 '17 at 23:10