Sum-of-squares optimization

Given a vector ${\displaystyle c\in \mathbb {R} ^{n}}$ and polynomials ${\displaystyle a_{k,j}}$ for ${\displaystyle k=1,\dots N_{s}}$, ${\displaystyle j=0,1,\dots ,n}$, a sum-of-squares optimization problem is written as

$${\displaystyle {\begin{aligned}{\underset {u\in \mathbb {R} ^{n}}{\text{maximize}}}\quad &c^{T}u\\{\text{subject to}}\quad &a_{k,0}(x)+a_{k,1}(x)u_{1}+\cdots +a_{k,n}(x)u_{n}\in {\text{SOS}}\quad (k=1,\ldots ,N_{s}).\end{aligned}}}$$

Here "SOS" represents the class of sum-of-squares (SOS) polynomials. The quantities ${\displaystyle u\in \mathbb {R} ^{n}}$ are the decision variables. SOS programs can be converted to semidefinite programs (SDPs) using the duality of the SOS polynomial program and a relaxation for constrained polynomial optimization using positive-semidefinite matrices, see the following section.

Suppose we have an ${\displaystyle n}$-variate polynomial ${\displaystyle p(x):\mathbb {R} ^{n}\to \mathbb {R} }$ , and suppose that we would like to minimize this polynomial over a subset ${\textstyle A\subseteq \mathbb {R} ^{n}}$. Suppose furthermore that the constraints on the subset ${\textstyle A}$ can be encoded using ${\textstyle m}$ polynomial equalities of degree at most ${\displaystyle 2d}$, each of the form ${\textstyle a_{i}(x)=0}$ where ${\displaystyle a_{i}:\mathbb {R} ^{n}\to \mathbb {R} }$ is a polynomial of degree at most ${\displaystyle 2d}$. A natural, though generally non-convex program for this optimization problem is the following:

$${\displaystyle \min _{x\in \mathbb {R} ^{n}}\langle C,x^{\leq d}(x^{\leq d})^{\top }\rangle }$$

subject to:

$${\displaystyle \langle A_{i},x^{\leq d}(x^{\leq d})^{\top }\rangle =0\qquad \forall \ i\in [m],}$$

(1)

$${\displaystyle x_{\emptyset }=1,}$$

where ${\textstyle x^{\leq d}}$ is the ${\displaystyle n^{O(d)}}$-dimensional vector with one entry for every monomial in ${\displaystyle x}$ of degree at most ${\displaystyle d}$, so that for each multiset ${\displaystyle S\subset [n],|S|\leq d,}$ ${\textstyle x_{S}=\prod _{i\in S}x_{i}}$, ${\textstyle C}$ is a matrix of coefficients of the polynomial ${\textstyle p(x)}$ that we want to minimize, and ${\textstyle A_{i}}$ is a matrix of coefficients of the polynomial ${\textstyle a_{i}(x)}$ encoding the ${\displaystyle i}$-th constraint on the subset ${\displaystyle A\subset \mathbb {R} ^{n}}$. The additional, fixed constant index in our search space, ${\displaystyle x_{\emptyset }=1}$, is added for the convenience of writing the polynomials ${\textstyle p(x)}$ and ${\textstyle a_{i}(x)}$ in a matrix representation.

This program is generally non-convex, because the constraints (1) are not convex. One possible convex relaxation for this minimization problem uses semidefinite programming to replace the rank-one matrix of variables ${\displaystyle x^{\leq d}(x^{\leq d})^{\top }}$ with a positive-semidefinite matrix ${\displaystyle X}$: we index each monomial of size at most ${\displaystyle 2d}$ by a multiset ${\displaystyle S}$ of at most ${\displaystyle 2d}$ indices, ${\displaystyle S\subset [n],|S|\leq 2d}$. For each such monomial, we create a variable ${\displaystyle X_{S}}$ in the program, and we arrange the variables ${\displaystyle X_{S}}$ to form the matrix ${\textstyle X\in \mathbb {R} ^{[n]^{\leq d}\times [n]^{\leq d}}}$, where ${\displaystyle \mathbb {R} ^{[n]^{\leq d}\times [n]^{\leq d}}}$is the set of real matrices whose rows and columns are identified with multisets of elements from ${\displaystyle n}$ of size at most ${\displaystyle d}$. We then write the following semidefinite program in the variables ${\displaystyle X_{S}}$:

