Adapted from: (3.6) and (3.19) of [BPT12]

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

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

\[ 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.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
@constraint(model, motzkin >= 0) # We constraint `motzkin` to be a sum of squares


We see that the problem is detected as infeasible...

INFEASIBLE::TerminationStatusCode = 2

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


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(solver)
@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares


Now the problem is declared feasible by the solver...

OPTIMAL::TerminationStatusCode = 1

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

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(solver)
X = monomials([x, y], 0:2)
6-element DynamicPolynomials.MonomialVector{true}:

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))

\[ (_[1])x^{2} + (_[2])xy + (_[3])y^{2} + (_[4])x + (_[5])y + (_[6]) \]

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)


FEASIBLE_POINT::ResultStatusCode = 1

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


\[ 12.292628687208005x^{2} + 12.292628753109426y^{2} + 23.251596478120113 \]

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
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
import Plots
    range(-2, stop=2, length=100),
    range(-2, stop=2, length=100),
    st = [:surface],
    clims = (-10, 80)

This page was generated using Literate.jl.