{ "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": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x^4y^2 + x^2y^4 - 3x^2y^2 + 1" ] }, "execution_count": 10, "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": 42, "metadata": {}, "outputs": [], "source": [ "using CSDP\n", "solver = CSDPSolver();" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using Mosek\n", "solver = MosekSolver();" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 22 \n", " Cones : 0 \n", " Scalar variables : 0 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator - tries : 0 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 22 \n", " Cones : 0 \n", " Scalar variables : 0 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 4 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 22\n", "Optimizer - Cones : 0\n", "Optimizer - Scalar variables : 0 conic : 0 \n", "Optimizer - Semi-definite variables: 1 scalarized : 36 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 253 after factor : 253 \n", "Factor - dense dim. : 0 flops : 7.61e+03 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 2.5e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.02 \n", "1 2.7e-01 1.1e-01 2.7e-01 1.21e+00 0.000000000e+00 2.126988719e-01 1.1e-01 0.02 \n", "2 4.9e-02 1.9e-02 3.9e-02 2.29e-01 0.000000000e+00 7.994925687e-01 1.9e-02 0.02 \n", "3 1.1e-04 4.3e-05 7.3e-07 -7.85e-01 0.000000000e+00 1.360221226e+04 4.3e-05 0.02 \n", "4 8.5e-12 3.0e-13 3.0e-13 -9.99e-01 0.000000000e+00 3.463598352e-01 3.0e-13 0.02 \n", "Optimizer terminated. Time: 0.03 \n", "\n", "\n", "Interior-point solution summary\n", " Problem status : PRIMAL_INFEASIBLE\n", " Solution status : PRIMAL_INFEASIBLE_CER\n", " Dual. obj: 3.4635983521e-01 nrm: 5e+00 Viol. con: 0e+00 barvar: 6e-13 \n" ] }, { "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": 44, "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": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 36 \n", " Cones : 0 \n", " Scalar variables : 0 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator - tries : 0 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 36 \n", " Cones : 0 \n", " Scalar variables : 0 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 4 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 36\n", "Optimizer - Cones : 0\n", "Optimizer - Scalar variables : 0 conic : 0 \n", "Optimizer - Semi-definite variables: 1 scalarized : 78 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 666 after factor : 666 \n", "Factor - dense dim. : 0 flops : 3.20e+04 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 2.5e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 \n", "1 4.0e-01 1.6e-01 4.2e-01 9.52e-01 0.000000000e+00 -4.673689580e-02 1.6e-01 0.00 \n", "2 6.6e-02 2.6e-02 1.4e-01 8.63e-01 0.000000000e+00 1.870385636e-02 2.6e-02 0.00 \n", "3 1.1e-02 4.5e-03 5.4e-02 9.17e-01 0.000000000e+00 6.770665391e-03 4.5e-03 0.00 \n", "4 2.8e-03 1.1e-03 2.7e-02 9.36e-01 0.000000000e+00 1.550400506e-03 1.1e-03 0.00 \n", "5 7.3e-04 2.9e-04 1.9e-02 1.27e+00 0.000000000e+00 -1.939464728e-04 2.9e-04 0.00 \n", "6 2.4e-04 9.6e-05 6.7e-03 5.73e-01 0.000000000e+00 2.366945923e-04 9.6e-05 0.00 \n", "7 4.8e-05 1.9e-05 7.2e-03 1.51e+00 0.000000000e+00 -5.200009967e-05 1.9e-05 0.00 \n", "8 9.8e-06 3.9e-06 3.0e-03 1.07e+00 0.000000000e+00 -9.294696566e-06 3.9e-06 0.01 \n", "9 2.6e-06 1.0e-06 8.7e-04 6.61e-01 0.000000000e+00 -3.670843815e-07 1.0e-06 0.01 \n", "10 4.6e-07 1.8e-07 6.7e-04 1.47e+00 0.000000000e+00 -4.560298380e-07 1.8e-07 0.01 \n", "11 1.0e-07 4.2e-08 2.0e-04 7.80e-01 0.000000000e+00 -3.923684580e-08 4.2e-08 0.01 \n", "12 2.3e-08 9.1e-09 1.2e-04 1.26e+00 0.000000000e+00 -1.444374153e-08 9.0e-09 0.01 \n", "Optimizer terminated. Time: 0.01 \n", "\n", "\n", "Interior-point solution summary\n", " Problem status : PRIMAL_AND_DUAL_FEASIBLE\n", " Solution status : OPTIMAL\n", " Primal. obj: 0.0000000000e+00 nrm: 4e+00 Viol. con: 5e-08 barvar: 0e+00 \n", " Dual. obj: -1.4443741647e-08 nrm: 4e+00 Viol. con: 0e+00 barvar: 2e-08 \n" ] }, { "data": { "text/plain": [ ":Optimal" ] }, "execution_count": 45, "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": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 45 \n", " Cones : 0 \n", " Scalar variables : 6 \n", " Matrix variables : 2 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator - tries : 0 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 45 \n", " Cones : 0 \n", " Scalar variables : 6 \n", " Matrix variables : 2 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 4 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 45\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 7 conic : 7 \n", "Optimizer - Semi-definite variables: 2 scalarized : 97 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 927 after factor : 927 \n", "Factor - dense dim. : 0 flops : 4.64e+04 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 2.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 \n", "1 4.0e-01 2.0e-01 2.1e-01 -3.77e-01 0.000000000e+00 5.278524367e-01 2.0e-01 0.00 \n", "2 5.1e-02 2.5e-02 1.2e-01 8.59e-01 0.000000000e+00 -2.053606671e-04 2.5e-02 0.00 \n", "3 7.8e-03 3.9e-03 4.5e-02 1.07e+00 0.000000000e+00 3.927299184e-04 3.9e-03 0.00 \n", "4 1.9e-03 9.6e-04 2.0e-02 9.60e-01 0.000000000e+00 4.160764371e-04 9.6e-04 0.00 \n", "5 3.5e-04 1.7e-04 8.8e-03 1.09e+00 0.000000000e+00 8.343664284e-05 1.7e-04 0.00 \n", "6 1.0e-04 5.2e-05 4.3e-03 9.48e-01 0.000000000e+00 5.014506528e-05 5.2e-05 0.01 \n", "7 1.9e-05 9.4e-06 2.0e-03 1.08e+00 0.000000000e+00 6.674228176e-06 9.4e-06 0.01 \n", "8 5.8e-06 2.9e-06 9.9e-04 9.61e-01 0.000000000e+00 3.271958410e-06 2.9e-06 0.01 \n", "9 1.1e-06 5.3e-07 4.6e-04 1.08e+00 0.000000000e+00 4.169073231e-07 5.3e-07 0.01 \n", "10 3.3e-07 1.6e-07 2.4e-04 9.62e-01 0.000000000e+00 1.967450815e-07 1.6e-07 0.01 \n", "11 6.0e-08 3.0e-08 1.1e-04 1.08e+00 0.000000000e+00 2.395130544e-08 3.0e-08 0.01 \n", "12 1.9e-08 9.4e-09 5.6e-05 9.62e-01 0.000000000e+00 1.121624481e-08 9.3e-09 0.01 \n", "13 3.4e-09 1.8e-09 2.6e-05 1.08e+00 0.000000000e+00 1.356679134e-09 1.7e-09 0.01 \n", "Optimizer terminated. Time: 0.01 \n", "\n", "\n", "Interior-point solution summary\n", " Problem status : PRIMAL_AND_DUAL_FEASIBLE\n", " Solution status : OPTIMAL\n", " Primal. obj: 0.0000000000e+00 nrm: 2e+00 Viol. con: 6e-09 var: 0e+00 barvar: 0e+00 \n", " Dual. obj: 1.3566791338e-09 nrm: 3e+00 Viol. con: 0e+00 var: 2e-09 barvar: 3e-09 \n" ] }, { "data": { "text/plain": [ ":Optimal" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using SumOfSquares\n", "using PolyJuMP\n", "using JuMP\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": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.8994524919313149x^2 - 8.417376223825856e-11xy + 0.8994524979159756y^2 + 6.987367755126592e-16x - 1.0014392160847468e-15y + 1.9999999943329119" ] }, "execution_count": 47, "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": 48, "metadata": { "collapsed": true }, "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": 49, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "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": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }