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 [10]:
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2
Out[10]:
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.

In [42]:
using CSDP
solver = CSDPSolver();
In [43]:
using Mosek
solver = MosekSolver();
In [44]:
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
Out[44]:
: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:

In [45]:
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  
Out[45]:
: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.

In [46]:
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  
Out[46]:
:Optimal

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

In [47]:
getvalue(deno)
Out[47]:
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.

In [48]:
using RecipesBase
using MultivariatePolynomials
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
end
In [49]:
using Plots
plot(linspace(-2, 2, 100), linspace(-2, 2, 100), motzkin, st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))
Out[49]:
In [ ]: