Multivariate polynomials implementations

Contributed by: Benoît Legat

The SumOfSquares package is built on top of the MultivariatePolynomials abstract interface. DynamicPolynomials is an implementation of this abstract interface so it can be used with SumOfSquares. Moreover, any other implementation can be used as well. To illustrate, we solve Examples 3.38 of [BPT12] with TypedPolynomials, another implementation of MultivariatePolynomials.

[BPT12] Blekherman, G. & Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

import TypedPolynomials
TypedPolynomials.@polyvar x y
using SumOfSquares
import CSDP
model = SOSModel(CSDP.Optimizer)
con_ref = @constraint(model, 2x^4 + 5y^4 - x^2*y^2 >= -2(x^3*y + x + 1))
optimize!(model)
solution_summary(model)

* Solver : CSDP

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Problem solved to optimality."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 0.00000e+00
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 6.60331e-02

We see that the problem is feasible. The Sum-of-Squares decomposition can be obtained as follows:

sos_decomposition(con_ref)

(-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

Why is there several implementations ? Depending in the use-case, one implementation may be more appropriate than another one. TypedPolynomials is faster than DynamicPolynomials but it requires new compilation whenever the list of variables changes. This means that TypedPolynomials is not appropriate when the number of variables is dynamic or too large. However, for a small number of variables, it can be faster. When solving Sum-of-Squares programs, the time is mostly taken by the Semidefinite programming solver. The time taken by SumOfSquares/JuMP/MathOptInterface are usually negligible or it time is taken by manipulation of JuMP or MathOptInterface functions therefore using TypedPolynomials over DynamicPolynomials may not make much difference in most cases.

One case for which using TypedPolynomials might be adequate is when using domain defined by equalities (possibly also with inequalities). Indeed, in that case, SumOfSquares computes the corresponding Gröbner basis which may take a non-negligible amount of time for large systems of equalities.

To illustrate this, consider the computation of Gröbner basis for the following system from [CLO05, p. 17]. The time taken by TypedPolynomials is below:

[CLO05] Cox, A. David & Little, John & O'Shea, Donal Using Algebraic Geometry. Graduate Texts in Mathematics, 2005. https://doi.org/10.1007/b138611

using BenchmarkTools
@btime let
    TypedPolynomials.@polyvar x y
    S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
    SemialgebraicSets.compute_gröbner_basis!(S.I)
end

true

The time taken by DynamicPolynomials is as follows:

import DynamicPolynomials
@btime let
    DynamicPolynomials.@polyvar x y
    S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
    SemialgebraicSets.compute_gröbner_basis!(S.I)
end

true

We see that TypedPolynomials is faster. The time is still negligible for this small system but for larger systems, choosing TypedPolynomials may be helpful. We can use this system in a Sum-of-Squares constraint as follows:

TypedPolynomials.@polyvar x y
S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
poly = -6x - 4y^3 + 2x*y^2 + 6x^3 - 3y^4 + 13x^2 * y^2
model = Model(CSDP.Optimizer)
con_ref = @constraint(model, poly in SOSCone(), domain = S)
optimize!(model)
solution_summary(model)

* Solver : CSDP

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Problem solved to optimality."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 0.00000e+00
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 1.85370e-01

We obtain the following decomposition:

dec = sos_decomposition(con_ref, 1e-6)

(-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

We can verify that it is correct as follows:

rem(dec - poly, S.I)

\[ 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} \]

Note that the difference between dec and poly is larger than between the full gram matrix because dec is obtained by dropping the lowest eigenvalues with the threshold 1e-6; see sos_decomposition.

rem(gram_matrix(con_ref) - poly, S.I)

\[ 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} \]

This page was generated using Literate.jl.