# Example: Multivariate polynomial optimization using symmetry reduction

In this example, we consider minimizing a multivariate polynomial with $S_3$ symmetries. The example is inspirated by example 7.1 from [7]. See the Examples folder for the file with the code. We consider the polynomial

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

for which we want to find the minimum value $f_{min} = \min_{x,y,z} f(x,y,z)$. Relaxing the problem with a sum-of-squares constraint gives

\begin{aligned} \max \quad& M &\\ \text{s.t.} \quad& f - M &\text{is a sum-of-squares,} \end{aligned}

Since the polynomial $f$ is invariant under permuting the variables (i.e., the group action of $S_3$), it is natural to consider only invariant sum-of-squares polynomials. From [7], we know that any sum-of-squares polynomial invariant under the action of $S_3$ can be written as

$$$\langle Y_1, ww^{\sf T} \rangle + \langle Y_2, \Pi_2 \otimes ww^{\sf T} \rangle + \langle Y_3, \Pi_3 \otimes ww^{\sf T} \rangle$$$

where $w(x,y,z)$ is a vector of basis polynomials of the invariant ring $\R[x,y,z]^{S_3} = \R[x+y+z, xy+yz+xz, xyz]$, and where $\Pi_2(x,y,z) = ((x-y)(y-z)(z-x))^2$. The matrix $\Pi_3$ is of rank 2 and has the decomposition $v_1 v_1^{\sf T} + v_2 v_2^{\sf T}$ with

$$$v_1 = \frac{1}{\sqrt{2}}\begin{pmatrix} 2x-y-z \\ 2yz-zx-xy \end{pmatrix} \quad \text{ and } \quad v_2 = \sqrt{\frac{3}{2}}\begin{pmatrix} y-z \\ zx-xy \end{pmatrix}$$$

Since we consider sum-of-squares polynomials of a certain degree $d$, we restrict to the elements of the matrices $\Pi_i \otimes ww^{\sf T}$ with degree at most $\lfloor d/2 \rfloor$.

To sample this constraint we need a three-variate minimal unisolvent set for invariant polynomials of degree at most $d$. In this case, one such example are the representatives of the orbits of $S_3$ of the rational points in the simplex with denominator $d$, since these points are invariant under $S_3$ and are minimal unisolvent (see [2]). However, if the polynomial space is more complicated it is unclear what a minimal unisolvent set is. To show one approach on this we instead make a grid of points which unisolvent but not minimal, and we use the approach of [2] through approximatefekete to choose a good subset of these points and a corresponding good basis for $w$.

As for the example of the Delsarte bound, the objective is simply one free variable $M$ with coefficient $1$. We also create a function to generate an invariant basis with variables $x,y,z$.

using ClusteredLowRankSolver, AbstractAlgebra

function invariant_basis(x,y,z, d)
# create a vector with a precise type
v = [(x*y*z)^0]
for deg=1:d, j=0:div(deg,3), i=0:div(deg-3j,2)
# monomials in the invariant elementary polynomials
# ordered by degree
push!(v, (x+y+z)^(deg-2i-3j) * (x*y+y*z+z*x)^i * (x*y*z)^j)
end
return v
end

function min_f(d)

obj = Objective(0, Dict(), Dict(:M => 1))

In this case, we need a three-variate polynomial ring, and a basis of invariant polynomials. We also use approximatefekete to find a subset of the sample points with a good basis.

    FF = RealField
R, (x,y,z) = polynomial_ring(FF, [:x, :y, :z])
# The polynomial f:
f =  x^4 + y^4 + z^4 - 4x*y*z + x + y + z

# An invariant basis up to degree d:
basis = invariant_basis(x, y, z, 2d)
# For the sum-of-squares polynomials we have to
# select elements of the basis based on the degree
degrees = [total_degree(p) for p in basis]

# generate samples and a good basis
cheb_points = [vcat(sample_points_chebyshev(2d+k)...) for k=0:2]
samples_grid = [[cheb_points[1][i+1], cheb_points[2][j+1], cheb_points[3][k+1]]
for i=0:2d for j=0:2d+1 for k=0:2d+2]
basis, samples = approximatefekete(basis, samples_grid)

Now we will construct the constraint matrices corresponding to the sum-of-squares parts. Although approximatefekete returns a polynomial basis only characterized by the evaluations on the sample points, basic operations on normal polynomials will also work with sampled polynomials.

    psd_dict = Dict()

equivariants = [[[R(1)]], [[(x-y)*(y-z)*(z-x)]], [[(2x-y-z), (2y*z-x*z-x*y)], [(y-z), (x*z-x*y)]]]
# The squares of the prefactors, in case we want to define the program over QQ
factors = [[1], [1], [1//2, 3//2]]
for eqi in eachindex(equivariants)
vecs = []
for r in eachindex(equivariants[eqi])
vec = []
for eq in equivariants[eqi][r], (q, qdeg) in zip(basis, degrees)
# only add terms for which the diagonal entry is of degree <=2d
if 2total_degree(eq) + 2qdeg <= 2d
push!(vec, eq * q)
end
end
if length(vec) > 0
push!(vecs, vec)
end
end
# only add variables for the non-empty matrices
if length(vecs) > 0
psd_dict[(:trivariatesos, eqi)] = LowRankMatPol(factors[eqi], vecs)
end
end

Now we can formulate the constraint and solve the problem:

    # the constraint is SOS + M = f
constr = Constraint(f, psd_dict, Dict(:M => 1), samples)
problem = Problem(Maximize(obj), [constr])

status, primalsol, dualsol, time, errorcode = solvesdp(problem)
return objvalue(problem, dualsol)
end

min_f(2)
-2.112913881423605414455367841281031182496061043438367972491981714926222891733139