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

$$(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.$$


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)

# quadratic forms
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

  • Pablo Parrilo, Semidefinite programming relaxations for semialgebraic problems, 2001.

  • Helfried Peyrl, Pablo Parrilo, Computing sum of squares decompositions with rational coefficients, Theoretical Computer Science, Volume 409, Issue 2, December 2008.


$$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{6}\left ( \sum\limits_{cyc}\,(\,-\,3\,yz+ 3\,xy+ y^{\,2}+ z^{\,2}- 2\,x^{\,2}\,)^{\,2} \right )$$ $$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{2}\left ( \sum\limits_{cyc}\,(\,-\,yz- zx+ 2\,xy+ z^{\,2}- x^{\,2}\,)^{\,2} \right )$$ After that, using $ab+ bc+ ca= 0\!\therefore\!a\equiv (a- b)(b- c),b\equiv (b- c)(c- a),c\equiv (c- a)(a- b)$, so $$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{14}\left ( \sum\limits_{cyc}\,(\,x^{\,2}+ 2\,y^{\,2}- 3\,z^{\,2}- 4\,xy- yz+ 5\,zx\,)^{\,2} \right )$$ $$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{14}\left ( \sum\limits_{cyc}\,(\,2\,x^{\,2}+ y^{\,2}- 3\,z^{\,2}- 5\,xy+ yz+ 4\,zx\,)^{\,2} \right )$$ $$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{14}\left ( \sum\limits_{cyc}\,(\,3\,x^{\,2}- y^{\,2}- 2\,z^{\,2}- 5\,xy+ 4\,yz+ zx\,)^{\,2} \right )$$ $$(\,\sum\limits_{cyc}\,x^{\,2}\,)^{\,2}- 3(\,\sum\limits_{cyc}\,x^{\,3}y\,)= \frac{1}{14}\left ( \sum\limits_{cyc}\,(\,3\,y^{\,2}- 2\,z^{\,2}- x^{\,2}- xy- 4\,yz+ 5\,zx\,)^{\,2} \right )$$ because we always have the following equality with $\sum\limits_{cyc}\left ( (\,a- b\,)(\,b- c\,) \right )\left ( (\,b- c\,)(\,c- a\,) \right )= 0$ $$\left ( \sum\limits_{cyc}\left ( (\,a- b\,)(\,b- c\,) \right )^{\,2} \right )(\,x^{\,2}+ y^{\,2}+ z^{\,2}\,)= \sum\limits_{cyc}\,\left ( \sum\limits_{cyc}\,\left ( x(\,a- b\,)(\,b- c\,) \right ) \right )^{\,2}$$


Furthermore, we have $\sum\limits_{cyc}\,(\,-\,3\,yz+ 3\,xy+ y^{\,2}+ z^{\,2}- 2\,x^{\,2}\,)= 0= \left ( \sum\limits_{cyc}(\,a- b\,) \right )$. To apply $$\sum\limits_{cyc}\,\frac{(\!a- b\!)^{2}}{2}\!=\!\frac{(\!a+ b- 2c\!)^{2}\!+\!3(\!a- b\!)^{2}}{4}\!=\!\frac{(\!b+ c- 2a\!)^{2}\!+\!3(\!b- c\!)^{2}}{4}\!=\!\frac{(\!c+ a- 2b\!)^{2}\!+\!3(\!c- a\!)^{2}}{4}$$ $$(\!\sum\limits_{cyc}\,x^{\,2}\!)^{\,2}- 3(\!\sum\limits_{cyc}\,x^{\,3}y\!)= \frac{(\!2\,x^{\,2}- y^{\,2}- z^{\,2}- 3\,xy+ 3\,yz\!)^{\,2}+ 3(\!y^{\,2}- z^{\,2}- xy- yz+ 2\,zx\!)^{\,2}}{4}$$