{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Polynomial Optimization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"Consider the polynomial optimization problem 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$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, 1//4, 0)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DynamicPolynomials\n",
"@polyvar x y\n",
"p = x^3 - x^2 + 2x*y -y^2 + y^3\n",
"using SemialgebraicSets\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)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but [Ipopt](https://github.com/JuliaOpt/Ipopt.jl/) only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\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 http://projects.coin-or.org/Ipopt\n",
"******************************************************************************\n",
"\n",
"status = :Optimal\n",
"getvalue(a) = 0.49999999667053496\n",
"getvalue(b) = 0.4999999966706694\n",
"getobjectivevalue(m) = 0.2499999950059033\n"
]
}
],
"source": [
"using JuMP\n",
"using Ipopt\n",
"m = Model(solver=IpoptSolver(print_level=0))\n",
"@variable m a >= 0\n",
"@variable m b >= 0\n",
"@constraint m a + b >= 1\n",
"@NLobjective(m, Min, a^3 - a^2 + 2a*b - b^2 + b^3)\n",
"status = solve(m)\n",
"@show status\n",
"@show getvalue(a)\n",
"@show getvalue(b)\n",
"@show getobjectivevalue(m);get("
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that an equivalent problem can be formulated reusing the polynomial `p` as follows:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"status = :Optimal\n",
"getvalue(a) = 0.499999995002878\n",
"getvalue(b) = 0.4999999950028773\n",
"getobjectivevalue(m) = 0.24999999250431656\n"
]
}
],
"source": [
"using JuMP\n",
"using Ipopt\n",
"m = Model(solver=IpoptSolver(print_level=0))\n",
"@variable m a >= 0\n",
"@variable m b >= 0\n",
"@constraint m a + b >= 1\n",
"peval(a, b) = p(x=>a, y=>b)\n",
"JuMP.register(m, :peval, 2, peval, autodiff=true)\n",
"@NLobjective(m, Min, peval(a, b))\n",
"status = solve(m)\n",
"@show status\n",
"@show getvalue(a)\n",
"@show getvalue(b)\n",
"@show getobjectivevalue(m);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sum-of-Squares approach\n",
"\n",
"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](http://www.juliaopt.org/) for a list of the available choices."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"using CSDP\n",
"solver = CSDPSolver(printlevel=0);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"using Mosek\n",
"solver = MosekSolver(LOG=0);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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 find zero."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"status = :Optimal\n",
"getobjectivevalue(m) = -2.0093032793155885e-10\n"
]
}
],
"source": [
"using JuMP\n",
"using SumOfSquares\n",
"m = SOSModel(solver = solver)\n",
"@variable m α\n",
"@objective m Max α\n",
"c3 = @constraint m p >= α domain = S\n",
"\n",
"status = solve(m)\n",
"@show status\n",
"@show getobjectivevalue(m);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"using MultivariateMoments\n",
"μ3 = -1 * getdual(c3) # /!\\ Need to multiply by -1 because the constraint is \">=\" and not \"<=\"\n",
"X3 = certificate_monomials(c3)\n",
"ν3 = matmeasure(μ3, X3)\n",
"extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fortunately, there is a hierarchy of programs with increasingly better programs 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",
"$$ p - \\alpha = s_0 + s_1 x + s_2 y + s_3 (x + y - 1). $$\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 freesom 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.\n",
"\n",
"### The maxdegree keyword argument\n",
"\n",
"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 m 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 different effect in the choice of the maximum total degree of $s_i$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"status = :Optimal\n",
"getobjectivevalue(m) = -8.670577589242612e-10\n"
]
}
],
"source": [
"using JuMP\n",
"using SumOfSquares\n",
"m = SOSModel(solver = solver)\n",
"@variable m α\n",
"@objective m Max α\n",
"c5 = @constraint m p >= α domain = S maxdegree = 5\n",
"\n",
"status = solve(m)\n",
"@show status\n",
"@show getobjectivevalue(m);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This time, 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)$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Atomic measure on the variables x, y with 2 atoms:\n",
" at [-0.00109071, 1.00109] with weight 0.4990820208096022\n",
" at [0.99992, 8.03849e-5] with weight 0.5006917404344113"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using MultivariateMoments\n",
"μ5 = -1 * getdual(c5) # /!\\ Need to multiply by -1 because the constraint is \">=\" and not \"<=\"\n",
"X5 = certificate_monomials(c5)\n",
"ν5 = matmeasure(μ5, X5)\n",
"extractatoms(ν5, 1e-3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A deeper look into atom extraction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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 matrix of moments and discard 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`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Algebraic Set defined by 1 equality\n",
" -x - 0.9999999999999996*y + 1.0000000000569829 == 0\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ν3 = matmeasure(μ3, X3)\n",
"MultivariateMoments.computesupport!(ν3, 1e-3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With `maxdegree = 5`, we obtain the system\n",
"\\begin{align}\n",
" x + y & = 1\\\\\n",
" y^2 & = y\\\\\n",
" xy & = 0\\\\\n",
" x^2 + y & = 1\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Algebraic Set defined by 4 equalities\n",
" -x - 1.0000000000000002*y + 1.0000000002403162 == 0\n",
" -y^2 + 1.0011710975756178*y - 8.047253705778523e-5 == 0\n",
" -x*y - 3.157029709751049e-12*y + 0.0001853598953271854 == 0\n",
" -x^2 - 1.0011710975546133*y + 1.001090625265065 == 0\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ν5 = matmeasure(μ5, X5)\n",
"MultivariateMoments.computesupport!(ν5, 1e-3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This system can be reduced to the equivalent system\n",
"\\begin{align}\n",
" x + y & = 1\\\\\n",
" y^2 & = y\n",
"\\end{align}\n",
"which has the solutions $(0, 1)$ and $(1, 0)$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Algebraic Set defined by 2 equalities\n",
" x + 1.0000000000000002*y - 1.0000000002403162 == 0\n",
" y^2 - 1.0011710975756178*y + 8.047253705778523e-5 == 0\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))\n",
"ν5.support"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function `extractatoms` then reuse 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 the 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.0",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}