$${\displaystyle \min _{X\in \mathbb {R} ^{[n]^{\leq d}\times [n]^{\leq d}}}\langle C,X\rangle }$$

subject to:

$${\displaystyle \langle A_{i},X\rangle =0\qquad \forall \ i\in [m],Q}$$

$${\displaystyle X_{\emptyset }=1,}$$

$${\displaystyle X_{U\cup V}=X_{S\cup T}\qquad \forall \ U,V,S,T\subseteq [n],|U|,|V|,|S|,|T|\leq d,{\text{ and}}\ U\cup V=S\cup T,}$$

$${\displaystyle X\succeq 0,}$$

where again ${\textstyle C}$ is the matrix of coefficients of the polynomial ${\textstyle p(x)}$ that we want to minimize, and ${\textstyle A_{i}}$ is the matrix of coefficients of the polynomial ${\textstyle a_{i}(x)}$ encoding the ${\displaystyle i}$-th constraint on the subset ${\displaystyle A\subset \mathbb {R} ^{n}}$.

The third constraint ensures that the value of a monomial that appears several times within the matrix is equal throughout the matrix, and is added to make ${\displaystyle X}$ respect the symmetries present in the quadratic form ${\displaystyle x^{\leq d}(x^{\leq d})^{\top }}$.

Duality

One can take the dual of the above semidefinite program and obtain the following program:

