{ "cells": [ { "cell_type": "markdown", "source": [ "# Multivariate polynomials implementations" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "**Contributed by**: Benoît Legat" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "The SumOfSquares package is built on top of the [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl)\n", "abstract interface. [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl/)\n", "is an implementation of this abstract interface so it can be used with\n", "SumOfSquares. Moreover, any other implementation can be used as well. To\n", "illustrate, we solve Examples 3.38 of [BPT12] with\n", "[TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl),\n", "another implementation of [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl).\n", "\n", "[BPT12] Blekherman, G. & Parrilo, P. A. & Thomas, R. R.\n", "*Semidefinite Optimization and Convex Algebraic Geometry*.\n", "Society for Industrial and Applied Mathematics, **2012**." ], "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.39e-01 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 9.2771541e+01 \n", "Iter: 2 Ap: 1.00e+00 Pobj: 0.0000000e+00 Ad: 1.00e+00 Dobj: 5.3405053e+01 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "* Solver : CSDP\n\n* Status\n Result count : 1\n Termination status : OPTIMAL\n Message from the solver:\n \"Problem solved to optimality.\"\n\n* Candidate solution (result #1)\n Primal status : FEASIBLE_POINT\n Dual status : FEASIBLE_POINT\n Objective value : 0.00000e+00\n Dual objective value : 0.00000e+00\n\n* Work counters\n Solve time (sec) : 1.44833e-01\n" }, "metadata": {}, "execution_count": 1 } ], "cell_type": "code", "source": [ "import TypedPolynomials\n", "TypedPolynomials.@polyvar x y\n", "using SumOfSquares\n", "import CSDP\n", "model = SOSModel(CSDP.Optimizer)\n", "con_ref = @constraint(model, 2x^4 + 5y^4 - x^2*y^2 >= -2(x^3*y + x + 1))\n", "optimize!(model)\n", "solution_summary(model)" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "We see that the problem is feasible. The Sum-of-Squares decomposition can be\n", "obtained as follows:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "(-0.679754396922155 - 0.09754507709792247*y - 0.41759505637134503*x + 2.1931518901396183*y^2 - 0.12696062446783402*x*y - 0.8025867547845225*x^2)^2 + (0.9043573827480677 - 0.6455003192597272*y + 0.2742159536999836*x - 0.10770174782358193*y^2 - 1.229313806880072*x*y - 0.9300160939847002*x^2)^2 + (-0.5739168240154697 - 1.385350502215481*y - 0.6219655378569914*x - 0.2040924100740059*y^2 - 0.15182215880489808*x*y + 0.4443832632766483*x^2)^2 + (0.5703021360913619 - 0.5184093413621604*y + 0.3436728787624475*x + 0.2555950605470667*y^2 + 0.7618441294650051*x*y - 0.0209056737676136*x^2)^2 + (0.04140083721776801 + 0.049110402726379475*y - 0.49990658891857376*x - 0.2479214593908115*y^2 + 0.297244568991742*x*y - 0.5054184235227421*x^2)^2 + (-0.252453178814133 - 0.06336647012459035*y + 0.25394257222405114*x - 0.1001874415304813*y^2 + 0.05960179598589515*x*y - 0.1938124125581194*x^2)^2" }, "metadata": {}, "execution_count": 2 } ], "cell_type": "code", "source": [ "sos_decomposition(con_ref)" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Why is there several implementations ?\n", "Depending in the use-case, one implementation may be more appropriate than\n", "another one. [TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl)\n", "is faster than [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl/)\n", "but it requires new compilation whenever the list of variables changes.\n", "This means that [TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl)\n", "is not appropriate when the number of variables is dynamic or too large.\n", "However, for a small number of variables, it can be faster.\n", "When solving Sum-of-Squares programs, the time is mostly taken by the Semidefinite programming solver.\n", "The time taken by SumOfSquares/JuMP/MathOptInterface are usually negligible\n", "or it time is taken by manipulation of JuMP or MathOptInterface functions\n", "therefore using TypedPolynomials over DynamicPolynomials may not make much difference in most cases." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "One case for which using TypedPolynomials might be adequate is when\n", "using domain defined by equalities (possibly also with inequalities).\n", "Indeed, in that case, SumOfSquares computes the corresponding Gröbner basis which\n", "may take a non-negligible amount of time for large systems of equalities." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "To illustrate this, consider the computation of Gröbner basis for the\n", "following system from [CLO05, p. 17].\n", "The time taken by TypedPolynomials is below:\n", "\n", "[CLO05] Cox, A. David & Little, John & O'Shea, Donal\n", "*Using Algebraic Geometry*.\n", "Graduate Texts in Mathematics, **2005**.\n", "https://doi.org/10.1007/b138611" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 870.232 μs (6614 allocations: 1.34 MiB)\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "true" }, "metadata": {}, "execution_count": 3 } ], "cell_type": "code", "source": [ "using BenchmarkTools\n", "@btime let\n", " TypedPolynomials.@polyvar x y\n", " S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y\n", " SemialgebraicSets.compute_gröbner_basis!(S.I)\n", "end" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "The time taken by DynamicPolynomials is as follows:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1.358 ms (9990 allocations: 1.26 MiB)\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "true" }, "metadata": {}, "execution_count": 4 } ], "cell_type": "code", "source": [ "import DynamicPolynomials\n", "@btime let\n", " DynamicPolynomials.@polyvar x y\n", " S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y\n", " SemialgebraicSets.compute_gröbner_basis!(S.I)\n", "end" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "We see that TypedPolynomials is faster.\n", "The time is still negligible for this small system but for larger systems, choosing TypedPolynomials may be helpful.\n", "We can use this system in a Sum-of-Squares constraint as follows:" ], "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: 2.23e-16 \n", "Relative dual infeasibility: 5.00e-11 \n", "Real Relative Gap: 0.00e+00 \n", "XZ Relative Gap: 3.04e-10 \n", "DIMACS error measures: 2.78e-16 0.00e+00 5.00e-11 0.00e+00 0.00e+00 3.04e-10\n", "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.64e-01 Pobj: 0.0000000e+00 Ad: 8.07e-01 Dobj: 2.6104524e+00 \n", "Iter: 2 Ap: 7.71e-01 Pobj: 0.0000000e+00 Ad: 7.40e-01 Dobj: -5.9092654e-02 \n", "Iter: 3 Ap: 8.08e-01 Pobj: 0.0000000e+00 Ad: 7.65e-01 Dobj: 5.5549062e-01 \n", "Iter: 4 Ap: 7.56e-01 Pobj: 0.0000000e+00 Ad: 7.79e-01 Dobj: 4.8506833e-01 \n", "Iter: 5 Ap: 6.95e-01 Pobj: 0.0000000e+00 Ad: 6.77e-01 Dobj: 4.3029685e-01 \n", "Iter: 6 Ap: 6.35e-01 Pobj: 0.0000000e+00 Ad: 6.83e-01 Dobj: 2.2137644e-01 \n", "Iter: 7 Ap: 6.21e-01 Pobj: 0.0000000e+00 Ad: 6.43e-01 Dobj: 1.2986372e-01 \n", "Iter: 8 Ap: 7.42e-01 Pobj: 0.0000000e+00 Ad: 6.82e-01 Dobj: 6.1463607e-02 \n", "Iter: 9 Ap: 7.31e-01 Pobj: 0.0000000e+00 Ad: 7.05e-01 Dobj: 2.6276582e-02 \n", "Iter: 10 Ap: 6.92e-01 Pobj: 0.0000000e+00 Ad: 7.09e-01 Dobj: 1.1328533e-02 \n", "Iter: 11 Ap: 6.62e-01 Pobj: 0.0000000e+00 Ad: 6.68e-01 Dobj: 5.1233630e-03 \n", "Iter: 12 Ap: 7.09e-01 Pobj: 0.0000000e+00 Ad: 6.82e-01 Dobj: 2.3589113e-03 \n", "Iter: 13 Ap: 7.23e-01 Pobj: 0.0000000e+00 Ad: 6.78e-01 Dobj: 1.0697197e-03 \n", "Iter: 14 Ap: 7.14e-01 Pobj: 0.0000000e+00 Ad: 7.19e-01 Dobj: 4.1892817e-04 \n", "Iter: 15 Ap: 6.78e-01 Pobj: 0.0000000e+00 Ad: 6.29e-01 Dobj: 2.1663343e-04 \n", "Iter: 16 Ap: 7.27e-01 Pobj: 0.0000000e+00 Ad: 6.95e-01 Dobj: 8.2602885e-05 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "* Solver : CSDP\n\n* Status\n Result count : 1\n Termination status : OPTIMAL\n Message from the solver:\n \"Problem solved to optimality.\"\n\n* Candidate solution (result #1)\n Primal status : FEASIBLE_POINT\n Dual status : FEASIBLE_POINT\n Objective value : 0.00000e+00\n Dual objective value : 0.00000e+00\n\n* Work counters\n Solve time (sec) : 4.42834e-01\n" }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "TypedPolynomials.@polyvar x y\n", "S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y\n", "poly = -6x - 4y^3 + 2x*y^2 + 6x^3 - 3y^4 + 13x^2 * y^2\n", "model = Model(CSDP.Optimizer)\n", "con_ref = @constraint(model, poly in SOSCone(), domain = S)\n", "optimize!(model)\n", "solution_summary(model)" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "We obtain the following decomposition:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "(-6.841213011209407e-5 + 0.0006527587940908881*y - 0.0001414825075640032*x + 7.044903864853802*y^2 - 28.22963294804742*x*y + 25.36362402970179*x^2)^2 + (-1.4193698655877476e-5 + 0.00042745768528049843*y + 0.0008136661561626277*x + 22.4793208244021*y^2 - 13.864881630766194*x*y - 21.675339966436*x^2)^2 + (-3.8438856467205195e-5 - 0.000728530340103983*y + 0.0005524266995899736*x + 0.12018765789331433*y^2 + 0.09016503612816339*x*y + 0.06697055310943034*x^2)^2 + (-6.125566301704854e-7 - 0.03996374501021478*y + 0.053403305256105094*x - 0.0002606954900837711*y^2 - 0.00019563699075496687*x*y - 0.0001440071613000224*x^2)^2" }, "metadata": {}, "execution_count": 6 } ], "cell_type": "code", "source": [ "dec = sos_decomposition(con_ref, 1e-6)" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "We can verify that it is correct as follows:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "6.35960154013951484141801108397780428749257453091559000313282012939453125e-09 - 2.168942626362357496871710423330011975642296653323006840088428595127160415659131e-07y - 1.62061044993978297358038477026499234713040865384615384615384615384614703022391e-07x - 1.304211826624705750088395461716572754085063934326171875e-05y² - 1.95310533733401847644728377417777664959430694580078125e-05xy - 7.300716056938584552771231983570032753050327301025390625e-06x² - 4.627175886895429357537068426609039306640625e-08y³ - 4.97990413350635208189487457275390625e-08xy² - 4.160468380499289802083862088721843974781222641468048095703125e-08x²y + 6.894101016559244061891849224384014423076923076923076923076923076917624179115685e-09y⁴", "text/latex": "$$ 6.35960154013951484141801108397780428749257453091559000313282012939453125 \\cdot 10^{-09} - 2.168942626362357496871710423330011975642296653323006840088428595127160415659131 \\cdot 10^{-07}y - 1.62061044993978297358038477026499234713040865384615384615384615384614703022391 \\cdot 10^{-07}x - 1.304211826624705750088395461716572754085063934326171875 \\cdot 10^{-05}y^{2} - 1.95310533733401847644728377417777664959430694580078125 \\cdot 10^{-05}xy - 7.300716056938584552771231983570032753050327301025390625 \\cdot 10^{-06}x^{2} - 4.627175886895429357537068426609039306640625 \\cdot 10^{-08}y^{3} - 4.97990413350635208189487457275390625 \\cdot 10^{-08}xy^{2} - 4.160468380499289802083862088721843974781222641468048095703125 \\cdot 10^{-08}x^{2}y + 6.894101016559244061891849224384014423076923076923076923076923076917624179115685 \\cdot 10^{-09}y^{4} $$" }, "metadata": {}, "execution_count": 7 } ], "cell_type": "code", "source": [ "rem(dec - poly, S.I)" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "Note that the difference between `dec` and `poly` is larger\n", "than between the full gram matrix because `dec` is obtained by dropping\n", "the lowest eigenvalues with the threshold `1e-6`; see `sos_decomposition`." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "7.413974285140153089766624073912726355928271004813723266124725341796875e-09 - 5.030138093290740859046486526160958983557345793367578433110163762005625482154789e-10y - 1.082308943322319604563885010205782376802884615384615384615384615377812741077394e-10x + 7.4139898807897008925493764763814397156238555908203125e-09y² - 5.897886345973546440291102044284343719482421875e-14xy + 7.4139714055510008705596192157827317714691162109375e-09x² - 7.712349277729420767476161321004231770833333333333333333333333333563631161469185e-13y³ - 2.565355335567195046072204907735188802083333333333333333333333333218184419265407e-13xy² + 1.4523104940877828994416631758213043212890625e-13x²y + 6.952334912896801072817582350510817307692307692307692307692307692302250192860886e-09y⁴", "text/latex": "$$ 7.413974285140153089766624073912726355928271004813723266124725341796875 \\cdot 10^{-09} - 5.030138093290740859046486526160958983557345793367578433110163762005625482154789 \\cdot 10^{-10}y - 1.082308943322319604563885010205782376802884615384615384615384615377812741077394 \\cdot 10^{-10}x + 7.4139898807897008925493764763814397156238555908203125 \\cdot 10^{-09}y^{2} - 5.897886345973546440291102044284343719482421875 \\cdot 10^{-14}xy + 7.4139714055510008705596192157827317714691162109375 \\cdot 10^{-09}x^{2} - 7.712349277729420767476161321004231770833333333333333333333333333563631161469185 \\cdot 10^{-13}y^{3} - 2.565355335567195046072204907735188802083333333333333333333333333218184419265407 \\cdot 10^{-13}xy^{2} + 1.4523104940877828994416631758213043212890625 \\cdot 10^{-13}x^{2}y + 6.952334912896801072817582350510817307692307692307692307692307692302250192860886 \\cdot 10^{-09}y^{4} $$" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "rem(gram_matrix(con_ref) - poly, S.I)" ], "metadata": {}, "execution_count": 8 }, { "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.9.3" }, "kernelspec": { "name": "julia-1.9", "display_name": "Julia 1.9.3", "language": "julia" } }, "nbformat": 4 }