{ "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/latex": [ "$$x^{4}y^{2} + x^{2}y^{4} - 3x^{2}y^{2} + 1$$" ], "text/plain": [ "x⁴y² + x²y⁴ - 3x²y² + 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://jump.dev/JuMP.jl/dev/installation/#Getting-Solvers-1) for a list of the available choices." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using SumOfSquares, CSDP\n", "factory = CSDP.Optimizer;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using SumOfSquares, MosekTools\n", "factory = Mosek.Optimizer;" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$(1)x^{4}y^{2} + (1)x^{2}y^{4} + (-3)x^{2}y^{2} + (1) \\text{ is SOS}$" ], "text/plain": [ "(1)x⁴y² + (1)x²y⁴ + (-3)x²y² + (1) is SOS" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = SOSModel(factory)\n", "@constraint(model, motzkin >= 0) # We constraint motzkin to be a sum of squares" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CSDP 6.2.0\n", "This is a pure primal feasibility problem.\n", "Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 \n", "Iter: 1 Ap: 8.21e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.2410827e+01 \n", "Iter: 2 Ap: 9.49e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.0226979e+01 \n", "Iter: 3 Ap: 8.52e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.8043991e+00 \n", "Iter: 4 Ap: 4.43e-01 Pobj: 0.0000000e+00 Ad: 9.70e-01 Dobj: -2.0131592e+01 \n", "Declaring primal infeasibility.\n" ] } ], "source": [ "optimize!(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the problem is detected as infeasible..." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "INFEASIBLE::TerminationStatusCode = 2" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "termination_status(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and that the dual solution is a certificate of the infeasibility of the problem." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dual_status(model)" ] }, { "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": 7, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$(1)x^{6}y^{2} + (2)x^{4}y^{4} + (1)x^{2}y^{6} + (-3)x^{4}y^{2} + (-3)x^{2}y^{4} + (1)x^{2} + (1)y^{2} \\text{ is SOS}$" ], "text/plain": [ "(1)x⁶y² + (2)x⁴y⁴ + (1)x²y⁶ + (-3)x⁴y² + (-3)x²y⁴ + (1)x² + (1)y² is SOS" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = SOSModel(factory)\n", "@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the (x^2 + y^2) * motzkin to be a sum of squares" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CSDP 6.2.0\n", "This is a pure primal feasibility problem.\n", "Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 \n", "Iter: 1 Ap: 7.84e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 6.1906404e+01 \n", "Iter: 2 Ap: 9.45e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 6.0656534e+01 \n", "Iter: 3 Ap: 9.05e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 2.7707890e+01 \n", "Iter: 4 Ap: 7.75e-01 Pobj: 0.0000000e+00 Ad: 8.23e-01 Dobj: 5.8982728e+00 \n", "Iter: 5 Ap: 8.15e-01 Pobj: 0.0000000e+00 Ad: 9.46e-01 Dobj: 1.3960890e+00 \n", "Iter: 6 Ap: 7.61e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.0730107e-01 \n", "Iter: 7 Ap: 7.54e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 8.8693718e-02 \n", "Iter: 8 Ap: 7.91e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.7142486e-02 \n", "Iter: 9 Ap: 7.58e-01 Pobj: 0.0000000e+00 Ad: 9.68e-01 Dobj: 8.4582270e-03 \n", "Iter: 10 Ap: 7.89e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 4.0179083e-03 \n", "Iter: 11 Ap: 7.66e-01 Pobj: 0.0000000e+00 Ad: 9.92e-01 Dobj: 7.5288196e-04 \n", "Iter: 12 Ap: 7.96e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 3.7224432e-04 \n", "Iter: 13 Ap: 7.62e-01 Pobj: 0.0000000e+00 Ad: 9.73e-01 Dobj: 6.9932371e-05 \n", "Iter: 14 Ap: 6.52e-01 Pobj: 0.0000000e+00 Ad: 6.07e-01 Dobj: 2.9520639e-05 \n" ] } ], "source": [ "optimize!(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the problem is declared feasible by the solver..." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "OPTIMAL::TerminationStatusCode = 1" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "termination_status(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FEASIBLE_POINT::ResultStatusCode = 1" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "primal_status(model)" ] }, { "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": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element MonomialVector{true}:\n", " x²\n", " xy\n", " y²\n", " x \n", " y \n", " 1 " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = SOSModel(factory)\n", "X = monomials([x, y], 0:2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a quadratic polynomial that is not necessarily a sum of squares since this is implied by the next constraint: deno >= 1." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$(noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname)$$" ], "text/plain": [ "(noname)x² + (noname)xy + (noname)y² + (noname)x + (noname)y + (noname)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@variable(model, deno, Poly(X))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$(noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname - 1) \\text{ is SOS}$" ], "text/plain": [ "(noname)x² + (noname)xy + (noname)y² + (noname)x + (noname)y + (noname - 1) is SOS" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(model, deno >= 1)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$(noname)x^{6}y^{2} + (noname)x^{5}y^{3} + (noname + noname)x^{4}y^{4} + (noname)x^{3}y^{5} + (noname)x^{2}y^{6} + (noname)x^{5}y^{2} + (noname)x^{4}y^{3} + (noname)x^{3}y^{4} + (noname)x^{2}y^{5} + (-3 noname + noname)x^{4}y^{2} + (-3 noname)x^{3}y^{3} + (-3 noname + noname)x^{2}y^{4} + (-3 noname)x^{3}y^{2} + (-3 noname)x^{2}y^{3} + (-3 noname)x^{2}y^{2} + (noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname) \\text{ is SOS}$" ], "text/plain": [ "(noname)x⁶y² + (noname)x⁵y³ + (noname + noname)x⁴y⁴ + (noname)x³y⁵ + (noname)x²y⁶ + (noname)x⁵y² + (noname)x⁴y³ + (noname)x³y⁴ + (noname)x²y⁵ + (-3 noname + noname)x⁴y² + (-3 noname)x³y³ + (-3 noname + noname)x²y⁴ + (-3 noname)x³y² + (-3 noname)x²y³ + (-3 noname)x²y² + (noname)x² + (noname)xy + (noname)y² + (noname)x + (noname)y + (noname) is SOS" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(model, deno * motzkin >= 0)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CSDP 6.2.0\n", "This is a pure primal feasibility problem.\n", "Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 \n", "Iter: 1 Ap: 7.95e-01 Pobj: 0.0000000e+00 Ad: 7.98e-01 Dobj: -1.4969375e+00 \n", "Iter: 2 Ap: 8.28e-01 Pobj: 0.0000000e+00 Ad: 7.89e-01 Dobj: 3.6682062e-01 \n", "Iter: 3 Ap: 8.25e-01 Pobj: 0.0000000e+00 Ad: 8.48e-01 Dobj: -4.2971015e-03 \n", "Iter: 4 Ap: 7.85e-01 Pobj: 0.0000000e+00 Ad: 8.48e-01 Dobj: 1.9519502e-04 \n", "Iter: 5 Ap: 7.02e-01 Pobj: 0.0000000e+00 Ad: 8.45e-01 Dobj: -1.7443054e-04 \n", "Iter: 6 Ap: 6.84e-01 Pobj: 0.0000000e+00 Ad: 7.45e-01 Dobj: 3.4354877e-05 \n", "Iter: 7 Ap: 6.59e-01 Pobj: 0.0000000e+00 Ad: 7.05e-01 Dobj: -1.2384593e-05 \n", "Iter: 8 Ap: 7.09e-01 Pobj: 0.0000000e+00 Ad: 5.96e-01 Dobj: -8.3434950e-06 \n", "Iter: 9 Ap: 6.92e-01 Pobj: 0.0000000e+00 Ad: 7.86e-01 Dobj: 1.3482976e-06 \n", "Iter: 10 Ap: 7.08e-01 Pobj: 0.0000000e+00 Ad: 6.19e-01 Dobj: -2.1547496e-07 \n", "Iter: 11 Ap: 6.29e-01 Pobj: 0.0000000e+00 Ad: 8.05e-01 Dobj: 1.1922145e-07 \n", "Iter: 12 Ap: 7.29e-01 Pobj: 0.0000000e+00 Ad: 6.16e-01 Dobj: -1.2524057e-08 \n", "Iter: 13 Ap: 8.49e-01 Pobj: 0.0000000e+00 Ad: 8.38e-01 Dobj: 8.7355197e-09 \n", "Iter: 14 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.32e-01 Dobj: -1.3938737e-09 \n", "Iter: 15 Ap: 8.46e-01 Pobj: 0.0000000e+00 Ad: 8.14e-01 Dobj: 5.3843217e-10 \n", "Iter: 16 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.60e-01 Dobj: -5.6655207e-11 \n", "Iter: 17 Ap: 9.05e-01 Pobj: 0.0000000e+00 Ad: 8.60e-01 Dobj: 1.1193794e-11 \n", "Iter: 18 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 7.25e-01 Dobj: -2.4333986e-12 \n", "Iter: 19 Ap: 6.67e-01 Pobj: 0.0000000e+00 Ad: 7.60e-01 Dobj: 9.9750431e-13 \n", "Iter: 20 Ap: 3.18e-02 Pobj: 0.0000000e+00 Ad: 5.86e-01 Dobj: -3.0148476e-13 \n", "Iter: 21 Ap: 8.80e-02 Pobj: 0.0000000e+00 Ad: 3.28e-01 Dobj: 3.1504853e-13 \n", "Iter: 22 Ap: 3.34e-03 Pobj: 0.0000000e+00 Ad: 3.27e-01 Dobj: -6.7360709e-14 \n", "Iter: 23 Ap: 1.21e-03 Pobj: 0.0000000e+00 Ad: 4.07e-01 Dobj: 1.4309873e-13 \n", "Iter: 24 Ap: 5.29e-04 Pobj: 0.0000000e+00 Ad: 1.98e-01 Dobj: -9.4979405e-14 \n", "Iter: 25 Ap: 2.24e-02 Pobj: 0.0000000e+00 Ad: 1.06e-01 Dobj: 1.0185252e-13 \n", "Iter: 26 Ap: 1.74e-02 Pobj: 0.0000000e+00 Ad: 1.20e-01 Dobj: -6.4880959e-14 \n", "Iter: 27 Ap: 2.36e-02 Pobj: 0.0000000e+00 Ad: 1.49e-01 Dobj: 8.2554804e-14 \n", "Iter: 28 Ap: 2.93e-03 Pobj: 0.0000000e+00 Ad: 1.65e-01 Dobj: -4.8411558e-14 \n", "Iter: 29 Ap: 8.78e-02 Pobj: 0.0000000e+00 Ad: 1.78e-01 Dobj: 4.5825457e-14 \n", "Iter: 30 Ap: 1.32e-03 Pobj: 0.0000000e+00 Ad: 3.17e-01 Dobj: -2.4389675e-14 \n", "Iter: 31 Ap: 4.83e-02 Pobj: 0.0000000e+00 Ad: 2.09e-01 Dobj: 3.0447835e-14 \n", "Iter: 32 Ap: 6.47e-03 Pobj: 0.0000000e+00 Ad: 2.76e-01 Dobj: -1.2564183e-14 \n", "Iter: 33 Ap: 2.78e-03 Pobj: 0.0000000e+00 Ad: 1.90e-01 Dobj: 4.5354455e-15 \n", "Iter: 34 Ap: 1.54e-05 Pobj: 0.0000000e+00 Ad: 7.77e-02 Dobj: 2.5022451e-16 \n", "Stuck at edge of primal feasibility, giving up. \n" ] } ], "source": [ "optimize!(model)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "ALMOST_OPTIMAL::TerminationStatusCode = 7" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "termination_status(model)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "NEARLY_FEASIBLE_POINT::ResultStatusCode = 2" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "primal_status(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check the denominator found by the program using JuMP.getvalue" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$11788.376578732263x^{2} + 11788.412205445507y^{2} + 20644.79368009344$$" ], "text/plain": [ "11788.376578732263x² + 11788.412205445507y² + 20644.79368009344" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "value(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": 19, "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": 20, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "plot(range(-2, stop=2, length=100), range(-2, stop=2, length=100), motzkin,\n", " 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 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" } }, "nbformat": 4, "nbformat_minor": 2 }