The first explicit example of nonnegative polynomial that is not a sum of squares was found by Motzkin in 1967. By the Arithmetic-geometric mean, $$ \frac{x^4y^2 + x^2y^4 + 1}{3} \ge \sqrt[3]{x^4y^2 \cdot x^2y^4 \cdot 1} = x^2y^2 $$ hence $$ x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \ge 0. $$ The code belows construct the Motzkin polynomial using DynamicPolynomials.
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2
The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows. We first need to pick an SDP solver, see here for a list of the available choices.
using SumOfSquares, CSDP
factory = CSDP.Optimizer;
using SumOfSquares, MosekTools
factory = Mosek.Optimizer;
model = SOSModel(factory)
@constraint(model, motzkin >= 0) # We constraint `motzkin` to be a sum of squares
optimize!(model)
CSDP 6.2.0 This is a pure primal feasibility problem. Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 8.21e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.2410827e+01 Iter: 2 Ap: 9.49e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.0226979e+01 Iter: 3 Ap: 8.52e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.8043991e+00 Iter: 4 Ap: 4.43e-01 Pobj: 0.0000000e+00 Ad: 9.70e-01 Dobj: -2.0131592e+01 Declaring primal infeasibility.
We see that the problem is detected as infeasible...
termination_status(model)
INFEASIBLE::TerminationStatusCode = 2
... and that the dual solution is a certificate of the infeasibility of the problem.
dual_status(model)
INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4
Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares. Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares. The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$. This can be verified numerically as follows:
model = SOSModel(factory)
@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares
optimize!(model)
CSDP 6.2.0 This is a pure primal feasibility problem. Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 7.84e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 6.1906404e+01 Iter: 2 Ap: 9.45e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 6.0656534e+01 Iter: 3 Ap: 9.05e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 2.7707890e+01 Iter: 4 Ap: 7.75e-01 Pobj: 0.0000000e+00 Ad: 8.23e-01 Dobj: 5.8982728e+00 Iter: 5 Ap: 8.15e-01 Pobj: 0.0000000e+00 Ad: 9.46e-01 Dobj: 1.3960890e+00 Iter: 6 Ap: 7.61e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.0730107e-01 Iter: 7 Ap: 7.54e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 8.8693718e-02 Iter: 8 Ap: 7.91e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.7142486e-02 Iter: 9 Ap: 7.58e-01 Pobj: 0.0000000e+00 Ad: 9.68e-01 Dobj: 8.4582270e-03 Iter: 10 Ap: 7.89e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.0179083e-03 Iter: 11 Ap: 7.66e-01 Pobj: 0.0000000e+00 Ad: 9.92e-01 Dobj: 7.5288196e-04 Iter: 12 Ap: 7.96e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.7224432e-04 Iter: 13 Ap: 7.62e-01 Pobj: 0.0000000e+00 Ad: 9.73e-01 Dobj: 6.9932371e-05 Iter: 14 Ap: 6.52e-01 Pobj: 0.0000000e+00 Ad: 6.07e-01 Dobj: 2.9520639e-05
Now the problem is declared feasible by the solver...
termination_status(model)
OPTIMAL::TerminationStatusCode = 1
... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial.
primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator. In fact, the search for the denominator can be carried out in parallel to the search of the numerator. In the example below, we search for a denominator with monomials of degrees from 0 to 2. If none is found, we can increase the maximum degree 2 to 4, 6, 8, ... This gives a hierarchy of programs to try in order to certify the nonnegativity of a polynomial by identifying it with a fraction of sum of squares polynomials. In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works.
model = SOSModel(factory)
X = monomials([x, y], 0:2)
6-element MonomialVector{true}: x² xy y² x y 1
We create a quadratic polynomial that is not necessarily a sum of squares since this is implied by the next constraint: deno >= 1
.
@variable(model, deno, Poly(X))
We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance.
@constraint(model, deno >= 1)
@constraint(model, deno * motzkin >= 0)
optimize!(model)
CSDP 6.2.0 This is a pure primal feasibility problem. Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 7.95e-01 Pobj: 0.0000000e+00 Ad: 7.98e-01 Dobj: -1.4969375e+00 Iter: 2 Ap: 8.28e-01 Pobj: 0.0000000e+00 Ad: 7.89e-01 Dobj: 3.6682062e-01 Iter: 3 Ap: 8.25e-01 Pobj: 0.0000000e+00 Ad: 8.48e-01 Dobj: -4.2971015e-03 Iter: 4 Ap: 7.85e-01 Pobj: 0.0000000e+00 Ad: 8.48e-01 Dobj: 1.9519502e-04 Iter: 5 Ap: 7.02e-01 Pobj: 0.0000000e+00 Ad: 8.45e-01 Dobj: -1.7443054e-04 Iter: 6 Ap: 6.84e-01 Pobj: 0.0000000e+00 Ad: 7.45e-01 Dobj: 3.4354877e-05 Iter: 7 Ap: 6.59e-01 Pobj: 0.0000000e+00 Ad: 7.05e-01 Dobj: -1.2384593e-05 Iter: 8 Ap: 7.09e-01 Pobj: 0.0000000e+00 Ad: 5.96e-01 Dobj: -8.3434950e-06 Iter: 9 Ap: 6.92e-01 Pobj: 0.0000000e+00 Ad: 7.86e-01 Dobj: 1.3482976e-06 Iter: 10 Ap: 7.08e-01 Pobj: 0.0000000e+00 Ad: 6.19e-01 Dobj: -2.1547496e-07 Iter: 11 Ap: 6.29e-01 Pobj: 0.0000000e+00 Ad: 8.05e-01 Dobj: 1.1922145e-07 Iter: 12 Ap: 7.29e-01 Pobj: 0.0000000e+00 Ad: 6.16e-01 Dobj: -1.2524057e-08 Iter: 13 Ap: 8.49e-01 Pobj: 0.0000000e+00 Ad: 8.38e-01 Dobj: 8.7355197e-09 Iter: 14 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.32e-01 Dobj: -1.3938737e-09 Iter: 15 Ap: 8.46e-01 Pobj: 0.0000000e+00 Ad: 8.14e-01 Dobj: 5.3843217e-10 Iter: 16 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.60e-01 Dobj: -5.6655207e-11 Iter: 17 Ap: 9.05e-01 Pobj: 0.0000000e+00 Ad: 8.60e-01 Dobj: 1.1193794e-11 Iter: 18 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.25e-01 Dobj: -2.4333986e-12 Iter: 19 Ap: 6.67e-01 Pobj: 0.0000000e+00 Ad: 7.60e-01 Dobj: 9.9750431e-13 Iter: 20 Ap: 3.18e-02 Pobj: 0.0000000e+00 Ad: 5.86e-01 Dobj: -3.0148476e-13 Iter: 21 Ap: 8.80e-02 Pobj: 0.0000000e+00 Ad: 3.28e-01 Dobj: 3.1504853e-13 Iter: 22 Ap: 3.34e-03 Pobj: 0.0000000e+00 Ad: 3.27e-01 Dobj: -6.7360709e-14 Iter: 23 Ap: 1.21e-03 Pobj: 0.0000000e+00 Ad: 4.07e-01 Dobj: 1.4309873e-13 Iter: 24 Ap: 5.29e-04 Pobj: 0.0000000e+00 Ad: 1.98e-01 Dobj: -9.4979405e-14 Iter: 25 Ap: 2.24e-02 Pobj: 0.0000000e+00 Ad: 1.06e-01 Dobj: 1.0185252e-13 Iter: 26 Ap: 1.74e-02 Pobj: 0.0000000e+00 Ad: 1.20e-01 Dobj: -6.4880959e-14 Iter: 27 Ap: 2.36e-02 Pobj: 0.0000000e+00 Ad: 1.49e-01 Dobj: 8.2554804e-14 Iter: 28 Ap: 2.93e-03 Pobj: 0.0000000e+00 Ad: 1.65e-01 Dobj: -4.8411558e-14 Iter: 29 Ap: 8.78e-02 Pobj: 0.0000000e+00 Ad: 1.78e-01 Dobj: 4.5825457e-14 Iter: 30 Ap: 1.32e-03 Pobj: 0.0000000e+00 Ad: 3.17e-01 Dobj: -2.4389675e-14 Iter: 31 Ap: 4.83e-02 Pobj: 0.0000000e+00 Ad: 2.09e-01 Dobj: 3.0447835e-14 Iter: 32 Ap: 6.47e-03 Pobj: 0.0000000e+00 Ad: 2.76e-01 Dobj: -1.2564183e-14 Iter: 33 Ap: 2.78e-03 Pobj: 0.0000000e+00 Ad: 1.90e-01 Dobj: 4.5354455e-15 Iter: 34 Ap: 1.54e-05 Pobj: 0.0000000e+00 Ad: 7.77e-02 Dobj: 2.5022451e-16 Stuck at edge of primal feasibility, giving up.
termination_status(model)
ALMOST_OPTIMAL::TerminationStatusCode = 7
primal_status(model)
NEARLY_FEASIBLE_POINT::ResultStatusCode = 2
We can check the denominator found by the program using JuMP.getvalue
value(deno)
Because a picture is worth a thousand words let's plot the beast.
We can easily extend Plots
by adding a recipe to plot bivariate polynomials.
using RecipesBase
using MultivariatePolynomials
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
x, y, (x, y) -> p(variables(p) => [x, y])
end
using Plots
plot(range(-2, stop=2, length=100), range(-2, stop=2, length=100), motzkin,
st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))