$${\displaystyle \max _{y\in \mathbb {R} ^{m'}}y_{0},}$$

subject to:

$${\displaystyle C-y_{0}e_{\emptyset }-\sum _{i\in [m]}y_{i}A_{i}-\sum _{S\cup T=U\cup V}y_{S,T,U,V}(e_{S,T}-e_{U,V})\succeq 0.}$$

We have a variable ${\displaystyle y_{0}}$ corresponding to the constraint ${\displaystyle \langle e_{\emptyset },X\rangle =1}$ (where ${\displaystyle e_{\emptyset }}$ is the matrix with all entries zero save for the entry indexed by ${\displaystyle (\varnothing ,\varnothing )}$), a real variable ${\displaystyle y_{i}}$for each polynomial constraint ${\displaystyle \langle X,A_{i}\rangle =0\quad s.t.i\in [m],}$ and for each group of multisets ${\displaystyle S,T,U,V\subset [n],|S|,|T|,|U|,|V|\leq d,S\cup T=U\cup V}$, we have a dual variable ${\displaystyle y_{S,T,U,V}}$for the symmetry constraint ${\displaystyle \langle X,e_{S,T}-e_{U,V}\rangle =0}$. The positive-semidefiniteness constraint ensures that ${\displaystyle p(x)-y_{0}}$ is a sum-of-squares of polynomials over ${\displaystyle A\subset \mathbb {R} ^{n}}$: by a characterization of positive-semidefinite matrices, for any positive-semidefinite matrix ${\textstyle Q\in \mathbb {R} ^{m\times m}}$, we can write ${\textstyle Q=\sum _{i\in [m]}f_{i}f_{i}^{\top }}$ for vectors ${\textstyle f_{i}\in \mathbb {R} ^{m}}$. Thus for any ${\textstyle x\in A\subset \mathbb {R} ^{n}}$,

$${\displaystyle {\begin{aligned}p(x)-y_{0}&=p(x)-y_{0}-\sum _{i\in [m']}y_{i}a_{i}(x)\qquad {\text{since }}x\in A\\&=(x^{\leq d})^{\top }\left(C-y_{0}e_{\emptyset }-\sum _{i\in [m']}y_{i}A_{i}-\sum _{S\cup T=U\cup V}y_{S,T,U,V}(e_{S,T}-e_{U,V})\right)x^{\leq d}\qquad {\text{by symmetry}}\\&=(x^{\leq d})^{\top }\left(\sum _{i}f_{i}f_{i}^{\top }\right)x^{\leq d}\\&=\sum _{i}\langle x^{\leq d},f_{i}\rangle ^{2}\\&=\sum _{i}f_{i}(x)^{2},\end{aligned}}}$$

where we have identified the vectors ${\textstyle f_{i}}$ with the coefficients of a polynomial of degree at most ${\displaystyle d}$. This gives a sum-of-squares proof that the value ${\textstyle p(x)\geq y_{0}}$ over ${\displaystyle A\subset \mathbb {R} ^{n}}$.

The above can also be extended to regions ${\displaystyle A\subset \mathbb {R} ^{n}}$ defined by polynomial inequalities.

The sum-of-squares hierarchy (SOS hierarchy), also known as the Lasserre hierarchy, is a hierarchy of convex relaxations of increasing power and increasing computational cost. For each natural number ${\textstyle d\in \mathbb {N} }$ the corresponding convex relaxation is known as the ${\textstyle d}$th level or ${\textstyle d}$-th round of the SOS hierarchy. The ${\textstyle 1}$st round, when ${\textstyle d=1}$, corresponds to a basic semidefinite program, or to sum-of-squares optimization over polynomials of degree at most ${\displaystyle 2}$. To augment the basic convex program at the ${\textstyle 1}$st level of the hierarchy to ${\textstyle d}$-th level, additional variables and constraints are added to the program to have the program consider polynomials of degree at most ${\displaystyle 2d}$.

The SOS hierarchy derives its name from the fact that the value of the objective function at the ${\textstyle d}$-th level is bounded with a sum-of-squares proof using polynomials of degree at most ${\textstyle 2d}$ via the dual (see "Duality" above). Consequently, any sum-of-squares proof that uses polynomials of degree at most ${\textstyle 2d}$ can be used to bound the objective value, allowing one to prove guarantees on the tightness of the relaxation.

In conjunction with a theorem of Berg, this further implies that given sufficiently many rounds, the relaxation becomes arbitrarily tight on any fixed interval. Berg's result[5][6] states that every non-negative real polynomial within a bounded interval can be approximated within accuracy ${\textstyle \varepsilon }$ on that interval with a sum-of-squares of real polynomials of sufficiently high degree, and thus if ${\textstyle OBJ(x)}$ is the polynomial objective value as a function of the point ${\textstyle x}$, if the inequality ${\textstyle c+\varepsilon -OBJ(x)\geq 0}$ holds for all ${\textstyle x}$ in the region of interest, then there must be a sum-of-squares proof of this fact. Choosing ${\textstyle c}$ to be the minimum of the objective function over the feasible region, we have the result.

A polynomial ${\displaystyle p}$ is a sum of squares (SOS) if there exist polynomials ${\displaystyle \{f_{i}\}_{i=1}^{m}}$ such that ${\textstyle p=\sum _{i=1}^{m}f_{i}^{2}}$. For example,

$${\displaystyle p=x^{2}-4xy+7y^{2}}$$

is a sum of squares since

$${\displaystyle p=f_{1}^{2}+f_{2}^{2}}$$

where

$${\displaystyle f_{1}=(x-2y){\text{ and }}f_{2}={\sqrt {3}}y.}$$

Note that if ${\displaystyle p}$ is a sum of squares then ${\displaystyle p(x)\geq 0}$ for all ${\displaystyle x\in \mathbb {R} ^{n}}$. Detailed descriptions of polynomial SOS are available.[7][8][9]

Quadratic forms can be expressed as ${\displaystyle p(x)=x^{T}Qx}$ where ${\displaystyle Q}$ is a symmetric matrix. Similarly, polynomials of degree ≤ 2d can be expressed as

$${\displaystyle p(x)=z(x)^{\mathsf {T}}Qz(x),}$$

where the vector ${\displaystyle z}$ contains all monomials of degree ${\displaystyle \leq d}$. This is known as the Gram matrix form. An important fact is that ${\displaystyle p}$ is SOS if and only if there exists a symmetric and positive-semidefinite matrix ${\displaystyle Q}$ such that ${\displaystyle p(x)=z(x)^{\mathsf {T}}Qz(x)}$. This provides a connection between SOS polynomials and positive-semidefinite matrices.

• SOSTOOLS, licensed under the GNU GPL. The reference guide is available at arXiv:1310.4716 [math.OC], and a presentation about its internals is available here.
• CDCS-sos, a package from CDCS, an augmented Lagrangian method solver, to deal with large scale SOS programs.
• The SumOfSquares extension of JuMP for Julia.
• TSSOS for Julia, a polynomial optimization tool based on the sparsity adapted moment-SOS hierarchies.
• For the dual problem of constrained polynomial optimization, GloptiPoly for MATLAB/Octave, Ncpol2sdpa for Python and MomentOpt for Julia.