{ "cells": [ { "cell_type": "markdown", "source": [ "# Polynomial Optimization" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "**Contributed by**: Benoît Legat" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Introduction" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "Consider the polynomial optimization problem [L09, Example 2.2] of\n", "minimizing the polynomial $x^3 - x^2 + 2xy - y^2 + y^3$\n", "over the polyhedron defined by the inequalities $x \\ge 0, y \\ge 0$ and $x + y \\geq 1$." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "[L09] Lasserre, J. B.\n", "*Moments, positive polynomials and their applications*.\n", "World Scientific, **2009**." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "(0, 1//4, 0)" }, "metadata": {}, "execution_count": 1 } ], "cell_type": "code", "source": [ "using DynamicPolynomials\n", "@polyvar x y\n", "p = x^3 - x^2 + 2x*y -y^2 + y^3\n", "using SumOfSquares\n", "S = @set x >= 0 && y >= 0 && x + y >= 1\n", "p(x=>1, y=>0), p(x=>1//2, y=>1//2), p(x=>0, y=>1)" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "The optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but [Ipopt](https://github.com/jump-dev/Ipopt.jl/) only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success: SDP solved\n", "Primal objective value: 0.0000000e+00 \n", "Dual objective value: 0.0000000e+00 \n", "Relative primal infeasibility: 1.53e-16 \n", "Relative dual infeasibility: 5.00e-11 \n", "Real Relative Gap: 0.00e+00 \n", "XZ Relative Gap: 1.06e-10 \n", "DIMACS error measures: 2.30e-16 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.06e-10\n", "\n", "******************************************************************************\n", "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", " For more information visit https://github.com/coin-or/Ipopt\n", "******************************************************************************\n", "\n", "termination_status(model) = MathOptInterface.LOCALLY_SOLVED\n", "value(a) = 0.4999999966705879\n", "value(b) = 0.49999999667061656\n", "objective_value(model) = 0.24999999500590336\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "0.24999999500590336" }, "metadata": {}, "execution_count": 2 } ], "cell_type": "code", "source": [ "import Ipopt\n", "model = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n", "@variable(model, a >= 0)\n", "@variable(model, b >= 0)\n", "@constraint(model, a + b >= 1)\n", "@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)\n", "optimize!(model)\n", "@show termination_status(model)\n", "@show value(a)\n", "@show value(b)\n", "@show objective_value(model)" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Note that the problem can be written equivalently as follows using [registered functions](https://jump.dev/JuMP.jl/v0.21/nlp/#User-defined-Functions-1)." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "termination_status(model) = MathOptInterface.LOCALLY_SOLVED\n", "value(a) = 0.499999995002878\n", "value(b) = 0.4999999950028773\n", "objective_value(model) = 0.24999999250431656\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "0.24999999250431656" }, "metadata": {}, "execution_count": 3 } ], "cell_type": "code", "source": [ "using Ipopt\n", "model = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n", "@variable(model, a >= 0)\n", "@variable(model, b >= 0)\n", "@constraint(model, a + b >= 1)\n", "peval(a, b) = p(x=>a, y=>b)\n", "register(model, :peval, 2, peval, autodiff=true)\n", "@NLobjective(model, Min, peval(a, b))\n", "optimize!(model)\n", "@show termination_status(model)\n", "@show value(a)\n", "@show value(b)\n", "@show objective_value(model)\n", "\n", "# Sum-of-Squares approach" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We will now see how to find the optimal solution using Sum of Squares Programming.\n", "We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "MathOptInterface.OptimizerWithAttributes(CSDP.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])" }, "metadata": {}, "execution_count": 4 } ], "cell_type": "code", "source": [ "import CSDP\n", "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "A Sum-of-Squares certificate that $p \\ge \\alpha$ over the domain `S`, ensures that $\\alpha$ is a lower bound to the polynomial optimization problem.\n", "The following program searches for the largest upper bound and finds zero." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "termination_status(model) = MathOptInterface.OPTIMAL\n", "objective_value(model) = -1.7799189899747603e-10\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "-1.7799189899747603e-10" }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "model = SOSModel(solver)\n", "@variable(model, α)\n", "@objective(model, Max, α)\n", "@constraint(model, c3, p >= α, domain = S)\n", "optimize!(model)\n", "@show termination_status(model)\n", "@show objective_value(model)" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "Using the solution $(1/2, 1/2)$ found by Ipopt of objective value $1/4$\n", "and this certificate of lower bound $0$ we know that the optimal objective value is in the interval $[0, 1/4]$\n", "but we still do not know what it is (if we consider that we did not try the solutions $(1, 0)$ and $(0, 1)$ as done in the introduction).\n", "If the dual of the constraint `c3` was atomic, its atoms would have given optimal solutions of objective value $0$ but that is not the case." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ν3 = moment_matrix(c3)\n", "extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "Fortunately, there is a hierarchy of programs with increasingly better bounds that can be solved until we get one with atom dual variables.\n", "This comes from the way the Sum-of-Squares constraint with domain `S` is formulated.\n", "The polynomial $p - \\alpha$ is guaranteed to be nonnegative over the domain `S` if there exists Sum-of-Squares polynomials $s_0$, $s_1$, $s_2$, $s_3$ such that\n", "$$\n", "p - \\alpha = s_0 + s_1 x + s_2 y + s_3 (x + y - 1).\n", "$$\n", "Indeed, in the domain `S`, $x$, $y$ and $x + y - 1$ are nonnegative so the right-hand side is a sum of squares hence is nonnegative.\n", "Once the degrees of $s_1$, $s_2$ and $s_3$ have been decided, the degree needed for $s_0$ will be determined but we have a freedom in choosing the degrees of $s_1$, $s_2$ and $s_3$.\n", "By default, they are chosen so that the degrees of $s_1 x$, $s_2 y$ and $s_3 (x + y - 1)$ match those of $p - \\alpha$ but this can be overwritten using the `maxdegree` keyword argument." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### The maxdegree keyword argument" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "The maximum total degree (i.e. maximum sum of the exponents of $x$ and $y$) of the monomials of $p$ is 3 so the constraint in the program above is equivalent to `@constraint(model, p >= α, domain = S, maxdegree = 3)`.\n", "That is, since $x$, $y$ and $x + y - 1$ have total degree 1, the sum of squares polynomials $s_1$, $s_2$ and $s_3$ have been chosen with maximum total degree $2$.\n", "Since these polynomials are sums of squares, their degree must be even so the next maximum total degree to try is 4.\n", "For this reason, the keywords `maxdegree = 4` and `maxdegree = 5` have the same effect in this example.\n", "In general, if the polynomials in the domain are not all odd or all even, each value of `maxdegree` has a different effect in the choice of the maximum total degree of some $s_i$." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "termination_status(model) = MathOptInterface.OPTIMAL\n", "objective_value(model) = -3.4520638569901507e-9\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "-3.4520638569901507e-9" }, "metadata": {}, "execution_count": 7 } ], "cell_type": "code", "source": [ "model = SOSModel(solver)\n", "@variable(model, α)\n", "@objective(model, Max, α)\n", "@constraint(model, c5, p >= α, domain = S, maxdegree = 5)\n", "optimize!(model)\n", "@show termination_status(model)\n", "@show objective_value(model)" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "As shown below, for `maxdegree = 5`, the dual variable is atomic as it is the moments of the measure\n", "$$0.5 \\delta(x-1, y) + 0.5 \\delta(x, y-1)$$\n", "where $\\delta(x, y)$ is the dirac measure centered at $(0, 0)$.\n", "Therefore the program provides both a certificate that $0$ is a lower bound and a certificate that it is also an upper bound since it is attained at the global minimizers $(1, 0)$ and $(0, 1)$." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Atomic measure on the variables x, y with 2 atoms:\n at [-0.00027755883462882424, 1.0002775595143842] with weight 0.4999650862222643\n at [0.9997947978356583, 0.00020520284409647233] with weight 0.5005925828244099" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "ν5 = moment_matrix(c5)\n", "extractatoms(ν5, 1e-3)" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "## A deeper look into atom extraction" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "The `extractatoms` function requires a `ranktol` argument that we have set to `1e-3` in the preceding section.\n", "This argument is used to transform the dual variable into a system of polynomials equations whose solutions give the atoms.\n", "This transformation uses the SVD decomposition of the moment matrix and discards the equations corresponding to a singular value lower than `ranktol`.\n", "When this system of equation has an infinite number of solutions, `extractatoms` concludes that the measure is not atomic.\n", "For instance, with `maxdegree = 3`, we obtain the system\n", "$$x + y = 1$$\n", "which contains a whole line of solution.\n", "This explains `extractatoms` returned `nothing`." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Algebraic Set defined by 1 equalitty\n -x - 1.0000000000000002*y + 1.00000000004984 = 0\n" }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "ν3 = moment_matrix(c3)\n", "SumOfSquares.MultivariateMoments.computesupport!(ν3, 1e-3)" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "With `maxdegree = 5`, we obtain the system\n", "$$\n", "\\begin{aligned}\n", " x + y & = 1\\\\\n", " y^2 & = y\\\\\n", " xy & = 0\\\\\n", " x^2 + y & = 1\n", "\\end{aligned}\n", "$$" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Algebraic Set defined by 3 equalities\n -x - 0.9999999999999994*y + 1.0000000006797547 = 0\n -x*y + 0.44029103700223565*y^2 - 0.44098635462040214*y + 0.00029563385033835106 = 0\n -x^2 + 1.0000005296698848*y^2 - 2.003158975989828*y + 1.0015792234970937 = 0\n" }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "ν5 = moment_matrix(c5)\n", "SumOfSquares.MultivariateMoments.computesupport!(ν5, 1e-3)" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "This system can be reduced to the equivalent system\n", "$$\n", "\\begin{aligned}\n", " x + y & = 1\\\\\n", " y^2 & = y\n", "\\end{aligned}\n", "$$\n", "which has the solutions $(0, 1)$ and $(1, 0)$." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Algebraic Set defined by 2 equalities\n x + 0.9999999999999994*y - 1.0000000006797547 = 0\n y^2 - 1.0004827623584807*y + 0.00020525980009822995 = 0\n" }, "metadata": {}, "execution_count": 11 } ], "cell_type": "code", "source": [ "SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))\n", "ν5.support" ], "metadata": {}, "execution_count": 11 }, { "cell_type": "markdown", "source": [ "The function `extractatoms` then reuses the matrix of moments to find the weights $1/2$, $1/2$ corresponding to the diracs centered respectively at $(0, 1)$ and $(1, 0)$.\n", "This details how the function obtained the result\n", "$$0.5 \\delta(x-1, y) + 0.5 \\delta(x, y-1)$$\n", "given in the previous section." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## HomotopyContinuation" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "As discussed in the previous section, the atom extraction relies on the solution\n", "of a system of algebraic equations. The `extractatoms` function takes an optional\n", "`solver` argument that is used to solve this system of equation.\n", "If no solver is provided, the default solver of SemialgebraicSets.jl is used which\n", "currently computes the Gröbner basis, then the multiplication matrices and\n", "then the Schur decomposition of a random combination of these matrices.\n", "As the system of equations is obtained from a numerical solution and is represented\n", "using floating point coefficients, homotopy continuation is recommended as it is\n", "more numerically robust than Gröbner basis computation.\n", "The following uses homotopy continuation to solve the system of equations." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\rTracking 4 paths... 50%|███████████████▌ | ETA: 0:00:15\u001b[K\r\n", " # paths tracked: 2\u001b[K\r\n", " # non-singular solutions (real): 0 (0)\u001b[K\r\n", " # singular endpoints (real): 0 (0)\u001b[K\r\n", " # total solutions (real): 0 (0)\u001b[K\r\u001b[A\r\u001b[A\r\u001b[A\r\u001b[A\n", "\n", "\n", "\n", "\r\u001b[K\u001b[A\r\u001b[K\u001b[A\r\u001b[K\u001b[A\r\u001b[K\u001b[A\rTracking 4 paths... 100%|███████████████████████████████| Time: 0:00:16\u001b[K\r\n", " # paths tracked: 4\u001b[K\r\n", " # non-singular solutions (real): 0 (0)\u001b[K\r\n", " # singular endpoints (real): 0 (0)\u001b[K\r\n", " # total solutions (real): 0 (0)\u001b[K\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Atomic measure on the variables x, y with 2 atoms:\n at [1.0011790325587666, -0.00033651778432433973] with weight 0.4993038885352588\n at [-0.00017537156118253115, 1.001163506389064] with weight 0.4993385810783572" }, "metadata": {}, "execution_count": 12 } ], "cell_type": "code", "source": [ "using HomotopyContinuation\n", "solver = SemialgebraicSetsHCSolver(; excess_residual_tol = 2e-2, real_tol = 2e-2, compile = false)\n", "extractatoms(ν5, 1e-3, solver)" ], "metadata": {}, "execution_count": 12 }, { "cell_type": "markdown", "source": [ "As the system has 3 equations for 2 variables and the coefficients of the equations\n", "are to be treated with tolerance since they originate from the solution of an SDP,\n", "we need to set `excess_residual_tol` and `real_tol` to a high tolerance otherwise,\n", "HomotopyContinuation would consider that there is no solution.\n", "Indeed, as the system is overdetermined (it has more equations than variables)\n", "HomotopyContinuation expects to have excess solution hence it filters out\n", "excess solution among the solution found. It determines which solution are in excess\n", "by comparing the infinity norm of the residuals of the equations at the solution with `excess_residual_tol`.\n", "It also filters out solution for which the absolute value of the imaginary part of one of the entry\n", "is larger than `real_tol` and strips out the imaginary part.\n", "The raw solutions obtained by HomotopyContinuation can be obtained as follows:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\rTracking 4 paths... 50%|███████████████▌ | ETA: 0:00:05\u001b[K\r\n", " # paths tracked: 2\u001b[K\r\n", " # non-singular solutions (real): 0 (0)\u001b[K\r\n", " # singular endpoints (real): 0 (0)\u001b[K\r\n", " # total solutions (real): 0 (0)\u001b[K\r\u001b[A\r\u001b[A\r\u001b[A\r\u001b[A\n", "\n", "\n", "\n", "\r\u001b[K\u001b[A\r\u001b[K\u001b[A\r\u001b[K\u001b[A\r\u001b[K\u001b[A\rTracking 4 paths... 100%|███████████████████████████████| Time: 0:00:04\u001b[K\r\n", " # paths tracked: 4\u001b[K\r\n", " # non-singular solutions (real): 0 (0)\u001b[K\r\n", " # singular endpoints (real): 0 (0)\u001b[K\r\n", " # total solutions (real): 0 (0)\u001b[K\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "4-element Vector{HomotopyContinuation.PathResult}:\n PathResult:\n • return_code → :excess_solution\n • solution → ComplexF64[1.0008057734171418 + 0.0002048385483798488im, 0.00028102311437224316 - 0.00042280351915176126im]\n • accuracy → 0.0011084\n • residual → 0.00043123\n • condition_jacobian → 27.8\n • steps → 45 / 0\n • extended_precision → false\n • path_number → 1\n\n PathResult:\n • return_code → :excess_solution\n • solution → ComplexF64[-0.40711828679793954 - 2.0709817730382825im, -0.010328770142141848 - 1.7903182953611223im]\n • accuracy → 4.1132\n • residual → 1.4464\n • condition_jacobian → 10.315\n • steps → 68 / 0\n • extended_precision → false\n • path_number → 2\n\n PathResult:\n • return_code → :excess_solution\n • solution → ComplexF64[0.2672505900648712 + 0.15449396539443505im, 0.6648120003182358 + 0.4339123228356712im]\n • accuracy → 0.59232\n • residual → 0.39383\n • condition_jacobian → 4.2396\n • steps → 44 / 0\n • extended_precision → false\n • path_number → 3\n\n PathResult:\n • return_code → :excess_solution\n • solution → ComplexF64[0.0010478559325373322 - 0.00032951295580495765im, 1.0011077722299275 + 0.0013306772319096646im]\n • accuracy → 0.0023768\n • residual → 0.0015729\n • condition_jacobian → 9.8373\n • steps → 44 / 0\n • extended_precision → false\n • path_number → 4\n" }, "metadata": {}, "execution_count": 13 } ], "cell_type": "code", "source": [ "F = HomotopyContinuation.System(ν5.support)\n", "res = HomotopyContinuation.solve(F, solver.options...)\n", "path_results(res)" ], "metadata": {}, "execution_count": 13 }, { "cell_type": "markdown", "source": [ "The printed `residual` above shows why `2e-2` allows to filter how the 2 actual\n", "solutions from the 2 excess solutions." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "---\n", "\n", "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" ], "metadata": {} } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.2" }, "kernelspec": { "name": "julia-1.7", "display_name": "Julia 1.7.2", "language": "julia" } }, "nbformat": 4 }