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.

In [1]:
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2
Out[1]:
$$ x^{4}y^{2} + x^{2}y^{4} - 3x^{2}y^{2} + 1 $$

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.

In [2]:
using SumOfSquares, CSDP
factory = CSDP.Optimizer;
In [ ]:
using SumOfSquares, MosekTools
factory = Mosek.Optimizer;
In [3]:
model = SOSModel(factory)
@constraint(model, motzkin >= 0) # We constraint `motzkin` to be a sum of squares
Out[3]:
$ (1)x^{4}y^{2} + (1)x^{2}y^{4} + (-3)x^{2}y^{2} + (1) \text{ is SOS} $
In [4]:
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...

In [5]:
termination_status(model)
Out[5]:
INFEASIBLE::TerminationStatusCode = 2

... and that the dual solution is a certificate of the infeasibility of the problem.

In [6]:
dual_status(model)
Out[6]:
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:

In [7]:
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
Out[7]:
$ (1)x^{6}y^{2} + (2)x^{4}y^{4} + (1)x^{2}y^{6} + (-3)x^{4}y^{2} + (-3)x^{2}y^{4} + (1)x^{2} + (1)y^{2} \text{ is SOS} $
In [8]:
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...

In [9]:
termination_status(model)
Out[9]:
OPTIMAL::TerminationStatusCode = 1

... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial.

In [10]:
primal_status(model)
Out[10]:
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.

In [11]:
model = SOSModel(factory)
X = monomials([x, y], 0:2)
Out[11]:
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.

In [12]:
@variable(model, deno, Poly(X))
Out[12]:
$$ (noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname) $$

We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance.

In [13]:
@constraint(model, deno >= 1)
Out[13]:
$ (noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname - 1) \text{ is SOS} $
In [14]:
@constraint(model, deno * motzkin >= 0)
Out[14]:
$ (noname)x^{6}y^{2} + (noname)x^{5}y^{3} + (noname + noname)x^{4}y^{4} + (noname)x^{3}y^{5} + (noname)x^{2}y^{6} + (noname)x^{5}y^{2} + (noname)x^{4}y^{3} + (noname)x^{3}y^{4} + (noname)x^{2}y^{5} + (-3 noname + noname)x^{4}y^{2} + (-3 noname)x^{3}y^{3} + (-3 noname + noname)x^{2}y^{4} + (-3 noname)x^{3}y^{2} + (-3 noname)x^{2}y^{3} + (-3 noname)x^{2}y^{2} + (noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname) \text{ is SOS} $
In [15]:
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. 
In [16]:
termination_status(model)
Out[16]:
ALMOST_OPTIMAL::TerminationStatusCode = 7
In [17]:
primal_status(model)
Out[17]:
NEARLY_FEASIBLE_POINT::ResultStatusCode = 2

We can check the denominator found by the program using JuMP.getvalue

In [18]:
value(deno)
Out[18]:
$$ 11788.376578732263x^{2} + 11788.412205445507y^{2} + 20644.79368009344 $$

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.

In [19]:
using RecipesBase
using MultivariatePolynomials
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
end
In [20]:
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))
Out[20]:
In [ ]: