{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean),\n",
"$$ \\frac{x^4y^2 + x^2y^4 + 1}{3} \\ge \\sqrt[3]{x^4y^2 \\cdot x^2y^4 \\cdot 1} = x^2y^2 $$\n",
"hence\n",
"$$ x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \\ge 0. $$\n",
"The code belows construct the Motzkin polynomial using [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"x^4y^2 + x^2y^4 - 3x^2y^2 + 1"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DynamicPolynomials\n",
"@polyvar x y\n",
"motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows.\n",
"We first need to pick an SDP solver, see [here](http://www.juliaopt.org/) for a list of the available choices."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using CSDP\n",
"solver = CSDPSolver();"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"using Mosek\n",
"solver = MosekSolver(LOG=0);"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[1m\u001b[33mWarning: \u001b[39m\u001b[22m\u001b[33mNot solved to optimality, status: Infeasible\u001b[39m\n"
]
},
{
"data": {
"text/plain": [
":Infeasible"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SumOfSquares\n",
"using JuMP\n",
"m = SOSModel(solver = solver)\n",
"@constraint m motzkin >= 0 # We constraint `motzkin` to be a sum of squares\n",
"solve(m) # Returns the status `:Infeasible`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares.\n",
"Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares.\n",
"The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$.\n",
"This can be verified numerically as follows:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
":Optimal"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SumOfSquares\n",
"using JuMP\n",
"m = SOSModel(solver = solver)\n",
"@constraint m (x^2 + y^2) * motzkin >= 0 # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares\n",
"solve(m) # Returns the status `:Optimal` which means that it is feasible"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator.\n",
"In fact, the search for the denominator can be carried out in parallel to the search of the numerator.\n",
"In the example below, we search for a denominator with monomials of degrees from 0 to 2.\n",
"If none is found, we can increase the maximum degree 2 to 4, 6, 8, ...\n",
"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.\n",
"In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
":Optimal"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SumOfSquares\n",
"using JuMP\n",
"using MultivariatePolynomials\n",
"m = SOSModel(solver = solver)\n",
"X = monomials([x, y], 0:2)\n",
"# We create a quadratic polynomial that is not necessarily a sum of squares\n",
"# since this is implied by the next constraint: `deno >= 1`\n",
"@variable m deno Poly(X)\n",
"# We want the denominator polynomial to be strictly positive,\n",
"# this prevents the trivial solution deno = 0 for instance.\n",
"@constraint m deno >= 1\n",
"@constraint m deno * motzkin >= 0\n",
"solve(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check the denominator found by the program using `JuMP.getvalue`"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.8994524919313149x^2 - 8.417376223825856e-11xy + 0.8994524979159756y^2 + 6.987367755126592e-16x - 1.0014392160847468e-15y + 1.9999999943329119"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"getvalue(deno)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because a picture is worth a thousand words let's plot the beast.\n",
"We can easily extend `Plots` by adding a recipe to plot bivariate polynomials."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using RecipesBase\n",
"using MultivariatePolynomials\n",
"@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)\n",
" x, y, (x, y) -> p(variables(p) => [x, y])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"plot(linspace(-2, 2, 100), linspace(-2, 2, 100), motzkin, st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.4",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}