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