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^4y^2 + x^2y^4 - 3x^2y^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 CSDP
solver = CSDPSolver();
using Mosek
solver = MosekSolver();
using SumOfSquares
using JuMP
m = SOSModel(solver = solver)
@constraint m motzkin >= 0 # We constraint `motzkin` to be a sum of squares
solve(m) # Returns the status `:Infeasible`
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 22 Cones : 0 Scalar variables : 0 Matrix variables : 1 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator - tries : 0 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 22 Cones : 0 Scalar variables : 0 Matrix variables : 1 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 22 Optimizer - Cones : 0 Optimizer - Scalar variables : 0 conic : 0 Optimizer - Semi-definite variables: 1 scalarized : 36 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 253 after factor : 253 Factor - dense dim. : 0 flops : 7.61e+03 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.5e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.02 1 2.7e-01 1.1e-01 2.7e-01 1.21e+00 0.000000000e+00 2.126988719e-01 1.1e-01 0.02 2 4.9e-02 1.9e-02 3.9e-02 2.29e-01 0.000000000e+00 7.994925687e-01 1.9e-02 0.02 3 1.1e-04 4.3e-05 7.3e-07 -7.85e-01 0.000000000e+00 1.360221226e+04 4.3e-05 0.02 4 8.5e-12 3.0e-13 3.0e-13 -9.99e-01 0.000000000e+00 3.463598352e-01 3.0e-13 0.02 Optimizer terminated. Time: 0.03 Interior-point solution summary Problem status : PRIMAL_INFEASIBLE Solution status : PRIMAL_INFEASIBLE_CER Dual. obj: 3.4635983521e-01 nrm: 5e+00 Viol. con: 0e+00 barvar: 6e-13
WARNING: Not solved to optimality, status: Infeasible
:Infeasible
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:
using SumOfSquares
using JuMP
m = SOSModel(solver = solver)
@constraint m (x^2 + y^2) * motzkin >= 0 # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares
solve(m) # Returns the status `:Optimal` which means that it is feasible
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 36 Cones : 0 Scalar variables : 0 Matrix variables : 1 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator - tries : 0 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 36 Cones : 0 Scalar variables : 0 Matrix variables : 1 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 36 Optimizer - Cones : 0 Optimizer - Scalar variables : 0 conic : 0 Optimizer - Semi-definite variables: 1 scalarized : 78 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 666 after factor : 666 Factor - dense dim. : 0 flops : 3.20e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.5e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 1 4.0e-01 1.6e-01 4.2e-01 9.52e-01 0.000000000e+00 -4.673689580e-02 1.6e-01 0.00 2 6.6e-02 2.6e-02 1.4e-01 8.63e-01 0.000000000e+00 1.870385636e-02 2.6e-02 0.00 3 1.1e-02 4.5e-03 5.4e-02 9.17e-01 0.000000000e+00 6.770665391e-03 4.5e-03 0.00 4 2.8e-03 1.1e-03 2.7e-02 9.36e-01 0.000000000e+00 1.550400506e-03 1.1e-03 0.00 5 7.3e-04 2.9e-04 1.9e-02 1.27e+00 0.000000000e+00 -1.939464728e-04 2.9e-04 0.00 6 2.4e-04 9.6e-05 6.7e-03 5.73e-01 0.000000000e+00 2.366945923e-04 9.6e-05 0.00 7 4.8e-05 1.9e-05 7.2e-03 1.51e+00 0.000000000e+00 -5.200009967e-05 1.9e-05 0.00 8 9.8e-06 3.9e-06 3.0e-03 1.07e+00 0.000000000e+00 -9.294696566e-06 3.9e-06 0.01 9 2.6e-06 1.0e-06 8.7e-04 6.61e-01 0.000000000e+00 -3.670843815e-07 1.0e-06 0.01 10 4.6e-07 1.8e-07 6.7e-04 1.47e+00 0.000000000e+00 -4.560298380e-07 1.8e-07 0.01 11 1.0e-07 4.2e-08 2.0e-04 7.80e-01 0.000000000e+00 -3.923684580e-08 4.2e-08 0.01 12 2.3e-08 9.1e-09 1.2e-04 1.26e+00 0.000000000e+00 -1.444374153e-08 9.0e-09 0.01 Optimizer terminated. Time: 0.01 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: 0.0000000000e+00 nrm: 4e+00 Viol. con: 5e-08 barvar: 0e+00 Dual. obj: -1.4443741647e-08 nrm: 4e+00 Viol. con: 0e+00 barvar: 2e-08
:Optimal
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.
using SumOfSquares
using PolyJuMP
using JuMP
m = SOSModel(solver = solver)
X = monomials([x, y], 0:2)
# We create a quadratic polynomial that is not necessarily a sum of squares
# since this is implied by the next constraint: `deno >= 1`
@variable m deno Poly(X)
# We want the denominator polynomial to be strictly positive,
# this prevents the trivial solution deno = 0 for instance.
@constraint m deno >= 1
@constraint m deno * motzkin >= 0
solve(m)
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 45 Cones : 0 Scalar variables : 6 Matrix variables : 2 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator - tries : 0 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 45 Cones : 0 Scalar variables : 6 Matrix variables : 2 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 45 Optimizer - Cones : 1 Optimizer - Scalar variables : 7 conic : 7 Optimizer - Semi-definite variables: 2 scalarized : 97 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 927 after factor : 927 Factor - dense dim. : 0 flops : 4.64e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 1 4.0e-01 2.0e-01 2.1e-01 -3.77e-01 0.000000000e+00 5.278524367e-01 2.0e-01 0.00 2 5.1e-02 2.5e-02 1.2e-01 8.59e-01 0.000000000e+00 -2.053606671e-04 2.5e-02 0.00 3 7.8e-03 3.9e-03 4.5e-02 1.07e+00 0.000000000e+00 3.927299184e-04 3.9e-03 0.00 4 1.9e-03 9.6e-04 2.0e-02 9.60e-01 0.000000000e+00 4.160764371e-04 9.6e-04 0.00 5 3.5e-04 1.7e-04 8.8e-03 1.09e+00 0.000000000e+00 8.343664284e-05 1.7e-04 0.00 6 1.0e-04 5.2e-05 4.3e-03 9.48e-01 0.000000000e+00 5.014506528e-05 5.2e-05 0.01 7 1.9e-05 9.4e-06 2.0e-03 1.08e+00 0.000000000e+00 6.674228176e-06 9.4e-06 0.01 8 5.8e-06 2.9e-06 9.9e-04 9.61e-01 0.000000000e+00 3.271958410e-06 2.9e-06 0.01 9 1.1e-06 5.3e-07 4.6e-04 1.08e+00 0.000000000e+00 4.169073231e-07 5.3e-07 0.01 10 3.3e-07 1.6e-07 2.4e-04 9.62e-01 0.000000000e+00 1.967450815e-07 1.6e-07 0.01 11 6.0e-08 3.0e-08 1.1e-04 1.08e+00 0.000000000e+00 2.395130544e-08 3.0e-08 0.01 12 1.9e-08 9.4e-09 5.6e-05 9.62e-01 0.000000000e+00 1.121624481e-08 9.3e-09 0.01 13 3.4e-09 1.8e-09 2.6e-05 1.08e+00 0.000000000e+00 1.356679134e-09 1.7e-09 0.01 Optimizer terminated. Time: 0.01 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: 0.0000000000e+00 nrm: 2e+00 Viol. con: 6e-09 var: 0e+00 barvar: 0e+00 Dual. obj: 1.3566791338e-09 nrm: 3e+00 Viol. con: 0e+00 var: 2e-09 barvar: 3e-09
:Optimal
We can check the denominator found by the program using JuMP.getvalue
getvalue(deno)
0.8994524919313149x^2 - 8.417376223825856e-11xy + 0.8994524979159756y^2 + 6.987367755126592e-16x - 1.0014392160847468e-15y + 1.9999999943329119
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(linspace(-2, 2, 100), linspace(-2, 2, 100), motzkin, st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))