Sums of squares in Macaulay2

TL;DR

The paper introduces a Macaulay2 package that efficiently decomposes polynomials into sums of squares using advanced methods, supporting various optimization techniques and external solvers.

Contribution

It presents a new Macaulay2 package with capabilities for sums-of-squares decomposition, integrating modern algorithms and solver support.

Findings

Supports external semidefinite programming solvers

Includes a data type for sums-of-squares polynomials

Enables optimization over algebraic varieties

Abstract

The Macaulay2 package SumsOfSquares decomposes polynomials as sums of squares. It is based on methods to rationalize sum-of-squares decompositions due to Parrilo and Peyrl. The package features a data type for sums-of-squares polynomials, support for external semidefinite programming solvers, and optimization over varieties.

Equations

$$f\text{ is SOS }⟺\exists Q\in {\mathbb{S}}_{+}^N\text{ such that }f=v^{T}Qv,$$

$$f=\frac{p_1\cdots p_s}{q_1\cdots q_t}\text{ pi,qi SOS}.$$

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x)\text{ such that }{h}_1(x)=\cdots =h_m(x)=0,$$

Full text

1. Introduction

Let $\mathbb{K}\!=\!\mathbb{Q}$ or $\mathbb{K}\!=\!\mathbb{R}$ be the rational or real numbers and $R=\mathbb{K}[x_{1},\dots,x_{n}]$ be the polynomial ring. An element $f\!\in\!R$ is nonnegative if $f(x)\!\geq\!0$ for all $x\!\in\!\mathbb{R}^{n}$, and $f$ is a sum of squares (SOS) if there are polynomials $f_{1},\dots,f_{m}\!\in\!R$ and positive scalars $\lambda_{1},\dots,\lambda_{m}\!\in\!\mathbb{K}$ such that $f=\sum_{i}\lambda_{i}f_{i}^{2}$. The scalars are not necessary when the field is $\mathbb{K}\!=\!\mathbb{R}$. Clearly, a sum of squares is nonnegative, but not every nonnegative polynomial is a sum of squares. Hilbert showed that the nonnegative polynomials of degree $d$ in $n$ variables are sums of squares if and only if: $n=1$; or $d=2$; or $n=2$ and $d=4$. For an introduction to the area we recommend [9, 2].

The SumsOfSquares package contains methods to compute sums-of-squares in Macaulay2 [6]. A particular focus is on trying to find rational sums-of-squares decompositions of polynomials with rational coefficients (whenever they exist).

Consider the basic problem of deciding whether a polynomial is a sum of squares. Let $f\!\in\!R$ of degree $2d$, and $v\!\in\!R^{N}$ a vector whose entries are the $N\!=\!\binom{n+d}{d}$ monomials of degree $\leq\!d$. The following fundamental result holds:

[TABLE]

where $\mathbb{S}_{+}^{N}$ is the cone of $N{\times}N$ symmetric positive semidefinite matrices; see [2, §3.1]. This reduces the problem to finding a Gram matrix $Q$ as above, which can be done efficiently with semidefinite programming (SDP).

The method solveSOS performs the computation above. We use it here to verify that $f=2x^{4}{+}5y^{4}{-}2x^{2}y^{2}{+}2x^{3}y$ is a sum of squares:

i1 : R = QQ[x,y]; i2 : f = 2x^4+5y^4-2x^2y^2+2x^3y; i3 : sol = solveSOS f; Executing CSDP Status: SDP solved, primal-dual feasible

The “Status” line indicates that a Gram matrix was found, so $f$ is indeed a sum of squares. In the example above the package called an external program to serve as semidefinite programming solver. The default solver is the open source program CSDP [3], which is included in Macaulay2. The output of solveSOS is an object of type SDPResult. It contains, in particular, the Gram matrix $Q$ and the monomial vector $v$.

