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.
On
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.
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.
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.
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}}$$
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.
On
Using Macaulay2 (version 1.16) with package SumsOfSquares:
i1 : needsPackage( "SumsOfSquares" );
i2 : R = QQ[x,y,z];
i3 : f = (x^2 + y^2 + z^2)^2 - 3*( x^3 * y + y^3 * z + z^3 * x);
i4 : sosPoly solveSOS f
1 2 1 1 1 2 2 9 1 2 2 2 1 2 2
o4 = (3)(- -x + x*y - -x*z - -y*z + -z ) + (-)(- -x + -y + x*z - y*z - -z )
2 2 2 2 4 3 3 3
o4 : SOSPoly
Printing:
i5 : tex o4
o5 = $\texttt{SOSPoly}\left\{\texttt{coefficients}\,\Rightarrow\,\left\{3,\,\frac{9}{4}\right\},\,\texttt{generators}\,\Rightarrow\,\left\{-\frac{1}{2}\,x^{2}+x\,y-\frac{1}{2}\,x\,z
-\frac{1}{2}\,y\,z+\frac{1}{2}\,z^{2},\,-\frac{1}{3}\,x^{2}+\frac{2}{3}\,y^{2}+x\,z-y\,z-\frac{1}{3}\,z^{2}\right\},\,\texttt{ring}\,\Rightarrow\,R\right\}$
In $\TeX$:
$$\texttt{SOSPoly}\left\{\texttt{coefficients}\,\Rightarrow\,\left\{3,\,\frac{9}{4}\right\},\,\texttt{generators}\,\Rightarrow\,\left\{-\frac{1}{2}\,x^{2}+x\,y-\frac{1}{2}\,x\,z -\frac{1}{2}\,y\,z+\frac{1}{2}\,z^{2},\,-\frac{1}{3}\,x^{2}+\frac{2}{3}\,y^{2}+x\,z-y\,z-\frac{1}{3}\,z^{2}\right\},\,\texttt{ring}\,\Rightarrow\,R\right\}$$
Hence,
$$\boxed{ f = 3 \left( -\frac{1}{2}\,x^{2}+x\,y-\frac{1}{2}\,x\,z -\frac{1}{2}\,y\,z+\frac{1}{2}\,z^{2} \right)^2 + \frac{9}{4} \left( -\frac{1}{3}\,x^{2}+\frac{2}{3}\,y^{2}+x\,z-y\,z-\frac{1}{3}\,z^{2} \right)^2 }$$
$$(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.$$