i4 : (Q,v) = ( sol#GramMatrix, sol#Monomials ) o4 = ( | 2 1 -83/40 |, | x2 | ) | 1 43/20 0 | | xy | | -83/40 0 5 | | y2 |

The result of the semidefinite programming solver is a floating point approximation of the Gram matrix. The SumsOfSquares package attempts to find a close enough rational Gram matrix by rounding its entries [8]. If this rounding procedure fails to find a feasible rational matrix, the method returns the floating point solution. The procedure is guaranteed to work when the floating point Gram matrix lies in the interior of $\mathbb{S}_{+}^{N}$. See Appendix A for more details about rational rounding.

The method sosPoly extracts the sum-of-squares decomposition from the returned SDPResult. This is done via an LDL factorization (a variant of Cholesky factorization) of the Gram matrix. For the function $f$ from above we get three squares:

i5 : s = sosPoly sol 83 2 2 2 43 20 2 2 231773 2 2 o5 = (5)(- ---x + y ) + (--)(--x + x*y) + (------)(x ) 200 20 43 344000

The output above is an object of type SOSPoly. An object of this type stores the coefficients $\lambda_{i}$ and polynomials (or generators) $f_{i}$ such that $f=\sum_{i}\lambda_{i}f_{i}^{2}$. We can extract the coefficients and generators as follows:

i5 : coefficients s o5 = {5, 43/20, 231773/344000} i6 : gens s o6 = {-83/200x^2 + y^2, 20/43x^2 + x*y, x^2}

The method solveSOS can also compute sums-of-squares decompositions in quotient rings. This can be useful to prove nonnegativity of a polynomial on a variety. We take an example from [7]. Consider proving that $f=10{-}x^{2}{-}y$ is nonnegative on the circle defined by $g=x^{2}{+}y^{2}{-}1$. To do this, we check if $f$ is a sum of squares in the quotient ring $\mathbb{Q}[x,y]/\langle g\rangle$. For such a computation, an even degree bound must be given by the user, as otherwise it is not obvious how to choose the monomial vector $v$. In the following example we use $2d=2$ as the degree bound.

i1 : R = QQ[x,y]/ideal(x^2 + y^2 - 1); i2 : f = 10-x^2-y; d = 1; i3 : sosPoly solveSOS (f, 2*d, TraceObj=>true) Executing CSDP Status: SDP solved, primal-dual feasible 1 2 35 2 o3 = (9)(- -- y + 1) + (--)(y) 18 36

In the computation above the option TraceObj=>true was used to reduce the number of squares in the SOS decomposition (see Section 6).

2. Sums of squares in ideals

Let $I\subset\mathbb{K}[x_{1},\dots,x_{n}]$ be an ideal. Given an even bound $2d$, consider the problem of finding a nonzero sum-of-squares polynomial of degree $\leq\!2d$ in the ideal $I$. If one of the generators of $I$ has degree $\leq\!d$, then the problem is trivial. But otherwise the problem can be hard. The method sosInIdeal can be used to solve it. One of the main motivations for this problem is that it reveals information about the real radical of the ideal $I$, i.e., the vanishing ideal of the real zeros of $I$. Indeed, if $f=\sum\lambda_{i}f_{i}^{2}\in I$ then each of the factors $f_{i}$ must lie in the real radical of $I$.

Given generators of the ideal $I=\langle h_{1},\dots,h_{m}\rangle$, we may solve this problem by looking for some polynomial multipliers $l_{i}(x)$ such that $\sum_{i}l_{i}(x)h_{i}(x)$ is a sum of squares. The method sosInIdeal can find these multipliers. The input is a matrix containing the generators, and the degree bound $2d$. We illustrate this for the ideal $I=\langle x^{2}{-}4x{+}2y^{2},2z^{2}{-}y^{2}{+}2\rangle$

i1 : R = QQ[x,y,z]; d = 1; i2 : h = matrix {{x^2-4x+2y^2, 2z^2-y^2+2}}; i3 : (sol,mult) = sosInIdeal (h, 2d); i4 : sosPoly sol 395 1 2 395 2 o4 = (---)(- - x + 1) + (---)(z) 2 2 2 i5 : h * mult == sosPoly sol o5 = true

Another way to approach this problem is to construct the quotient $S=\mathbb{K}[x_{1},\dots,x_{n}]/I$ and then write $0\in S$ as a sum of squares. In this case the input to sosInIdeal is simply the quotient ring $S$.

i6 : S = R/ideal h; i7 : sosPoly sosInIdeal (S, 2*d); 1031833 1 2 1031833 2 o7 = (-------)(- - x + 1) + (-------)(z) 2048 2 2048

In both cases we obtained a multiple of the sum-of-squares polynomial $(\frac{1}{2}x{-}1)^{2}{+}z^{2}$. This computation reveals that $x{-}2,z$ lie in the real radical of $I$. Indeed, we have $\sqrt[\mathbb{R}]{I}=\langle x{-}2,z,y^{2}{-}2\rangle$.

3. SOS decompositions of ternary forms

Hilbert showed that any nonnegative form $f\in\mathbb{K}[x,y,z]$ can be decomposed as a quotient of sums of squares. We can obtain this decomposition by iteratively calling sosInIdeal. Specifically, one can first find a multiplier $q_{1}$ such that $q_{1}f$ is a sum of squares. Since $q_{1}$ is also nonnegative, we can then search for a multiplier $p_{1}$ such that $p_{1}q_{1}$ is a sum of squares, and so on. The main observation is that the necessary degree of $p_{1}$ is lower than that of $q_{1}$ [5]. Hence this procedure terminates, and we can write

[TABLE]

As an illustration, we write the Motzkin polynomial as a quotient of sums of squares. We first use the function library, which contains a small library of interesting nonnegative forms.

i1 : R = QQ[x,y,z] i2 : f = library ("Motzkin", {x,y,z}) 4 2 2 4 2 2 2 6 o2 = x y + x y - 3x y z + z

We now apply the function sosdecTernary, which implements the iterative algorithm from above.

i3 : (Nums,Dens) = sosdecTernary f; Executing CSDP i4 : num = first Nums 2267 2 2 4 2 2003 1013 3 990 3 2 2 o4 = (----)(x y - z ) + (----)(- ----x y - ----xy + xy*z ) + ... 64 64 2003 2003 i5 : den = first Dens 2267 2 1079 2 33 2 o5 = (----)(z) + (----)(x) + (--)(y) 64 64 2

The result consists of two sums of squares, the second being the denominator. We can check the computation as follows.

i6 : f*value(den) == value(num) o6 = true

4. Parametric SOS problems

The SumsOfSquares package can also solve parametric problems. Assume now that $x\mapsto f(x;t)$ is a polynomial function, that depends affinely on some parameters $t$. The command solveSOS can be used to search for values of the parameters such that the polynomial is a sum of squares. In the following example, we change two coefficients of the Robinson polynomial so that it becomes a sum of squares.

i1 : R = QQ[x,y,z][s,t]; i2 : g = library("Robinson", {x,y,z}) + sx^6 + ty^6; i3 : sol = solveSOS g; Executing CSDP Status: SDP solved, primal-dual feasible i4 : sol#Parameters o4 = | 34 | | 34 |

In the code above, the ring construction (first line) indicates that $s,t$ should be treated as parameters. The values obtained were $s=t=34$.

It is also possible find the values of the parameters that optimize a given linear function. This allows us to find lower bounds for a polynomial function $f(x)$, by finding the largest $t$ such that $f(x)-t$ is a sum of squares. Here we apply this method to the dehomogenized Motzkin polynomial.

i1 : R = QQ[x,z][t]; i2 : f = library ("Motzkin", {x,1,z}); i3 : sol = solveSOS (f-t, -t, RoundTol=>12); Executing CSDP Status: SDP solved, primal-dual feasible i4 : sol#Parameters o4 = | -729/4096 |

Alternatively, the method lowerBound can be called with input $f(x)$. The method internally declares a new parameter $t$ and optimizes $f(x)-t$.

i1 : R = QQ[x,z]; i2 : f = library ("Motzkin", {x,1,z}); i3 : (t,sol) = lowerBound (f, RoundTol=>12); Executing CSDP Status: SDP solved, primal-dual feasible i4 : t o4 = - 729/4096

5. Polynomial optimization

In applications one often needs to find lower bounds for polynomials subject to some polynomial constraints. More precisely, consider the problem

[TABLE]

where $f,h_{1},\dots,h_{m}$ are polynomials. The SumsOfSquares package provides two ways to compute a lower bound for such a problem. The most elegant approach is to construct the associated quotient ring, and then call lowerBound. This will look for the largest $t$ such that $f(x)-t$ is a sum of squares (in the quotient ring). A degree bound $2d$ must be given by the user.

i1 : R = QQ[x,y]/ideal(x^2 - x, y^2 - y); i2 : f = x - y; d = 1; i3 : (t,sol) = lowerBound(f,2*d); Executing CSDP Status: SDP solved, primal-dual feasible i4 : t o4 = -1 i5 : f - t == sosPoly sol o5 = true

Calling lowerBound as above is conceptually simple, but requires knowledge of a Gröbner basis, which is computed when constructing the quotient ring. If no Gröbner basis is available there is an alternative way to call lowerBound with just the equations $h_{1},\dots,h_{m}$ as the input. The method will then look for polynomial multipliers $l_{i}(x)$ such that $f(x)-t+\sum_{i}l_{i}(x)h_{i}(x)$ is a sum of squares. This may result in larger semidefinite programs and weaker bounds.

i1 : R = QQ[x,y]; i2 : f = x - y; d = 1; i3 : h = matrix{{x^2 - x, y^2 - y}}; i4 : (t,sol,mult) = lowerBound (f, h, 2d); Executing CSDP Status: SDP solved, primal-dual feasible i5 : t o5 = -1 i6 : f - t + hmult == sosPoly sol o6 = true

Lower bounds for polynomial optimization problems critically depend on the degree bound chosen. While higher degree bounds lead to better bounds, the computational complexity escalates quite rapidly. Nonetheless, low degree SOS lower bounds often perform very well in applications. In some cases, the minimizer might also be recovered from the SDPResult with the method recoverSolution.

i7 : recoverSolution sol o7 = {x => 1.77345e-9, y => 1}

6. Optional arguments

SDP Solver

The optional argument Solver is available for many package methods and a particular semidefinite programming solver can be picked by setting it. These solvers are interfaced via the auxiliary Macaulay2 package SemidefiniteProgramming [4]. The package provides interfaces to the open source solvers CSDP [3] and SDPA [10], and the commercial solver MOSEK [1]. There is also a built-in solver in the Macaulay2 language. In our experience CSDP and MOSEK give the best results. CSDP is provided as part of Macaulay2 and configured as the default.

Rounding tolerance

The method lowerBound has the optional argument RoundTol, which specifies the precision of the rational rounding. Smaller values of RoundTol lead to rational matrices with smaller denominators but farther from the numerical solution. The rational rounding may be skipped by setting it to infinity.

Trace objective

The option TraceObj tells the semidefinite programming solver to minimize the trace of the Gram matrix. This is a known heuristic to reduce the number of squares in the SOS decomposition.

Appendix A Rational rounding

Sums-of-squares problems are solved numerically using an semidefinite programming solver, and afterwards the package attempts to round the floating point solution to rational numbers. We briefly describe the rounding procedure, which was proposed in [8].

Let $f\in\mathbb{Q}[x_{1},\dots,x_{n}]$ and consider the affine space $\mathcal{L}:=\{Q:v^{T}\!Qv\!=\!f\}$, where $v$ is a given monomial vector. A Gram matrix is an element of $\mathcal{L}\cap\mathbb{S}_{+}^{N}$. The semidefinite programming solver returns a numerical matrix $Q_{n}$, an “approximate” Gram matrix, which may not lie exactly on $\mathcal{L}$. The rounding problem consists in finding a nearby Gram matrix $Q_{r}$ with rational entries.

The procedure from [8] consists of two steps. First, the entries of $Q_{n}$ are rounded to a rational matrix $Q_{r}^{\prime}$. Then $Q_{r}$ is obtained as the orthogonal projection of $Q_{r}^{\prime}$ onto $\mathcal{L}$. The image of the projection is rational, lies in $\mathcal{L}$, but need not be positive semidefinite. We may ensure that $Q_{r}\in\mathbb{S}_{+}^{N}$ if the numerical matrix $Q_{n}$ is in the interior of $\mathbb{S}_{+}^{N}$ and sufficiently close to $\mathcal{L}$. More precisely, assume that $\lambda$, the smallest eigenvalue of $Q_{n}$, is greater than the distance $\delta:=\textrm{dist}(Q_{n},\mathcal{L})$. Then setting the rounding tolerance $\textrm{dist}(Q_{n},Q_{r}^{\prime})$ smaller than $\sqrt{\lambda^{2}-\delta^{2}}$ guarantees that $Q_{r}\in\mathbb{S}_{+}^{N}$; see [8, Prop. 8].

Acknowledgment

The authors would like to thank Bernd Sturmfels and the Max-Planck-Institute für Mathematik in den Naturwissenschaften in Leipzig for hosting the Macaulay2 workshop in May 2018. We thank Ilir Dema, Nidhi Kainsa and Anton Leykin, who contributed to the code. The package code contains parts of a proof-of-concept implementation of the methods in [8], originally written by Pablo Parrilo and Helfried Peyrl. Parts of this work were done while Diego Cifuentes visited the Max-Planck-Institute MiS, and while the first two authors visited ICERM supported by NSF grant No. DMS-1439786. Thomas Kahle is supported by the German Research Foundation under grant 314838170, GRK 2297 MathCoRe. Pablo Parrilo is supported in part by the National Science Foundation through NSF Grant CCF-1565235.