DSOS AND SDSOS OPTIMIZATION:
MORE TRACTABLE ALTERNATIVES TO
SUM OF SQUARES AND SEMIDEFINITE OPTIMIZATION∗

arXiv:1706.02586v2 [math.OC] 9 Oct 2017

AMIR ALI AHMADI† AND ANIRUDHA MAJUMDAR

‡

Abstract. In recent years, optimization theory has been greatly impacted by the advent of sum
of squares (SOS) optimization. The reliance of this technique on large-scale semidefinite programs
however, has limited the scale of problems to which it can be applied. In this paper, we introduce
DSOS and SDSOS optimization as more tractable alternatives to sum of squares optimization that
rely instead on linear and second order cone programs respectively. These are optimization problems
over certain subsets of sum of squares polynomials (or equivalently subsets of positive semidefinite
matrices), which can be of interest in general applications of semidefinite programming where scalability is a limitation. We show that some basic theorems from SOS optimization which rely on results
from real algebraic geometry are still valid for DSOS and SDSOS optimization. Furthermore, we
show with numerical experiments from diverse application areas—polynomial optimization, statistics and machine learning, derivative pricing, and control theory—that with reasonable tradeoffs
in accuracy, we can handle problems at scales that are currently far beyond the reach of sum of
squares approaches. Finally, we provide a review of recent techniques that bridge the gap between
our DSOS/SDSOS approach and the SOS approach at the expense of additional running time. The
appendix of the paper introduces an accompanying MATLAB package for DSOS and SDSOS optimization.
Key words. Sum of squares optimization, polynomial optimization, nonnegative polynomials,
semidefinite programming, linear programming, second order cone programming.
AMS subject classifications. 65K05, 90C25, 90C22, 90C05, 90C90, 12Y05, 93C85, 90C27

1. Introduction. For which values of the real coefficients c1 , c2 , c3 , is the polynomial
(1)

p(x1 , x2 ) = c1 x41 − 6x31 x2 − 4x31 + c2 x21 x22 + 10x21 + 12x1 x22 + c3 x42

nonnegative, i.e., satisfies p(x1 , x2 ) ≥ 0 for all (x1 , x2 ) ∈ R2 ? The problem of optimizing over nonnegative polynomials—of which our opening question is a toy example—is
fundamental to many problems of applied and computational modern mathematics.
In such an optimization problem, one would like to impose constraints on the coefficients of an unknown polynomial so as to make it nonnegative, either globally in
Rn (as in the above example), or on a certain basic semialgebraic set, i.e., a subset
of Rn defined by a finite number of polynomial inequalities. We will demonstrate
shortly (see Section 1.1) why optimization problems of this kind are ubiquitous in
applications and universal in the study of questions dealing in one way or another
with polynomial equations and inequalities.
Closely related to nonnegative polynomials are polynomials that are sums of
squares.
P We say that a polynomial p is a sum of squares (sos) if it can be written as
p = i qi2 for some (finite number of) polynomials qi . For example, the polynomial
∗ Submitted

to the editors on June 9, 2017.
Funding: Amir Ali Ahmadi was partially supported by a DARPA Faculty Award, a Sloan
Fellowship, an NSF CAREER Award, an AFOSR Young Investigator Program Award, and a Google
Research Award.
† Department of Operations Research and Financial Engineering, Princeton University.
(a a a@princeton.edu, http://aaa.princeton.edu/).
‡ Department of Mechanical and Aerospace Engineering, Princeton University (starting Aug. 1,
2017). (ani.majumdar@princeton.edu, http://web.stanford.edu/∼anirudha/).
1

2

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

in (1) with c1 = 13, c2 = 1, c3 = 4 is sos since it admits the decomposition
p(x1 , x2 ) = (x1 − 2x21 )2 + (3x1 + 2x22 )2 + (x1 x2 − 3x21 )2 .
The relationship between nonnegative and sum of squares polynomials has been a
classical subject of study in real algebraic geometry. For example, a result of Hilbert
from 1888 [37] states that all nonnegative bivariate polynomials of degree four are
sums of squares. It follows, as a special case, that the sets of coefficients for which
the polynomial p in (1) is nonnegative or a sum of squares in fact coincide. In
the general situation, however, while sum of squares polynomials are clearly always
nonnegative, it is not true that the converse always holds. This was shown in the
same paper of Hilbert [37], where he gave a non-constructive proof of existence of
nonnegative polynomials that are not sums of squares. Explicit examples of such
polynomials appeared many years later, starting with the work of Motzkin [59] in
the 1960s. Hilbert’s interest in this line of research is also showcased in his 17th
problem, which asks whether every nonnegative polynomial is a sum of squares of
rational functions. We refer the interested reader to an outstanding survey paper of
Reznick [74], which covers many historical aspects around Hilbert’s 17th problem,
including the affirmative solution by Artin [6], as well as several later developments.
The classical questions around nonnegative and sum of squares polynomials have
been revisited quite extensively in the past 10-15 years in different communities among
applied and computational mathematicians. The reason for this renewed interest is
twofold: (i) the discovery that many problems of modern practical interest can be cast
as optimization problems over nonnegative polynomials, and (ii) the observation that
while optimizing over nonnegative polynomials is generally NP-hard, optimization
over the set of sum of squares polynomials can be done via semidefinite programming
(SDP); see Theorem 1 in Section 2. The latter development, originally explored in the
pioneering works of Shor [79], Nesterov [61], Parrilo [65],[66], and Lasserre [44], has
led to the creation of sum of squares optimization—a computational framework, with
semidefinite programming as its underlying engine, that can tackle many fundamental
problems of real algebra and polynomial optimization.
The dependence of sum of squares approaches on semidefinite programming can
be viewed as both a strength and a weakness depending on one’s perspective. From
a computational complexity viewpoint, semidefinite programs can be solved with arbitrary accuracy in polynomial time using interior point methods (see [85] for a comprehensive survey). As a result, sum of squares techniques offer polynomial time algorithms that approximate a very broad class of NP-hard problems of interest. From
a more practical viewpoint however, SDPs are among the most expensive convex relaxations to solve. The speed and reliability of the current SDP solvers lag behind
those for other more restricted classes of convex programs (such as linear or second
order cone programs) by a wide margin. With the added complication that the SDPs
generated by sos problems are large (see Section 2), scalability has become the single
most outstanding challenge for sum of squares optimization in practice.
In this paper, we focus on the latter of the two viewpoints mentioned above
and offer alternatives to sum of squares optimization that while more conservative
in general, are significantly more scalable. Our hope is that by doing so, we expand
the use and applicability of algebraic techniques in optimization to new areas and
share its appeal with a broader audience. We call our new computational frameworks, which rely on linear and second order cone programming, DSOS and SDSOS

DSOS AND SDSOS OPTIMIZATION

3

optimization 1 . These are short for diagonally-dominant-sum-of-squares and scaleddiagonally-dominant-sum-of-squares; see Section 3 for precise definitions. While these
tools are primarily designed for sum of squares optimization, they are also applicable
to general applications of semidefinite programming where tradeoffs between scalability and performance may be desirable. In the interest of motivating our contributions
for a diverse audience, we delay a presentation of these contributions until Section 3
and start instead with a portfolio of problem areas involving nonnegative polynomials.
Any such problem area is one to which sum of squares optimization, as well as its new
DSOS and SDSOS counterparts, are directly applicable.
1.1. Why optimize over nonnegative polynomials?. We describe several
motivating applications in this section at a high level. These will be revisited later in
the paper with concrete computational examples; see Section 4.
Polynomial optimization. A polynomial optimization problem (POP) is an
optimization problem of the form
minimize p(x)
subject to x ∈ K := {x ∈ Rn | gi (x) ≥ 0, hi (x) = 0},

(2)

where p, gi , and hi are multivariate polynomials. The special case of problem (2)
where the polynomials p, gi , hi all have degree one is of course linear programming,
which can be solved very efficiently. For higher degrees, POP contains as special cases
many important problems in operations research; e.g., the optimal power flow problem
in power engineering [39], the computation of Nash equilibria in game theory [43], [67],
problems of Euclidean embedding and distance geometry [50], and a host of problems
in combinatorial optimization. We observe that if we could optimize over the set of
polynomials that are nonnegative on a basic semialgebraic set, then we could solve
the POP problem to global optimality. To see this, note that the optimal value of
problem (2) is equal to the optimal value of the following problem:
maximize γ
subject to p(x) − γ ≥ 0, ∀x ∈ K.

(3)

Here, we are trying to find the largest constant γ such that the polynomial p(x) − γ
is nonnegative on the set K; i.e., the largest lower bound on problem (2).
Combinatorial optimization. Proofs of nonnegativity of polynomials of degree
as low as four provide infeasibility certificates for all decision problems in the complexity class NP, including many well-known combinatorial optimization problems.
Consider, e.g., the simple-to-describe, NP-complete problem of PARTITION [29]:
Given a set of integers a1 , . . . , an , decide if they can be split into two sets with equal
sums. It is straightforward to see that a PARTITION instance is infeasible if and
only if the degree-four polynomial
n
n
X
X
2
2
p(x) =
(xi − 1) + (
xi ai )2
i=1

i=1

satisfies p(x) > 0, ∀x ∈ Rn . Indeed, p is by construction a sum of squares and hence
nonnegative. If it were to have a zero, each square in p would have to evaluate to zero;
but this can happen if and only if there is a binary vector x ∈ {−1, 1}n , which makes
1 We

use and recommend the pronunciation “d-sauce” and “s-d-sauce”.

4

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

Pn

i=1 xi ai = 0, giving a yes answer to the PARTITION instance. Suppose now that
for a given instance, and for some  > 0, we could prove that p(x) −  is nonnegative.
Then we would have proven that our PARTITION instance is infeasible, potentially
without trying out all 2n ways of splitting the integers into two sets.
Control systems and robotics. Numerous fundamental problems in nonlinear
dynamics and control, such as stability, robustness, collision avoidance, and controller
design can be turned into problems about finding special functions—the so-called
Lyapunov functions (see, e.g., [41])—that satisfy certain sign conditions. For example,
given a differential equation ẋ = f (x), where f : Rn → Rn , and with the origin as an
equilibrium point (i.e., satisfying f (0) = 0), consider the “region of attraction (ROA)
problem”: Determine the set of initial states in Rn from which the trajectories flow
to the origin. Lyapunov’s stability theorem (see, e.g., [41, Chap. 4]) tells us that if
we can find a (Lyapunov) function V : Rn → R, which together with its gradient ∇V
satisfies

(4)

V (x) > 0

∀x 6= 0,

and h∇V (x), f (x)i < 0

∀x ∈ {x| V (x) ≤ β, x 6= 0},

then the set {x ∈ Rn | V (x) ≤ β} is part of the region of attraction. If f is a polynomial
function (a very important case in applications [65]), and if we parameterize V as a
polynomial function, then the search for the coefficients of V satisfying the conditions
in (4) is an optimization problem over the set of nonnegative polynomials. Designing
stable equilibrium points with large regions of attraction is a fundamental problem in
control engineering and robotics; see, e.g., [40], [54].
Statistical regression with shape constraints. A problem that arises frequently in statistics is that of fitting a function to a set of data points with minimum
error, while ensuring that it meets some structural properties, such as nonnegativity,
monotonicity, or convexity [32], [34]. Requirements of this type are typically imposed
either as regularizers (to avoid overfitting), or more importantly as a result of prior
knowledge that the true function to be estimated satisfies the same structural properties. In economics, for example, a regression problem for estimating the utility of
consumers from sample measurements would come with a natural concavity requirement on the utility function. When the regression functions are polynomials, many
such structural properties lead to polynomial nonnegativity constraints. For example,
a polynomial p(x) := p(x1 , . . . , xn ) can be constrained to be concave by requiring its
Hessian matrix ∇2 p(x), which is an n × n symmetric matrix with polynomial entries,
to be negative semidefinite for all x. This condition is easily seen to be equivalent to
the polynomial −y T H(x)y in the 2n variables x := (x1 , . . . , xn ) and y := (y1 , . . . , yn )
being nonnegative.
Probability bounds given moments. The well-known inequalities of Markov
and Chebyshev in probability theory bound the probability that a univariate random
variable takes values in a subset of the real line given the first, or the first and second
moments of its distribution [10]. Consider a vast generalization of this problem where
one is given a finite number of moments {µk1 ,...,kn := E[X1k1 . . . Xnkn ], k1 +· · · , kn ≤ d}
of a multivariate multivariate random variable X := (X1 , . . . , Xn ) and is interested in
bounding the probability of an event X ∈ S, where S is a general basic semialgebraic
subset of the sample space Ω. A sharp upper bound on this probability can be
obtained by solving the following optimization problem for the coefficients ck1 ...kn of
an n-variate degree-d polynomial p (see [14], [70]):

DSOS AND SDSOS OPTIMIZATION

(5)

5

P
minimize E[p(X)] = ck1 ...kn µk1 ...kn
subject to p(x) ≥ 1, ∀x ∈ S
p(x) ≥ 0, ∀x ∈ Ω.

It is easy to see that any feasible solution of problem (5) gives an upper bound
on the probability of the event X ∈ S as the polynomial p is by construction placed
above the indicator function of the set S and the expected value of this indicator
function is precisely the probability of the event X ∈ S. Note that the constraints in
(5) are polynomial nonnegativity conditions.
Other applications. The application areas outlined above are only a few examples on a long list of problems that can be tackled by optimizing over the set of
nonnegative polynomials. Other interesting problems on this list include: computation of equilibria in games with continuous strategy spaces [67], distinguishing separable and entangled states in quantum information theory [25], computing bounds
on sphere packing problems in geometry [7], software verification [77], robust and
stochastic optimization [12], filter design [82], and automated theorem proving [36].
A great reference for many applications in this area, as well as the related theoretical
and computational background, is a recent volume edited by Blekherman, Parrilo,
and Thomas [15].
Organization of the paper. The remainder of this paper is structured as follows. In Section 2, we briefly review the sum of squares approach and its relation to
semidefinite programming for the general reader. We also comment on its challenges
regarding scalability and prior work that has tackled this issue. The familiar reader
can skip to Section 3, where our contributions begin. In Section 3.1, we introduce the
cone of dsos and sdsos polynomials and demonstrate that we can optimize over them
using linear and second order cone programming respectively. Section 3.2 introduces
hierarchies of cones based on dsos and sdsos polynomials that can better approximate
the cone of nonnegative polynomials. Asymptotic guarantees on these cones based on
classical Positivstellensätze are also presented. Section 4 presents numerical experiments that use our approach on large-scale examples from polynomial optimization
(Section 4.1), combinatorial optimization (Section 4.2), statistics and machine learning (Sections 4.3 and 4.5), derivative pricing (Section 4.4), and control theory and
robotics (Section 4.6). Section 5 provides a review of recent techniques for bridging
the gap between the DSOS/SDSOS approach and the SOS approach at the expense
of additional computational costs. Section 6 concludes the paper. The Appendix to
the paper introduces iSOS (“inside SOS”), an accompanying MATLAB package for
DSOS and SDSOS optimization written using the SPOT toolbox [57].
2. Review of the semidefinite programming-based approach and computational considerations. As mentioned earlier, sum of squares optimization handles a (global) nonnegativity constraint on a polynomial p by replacing it with the
stronger requirement that p be a sum of squares. The situation where p is only
constrained to be nonnegative on a certain basic semialgebraic set2
S := {x ∈ Rn | g1 (x) ≥ 0, . . . , gm (x) ≥ 0}
2 In this formulation, we have avoided equality constraints for simplicity. Obviously, there is no
loss of generality in doing this as an equality constraint h(x) = 0 can be imposed by the pair of
inequality constraints h(x) ≥ 0, −h(x) ≥ 0.

6

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

can also be handled with the help of appropriate sum of squares multipliers. For
example, if we succeed in finding sos polynomials s0 , s1 , . . . , sm , such that

(6)

p(x) = s0 (x) +

m
X

si (x)gi (x),

i=1

then we have found a certificate of nonnegativity of p on the set S. Indeed, if we evaluate the above expression at any x ∈ S, nonnegativity of the polynomials s0 , s1 . . . , sm
imply that p(x) ≥ 0. A Positivstellensatz from real algebraic geometry due to Putinar [72] states that if the set S satisfies the so-called Archimedean property, a property
only slightly stronger than compactness (see, e.g., [47] for a precise definition), then
every polynomial positive on S has a representation of the type (6), for some sos
polynomials s0 , s1 , . . . , sm of high enough degree (see also [62] for degree bounds).
Even with no qualifications whatsoever regarding the set S, there are other Positivstellensätze (e.g., due to Stengle [80]) that certify nonnegativity of a polynomial
on a basic semialgebraic set using sos polynomials. These certificates are only slightly
more complicated than (6) and involve sos multipliers associated with products among
polynomials gi that define S [66]. A great reference for the interested reader is the
survey paper by Laurent [47].
The computational advantage of a certificate of (global or local) nonnegativity
via sum of squares polynomials is that it can be automatically found by semidefinite
programming. The following well-known theorem establishes the link between sos
polynomials and SDP. Let us denote the set of n × n symmetric matrices by Sn .
Recall that a matrix A ∈ Sn is positive semidefinite (psd) if xT Ax ≥ 0, ∀x ∈ Rn , and
that semidefinite programming is the problem of optimizing a linear function over psd
matrices subject to affine inequalities on their entries [85]. We denote the positive
semidefiniteness of a matrix A with the standard notation A  0.
Theorem 1 (see, e.g., [65],[66]). A multivariate polynomial p := p(x) in n
variables and of degree 2d is a sum of squares if and only if there exists a symmetric
matrix Q (often called the Gram matrix) such that
(7)

p(x) =
Q 

z T Qz,
0,

where z is the vector of monomials of degree up to d:
z = [1, x1 , x2 , . . . , xn , x1 x2 , . . . , xdn ].
Searching for a matrix Q satisfying the positive semidefiniteness constraint, as well
as the linear equality constraints coming from (7), amounts to solving a semidefinite
program. The size of the matrix Q in this theorem is

 

n+d
n+d
×
,
d
d
which approximately equals nd × nd . While this number is polynomial in n for fixed
d, it can grow rather quickly even for low-degree polynomials. For example, a degree4 polynomial (d = 2) in 50 variables has 316251 coefficients and its Gram matrix,
which would need to be positive semidefinite, contains 879801 decision variables. A
semidefinite constraint of this size is quite expensive, and in fact a problem with 50

DSOS AND SDSOS OPTIMIZATION

7

variables is far beyond the current realm of possibilities in SOS optimization. In
the absence of problem structure, sum of squares problems involving degree-4 or 6
polynomials are currently limited, roughly speaking, to a handful or a dozen variables.
There have been many contributions already to improvements in scalability of
sum of squares techniques. One approach has been to develop systematic techniques
for taking advantage of problem structure, such as sparsity or symmetry of the underlying polynomials, to reduce the size of the SDPs [30], [84], [23], [76]. These techniques
have proven to be very useful as problems arising in practice sometimes do come with
a lot structure. Another approach which holds promise has been to design customized
solvers for special classes of SOS programs and avoid restoring to an off-the-shelf interior point solver. Examples in this direction include the works in [11], [63], [86], [49].
There has also been recent work by Lasserre et al. that increases scalability of sum of
squares optimization problems at the cost of accuracy of the solutions obtained. This
is done by bounding the size of the largest SDP constraint appearing in the sos formulation, and leads to what the authors call the BSOS (bounded SOS) hierarchy [46].
The approach we take in this paper for enhancing scalability is orthogonal to the
ones mentioned above (and can potentially later be combined with them). We propose
to eschew sum of squares decompositions to begin with. In our view, almost all application areas of this field use sum of squares decompositions not for their own sake,
but because they provide a means to polynomial nonnegativity. Hence, this paper is
motivated by a natural question: Can we give other sufficient conditions for polynomial nonnegativity that are perhaps stronger than a sum of squares decompostion,
but cheaper to work with?
In the next section, we identify some subsets of the cone of nonnegative polynomials that one can optimize over using linear programming (LP) and second order
cone programming (SOCP) [51], [5]. Not only are these much more efficiently solvable classes of convex programs, but they are also superior to semidefinite programs
in terms of numerical conditioning. In addition, working with these classes of convex programs allows us to take advantage of high-performance LP and SOCP solvers
(e.g., CPLEX [20], Mosek [58]) that have been matured over several decades because
of industry applications. We remark that while there have been other approaches
to produce LP hierarchies for polynomial optimization problems (e.g., based on the
Krivine-Stengle certificates of positivity [42, 80, 45]), these LPs, though theoretically
significant, are typically quite weak in practice and often numerically ill-conditioned
[46].
3. DSOS and SDSOS Optimization. We start with some basic definitions.
A monomial m : Rn → R in variables x := (x1 , . . . , xn ) is a function of the form
αn
1
m(x) = xα
1 . . . xn , where each αi is a nonnegative integer. The degree of such
n
a monomial is by definition α1 + · · · + αP
n . A polynomial p : R → R is a finite
α1
αn
linear combination of monomials p(x) = α1 ...αn cα1 ...αn x1 . . . xn . The degree of
a polynomial is the maximum degree of its monomials. A form or a homogeneous
polynomial is a polynomial whose monomials all have the same degree. For any
polynomial p := p(x1 , . . . , xn ) of degree d, one can define its homogenized version
ph by introducing a new variable y and defining ph (x1 , . . . , xn , y) = y d p( xy ). One
can reverse this operation (dehomogenize) by setting the homogenizing variable y
equal to one: p(x1 , . . . , xn ) = ph (x1 , . . . , xn , 1). The properties of being nonnegative
and a sum of squares (as defined earlier) are preserved under homogenization and
dehomogenization [74]. We say that a form f is positive definite if f (x) > 0 for all
x 6= 0. We denote by P SDn,2d and SOSn,2d the set of nonnegative (a.k.a. positive

8

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

semidefinite) and sum of squares polynomials in n variables and degree 2d respectively.
(Note that odd-degree polynomials can never be nonnegative.) Both these sets form
n+2d
proper (i.e., convex, closed, pointed, and solid) cones in R( 2d ) with the obvious
inclusion relationship SOSn,2d ⊆ P SDn,2d .
The basic idea behind our approach is to approximate SOSn,2d from the inside
with new sets that are more tractable for optimization purposes. Towards this goal,
one may think of several natural sufficient conditions for a polynomial to be a sum of
squares. For example, the following may be natural first candidates:
• The set P
of polynomials that are sums of 4-th powers of polynomials:
{p| p =
qi4 },
• The set of polynomials that are a sum of a constant number of squares of
polynomials; e.g., a sum of three squares: {p| p = q12 + q22 + q32 }.
Even though both of these sets clearly reside inside the set of sos polynomials,
they are not any easier to optimize over. In fact, they are much harder: testing
whether a polynomial is a sum of 4-th powers is NP-hard already for quartics [38]
(in fact, the cone of 4-th powers of linear forms is dual to the cone of nonnegative
quartic forms [75]) and optimizing over polynomials that are sums of three squares
is intractable (as this task even for quadratics subsumes the NP-hard problem of
positive semidefinite matrix completion with a rank constraint [48]). These examples
illustrate that in general, an inclusion relationship does not imply anything with
respect to optimization complexity. Hence, we need to take some care when choosing
which subsets of SOSn,2d to work with. On the one hand, these subsets have to
be “big enough” to be useful in practice; on the other hand, they should be more
tractable to optimize over.
3.1. The cone of dsos and sdsos polynomials. We now describe two interesting cones inside SOSn,2d that lend themselves to linear and second order cone
representations and are hence more tractable for optimization purposes. These cones
will also form the building blocks of some more elaborate cones (with improved performance) that we will present in Subsection 3.2 and Section 5.
Definition 2. A polynomial p := p(x) is diagonally-dominant-sum-of-squares
(dsos) if it can be written as
X
X
X
+
−
(8)
p(x) =
αi m2i (x) +
βij
(mi (x) + mj (x))2 +
βij
(mi (x) − mj (x))2 ,
i

i,j

i,j

+
−
for some monomials mi (x), mj (x) and some nonnegative scalars αi , βij
, βij
. We denote the set of polynomials in n variables and degree 2d that are dsos by DSOSn,2d .

Definition 3. A polynomial p := p(x) is scaled-diagonally-dominant-sum-ofsquares (sdsos) if it can be written as
X
X
X
+
+
−
−
(9) p(x) =
αi m2i (x) +
(β̂ij
mi (x) + β̃ij
mj (x))2 +
(β̂ij
mi (x) − β̃ij
mj (x))2 ,
i

i,j

i,j

+
+
−
−
for some monomials mi (x), mj (x) and some scalars αi , β̂ij
, β̃ij
, β̂ij
, β̃ij
with αi ≥ 0.
We denote the set of polynomials in n variables and degree 2d that are sdsos by
SDSOSn,2d .

From the definition, the following inclusion relationships should be clear:
DSOSn,2d ⊆ SDSOSn,2d ⊆ SOSn,2d ⊆ P OSn,2d .

DSOS AND SDSOS OPTIMIZATION

9

Fig. 1. A comparison of the set of dsos/sdsos/sos/psd polynomials on a parametric family of
bivariate quartics given in (10).

In general, all these containment relationships are strict. Figure 1 shows these sets
for a family of bivariate quartics parameterized by two coefficients a, b:3

(10)

1
1
p(x1 , x2 ) = x41 + x42 + ax31 x2 + (1 − a − b)x21 x22 + 2bx1 x32 .
2
2

Some elementary remarks about the two definitions above are in order. It is not
hard to see that if p has degree 2d, the monomials mi , mj in (8) or (9) never need
to have degree higher than d. In the special case where p is homogeneous of degree
2d (i.e., a form), the monomials mi , mj can be taken to have degree exactly d. We
also note that decompositions given in (8) or (9) are not unique; there can even be
infinitely many such decompositions but the definition requires just one.
Our terminology in Definitions 2 and 3 comes from a connection (which we will
soon establish) to the following classes of symmetric matrices.
Definition 4. A symmetric matrix A = (aij ) is diagonally dominant (dd) if
X
aii ≥
|aij |
j6=i

for all i. A symmetric matrix A is scaled diagonally dominant (sdd) if there exists a
diagonal matrix D, with positive diagonal entries, such that DAD is dd.4 We denote
the set of n × n dd and sdd matrices with DDn and SDDn respectively.
It follows from Gershgorin’s circle theorem [31] that diagonally dominant matrices
are positive semidefinite. This implies that scaled diagonally dominant matrices are
also positive semidefinite as the eigenvalues of DAD have the same sign as those of
A. Hence, if we denote the set of n × n positive semidefinite matrices by Pn , the
following inclusion relationships are evident:
3 It is well known that all psd bivariate forms are sos [74] and hence the outer set exactly characterizes all values of a and b for which this polynomial is nonnegative.
4 Checking whether a given matrix A = (a ) is sdd can be done by solving a linear program since
ij
the definition is equivalent to existence of an element-wise positive vector d such that
X
aii di ≥
|aij |dj , ∀i.
j6=i

10

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

Fig. 2. A section of the cone of 5 × 5 diagonally dominant, scaled diagonally dominant, and
positive semidefinite matrices. Optimization over these sets can respectively be done by LP, SOCP,
and SDP.

DDn ⊆ SDDn ⊆ Pn .
These containments are strict for n > 2. In Figure 2, we illustrate the set of x and y
for which the matrix
I + xA + yB
is diagonally dominant, scaled diagonally dominant, and positive semidefinite. Here,
I is the identity matrix and the 5 × 5 symmetric matrices A and B were generated
randomly with iid entries sampled from the standard normal distribution. Our interest
in these inner approximations to the set of psd matrices stems from the fact that
optimization over them can be done by linear and second order cone programming
respectively (see Theorem 10 below, which is more general). For now, let us relate
these matrices back to dsos and sdsos polynomials.
Theorem 5. A polynomial p of degree 2d is dsos if and only if it admits a representation as p(x) = z T (x)Qz(x), where z(x) is the standard monomial vector of
degree ≤ d and Q is a dd matrix.
The following lemma by Barker and Carlson provides an extreme ray characterization of the set of diagonally dominant matrices and will be used in the proof of the
above theorem.
2

Lemma 6 (Barker and Carlson [8]). Let {vi }ni=1 be the set of all nonzero vectors
in Rn with at most 2 nonzero components, each equal to ±1. Then, a symmetric n × n
matrix M is diagonally dominant if and only if it can be written as
2

M=

n
X

ηi vi viT ,

i=1

for some ηi ≥ 0.
Proof of Theorem 5. Suppose first that p(x) = z T (x)Qz(x) with Q diagonally
dominant. Lemma 6 implies that
2

p(x) =

n
X
i=1

2

ηi z

T

(x)vi viT z(x)



=

n
X
i=1

ηi (viT z(x))2 ,

DSOS AND SDSOS OPTIMIZATION

11

where ηi ≥ 0 and {vi } is the set of all nonzero vectors in Rn with at most 2 nonzero
components, each equal to ±1. In this sum, the vectors vi that have only one nonzero
entry lead to the first sum in the dsos expansion (8). The vectors vi with two 1s or
two −1s lead to the second sum, and those with one 1 and one −1 lead to the third.
For the reverse direction, suppose p is dsos, i.e., has the representation in (8). We
will show that
X
p(x) = z T (x)
Qk z(x),
k

where each Qk is dd and corresponds to a single term in the expansion (8). As the
sum of dd matrices is dd, this would establish the proof. For each term of the form
αi m2i (x) we can take Qk to be a matrix of all zeros except for a single diagonal entry
corresponding to the monomial mi in z(x), with this entry equaling αi . For each term
+
of the form βij
(mi (x) + mj (x))2 , we can take Qk to be a matrix of all zeros except for
four entires corresponding to the monomials mi (x) and mj (x) in z(x). All these four
+
−
entries will be set to βij
. Similarly, for each term of the form βij
(mi (x) − mj (x))2 ,
we can take Qk to be a matrix of all zeros except for four entires corresponding to
the monomials mi and mj (x) in z(x). In this 2 × 2 submatrix, the diagonal elements
−
−
will be βij
and the off-diagonal elements will be −βij
. Clearly, all the Qk matrices
we have constructed are dd.
Theorem 7. A polynomial p of degree 2d is sdsos if and only if it admits a
representation as p(x) = z T (x)Qz(x), where z(x) is the standard monomial vector of
degree ≤ d and Q is a sdd matrix.
Proof. Suppose first that p(x) = z T (x)Qz(x) with Q scaled diagonally dominant.
Then, there exists a diagonal matrix D, with positive diagonal entries, such that DQD
is diagonally dominant and
p(x) = (D−1 z(x))T DQD(D−1 z(x)).
Now an argument identical to that in the first part of the proof of Theorem 5 (after
replacing z(x) with D−1 z(x)) gives the desired sdsos representation in (9).
For the reverse direction, suppose p is sdsos, i.e., has the representation in (9).
We show that
X
p(x) = z T (x)
Qk z(x),
k

where each Qk is sdd and corresponds to a single term in the expansion (9). As the
sum of sdd matrices is sdd,5 this would establish the proof. Indeed, for each term
of the form αi m2i (x) we can take Qk (once again) to be a matrix of all zeros except
for a single diagonal entry corresponding to the monomial mi in z(x), with this entry
equaling αi .
+
+
For each term of the form (β̂ij
mi (x) + β̃ij
mj (x))2 , we can take Qk to be a matrix
of all zeros except for four entires corresponding
to the monomials mi (x) and mj (x)
!

+ 
β̂ij
+
+
. Similarly, for each term of
in z(x). This 2 × 2 block will then be
β̂ij
β̃ij
+
β̃ij
−
−
the form (β̂ij
mi (x) − β̃ij
mj (x))2 , we can take Qk to be a matrix of all zeros except for
four entires corresponding to the monomials mi and mj (x) in z(x). This 2 × 2 block
5 This

claim is not obvious from the definition but is apparent from Lemma 9 below which implies
that SDDn is a convex cone.

12

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR
−
β̂ij
−
−β̃ij

!


−
. It remains to show that a 2 × 2 rank-1 positive
−β̃ij
 

a
6
a b be such a matrix. Observe that
semidefinite matrix is sdd. Let A =
b
 1
 
1

0
0
1 1
a
a
A
=
1 1
0 1b
0 1b
will then be



−
β̂ij

is dd and hence A is by definition sdd. (Note that if a or b are zero, then A is already
dd.)
The following characterizations of sdd matrices will be important for us.
Theorem 8 (see theorems 8 and 9 by Boman et al. [16]). A symmetric matrix
Q is sdd if and only if it has “factor width” at most 2; i.e., a factorization Q = V V T ,
where each column of V has at most two nonzero entries.
Lemma 9. A symmetric n × n matrix Q is sdd if and only if it can be expressed
as
(11)

Q=

X

M ij ,

i<j

where each M ij is an n × n matrix with zeros everywhere except
four entries
 for

(M ij )ii (M ij )ij
ij
ij
ij
ij
(M )ii , (M )ij , (M )ji , (M )jj , which make the 2 × 2 matrix
(M ij )ji (M ij )jj
symmetric and positive semidefinite.
Proof. This is almost
from Theorem 8: If Q = V V T for some n × k
Pk immediate
T
matrix V , then Q = l=1 vl vl , where vl denotes the l-th column of V . Since each
column vl has at most two nonzero entries, say in positions i and j, each matrix vl vlT
will be exactly of the desired form M ij in the statement of the lemma. Conversely,
T
if Q can be written as in (11), then we can write each M ij as M ij = wij,1 wij,1
+
T
wij,2 wij,2 , where the vectors wij,1 and wij,2 have at most two nonzero entries. If we
then construct a matrix V which has the collection of the vectors wij,1 and wij,2 as
columns, we will have Q = V V T .
Theorem 10. For any fixed d, optimization over DSOSn,2d (resp. SDSOSn,2d )
can be done with a linear program (resp. second order cone program) of size polynomial
in n.
Proof. We will use the characterizations of dsos/sdsos polynomials given in theorems 5 and 7. In both cases, the equality p(x) = z T (x)Qz(x), ∀x can be imposed by a
finite set of linear equations in the coefficients of p and the entries of Q (these match
the coefficients of p(x) with those of z T (x)Qz(x)). The constraint that Q be dd can
be imposed, e.g., by a set of linear inequalities
P
Qii ≥ j6=i zij , ∀i,
−zij ≤ Qij ≤ zij , ∀i, j, i 6= j
in variables Qij and zij . This gives a linear program.
The constraint that Q be sdd can be imposed via Lemma 9 by a set of equality constraints that enforce equation (11). The constraint that each 2 × 2 matrix
6 This

also implies that any 2 × 2 positive semidefinite matrix is sdd.

DSOS AND SDSOS OPTIMIZATION

13


(M ij )ii (M ij )ij
be psd is a “rotated quadratic cone” constraint and can be im(M ij )ji (M ij )jj
posed using SOCP [5, 51]:


2(M ij )ij
ij
ij
≤ (M ij )ii + (M ij )jj .
(M )ii + (M )jj ≥ 0,
(M ij )ii − (M ij )jj


In both cases, the final
 LP and SOCP are of polynomial size because the size of the
Gram matrix is n+d
d , which is polynomial in n for fixed d.
We have written publicly-available code that automates the process of generating
LPs (resp. SOCPs) from dsos (resp. sdsos) constraints on a polynomial. More
information about this can be found in the appendix.
The ability to replace SDPs with LPs and SOCPs is what results in significant
speedups in our numerical experiments (see Section 4). For the special case where the
polynomials involved are quadratic forms, we get a systematic way of inner approximating semidefinite programs. A generic SDP that minimizes Tr(CX) subject to the
constraints Tr(Ai X) = bi , i = 1, . . . , m, and X  0, can be replaced by
• a diagonally dominant program (DDP); i.e., a problem of the form
(12)

minimize
X∈Sn

s.t.

Tr(CX)
Tr(Ai X) = bi ,
X

i = 1, . . . , m,

dd,

which is an LP, or
• a scaled diagonally dominant program (SDDP); i.e., a problem of the form
(13)

minimize
X∈Sn

s.t.

Tr(CX)
Tr(Ai X) = bi ,
X

i = 1, . . . , m,

sdd,

which is an SOCP.
DDP and SDDP fit nicely into the framework of conic programming as they are
optimization problems over the intersection of an affine subspace with a proper cone
(see, e.g., [9, Chap. 2] for a review of conic programming). In sections 4.4 and 4.5,
we show applications of this idea to problems in finance and statistics.
3.2. The cone of r-dsos and r-sdsos polynomials and asymptotic guarantees. We now present a hierarchy of cones based on the notions of dsos and sdsos
polynomials that can be used to better approximate the cone of nonnegative polynomials.
Definition 11. For an integer r ≥ 0, we say that a polynomial p := p(x1 , . . . , xn )
is r-dsos (resp. r-sdsos) if
n
X
p(x) · (
x2i )r
i=1

is dsos (resp. sdsos). We denote the set of polynomials in n variables and degree 2d
that are r-dsos (resp. r-sdsos) by rDSOSn,2d (resp. rDSOSn,2d ).
Note that
Pnfor r = 0 we recover our dsos/sdsos definitions. Moreover, because the
multiplier ( i=1 x2i )r is nonnegative, we see that for any r, the property of being

14

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

(a) The LP-based r-dsos hierarchy.

(b) The SOCP-based r-sdsos hierarchy.

Fig. 3. The improvement obtained by going one level up in the r-dsos/r-sdsos hierarchies for
the family of bivariate quartics in (10).

r-dsos or r-sdsos is a sufficient condition for nonnegativity:
rDSOSn,2d ⊆ rSDSOSn,2d ⊆ P SDn,2d .
Optimizing over r-dsos (resp. r-sdsos) polynomials is still an LP (resp. SOCP).
The proof of the following theorem is identical to that of Theorem 10 and hence
omitted.
Theorem 12. For any fixed d and r, optimization over the set rDSOSn,d (resp.
rSDSOSn,d ) can be done with linear programming (resp. second order cone programming) of size polynomial in n.
If we revisit the parametric family of bivariate quartics in (10), the improvements
obtained by going one level higher in the hierarchy are illustrated in Figure 3. Interestingly, for r ≥ 1, r-dsos and r-sdsos tests for nonnegativity can sometimes outperform
the sos test. The following two examples illustrate this.
Example 3.1. Consider the Motzkin polynomial [59], M (x) = x41 x22 + x21 x42 −
which historically is the first known example of a nonnegative polynomial
that is not a sum of squares. One can give an LP-based nonnegativity certificate of
this polynomial by showing that M ∈ 2DSOS3,6 . Hence, 2DSOS3,6 * SOS3,6 .
3x21 x22 x23 +x63 ,

Example 3.2. Consider the polynomial p(x) = x41 x22 + x42 x23 + x43 x21 − 3x21 x22 x23 .
Once again this polynomial is nonnegative but not a sum of squares [74]. However,
by solving an LP feasibility problem, one can check that p ∈ 1DSOS3,6 , which proves
that p is nonnegative, and that 1DSOS3,6 * SOS3,6 .
It is natural to ask whether every nonnegative polynomial is r-dsos (or r-sdsos)
for some r? The following theorems deal with this question and provide asymptotic
guarantees on r-dsos (and hence r-sdsos) hierarchies. Their proofs rely on classical
Positivstellensätze from real algebraic geometry.
Theorem 13. Let p be an even7 positive definite form. Then, there exists an
integer r for which p is r-dsos (and hence r-sdsos).
Proof. A well-known theorem of Pólya [69] states that ifP
a form f (x1 , . . . , xn ) is
positive on the simplex (i.e., the set ∆n := {x ∈ Rn |xi ≥ 0, i xi = 1}), then there
exists an integer r such that all coefficients of
(x1 + · · · + xn )r f (x1 , . . . , xn )
7 An

even polynomial is a polynomial whose individual variables are raised only to even degrees.

DSOS AND SDSOS OPTIMIZATION

15

are positive. Under the assumptions of the theorem, this implies that there is an r
√
√
for which the form p( x1 , . . . , xn )(x1 + · · · + xn )r and hence the form
p(x1 , . . . , xn )(x21 + · · · + x2n )
P
have positive coefficients. But this means that p(x)( i x2i ) is a nonnegative weighted
sum of squared monomials
and therefore clearly dsos (see (8)). (In a Gram matrix
P
representation p(x)( i x2i ) = z T (x)Qz(x), the argument we just gave shows that we
can take Q to be diagonal.)
We remark that even forms already constitute a very interesting class of polynomials since they include, e.g., all polynomials coming from copositive programming;
see Subsection 4.2. In fact, any NP-complete problem can be reduced in polynomial
time to the problem of checking whether a degree-4 even form is nonnegative [60]. An
application of this idea to the independent set problem in combinatorial optimization
is presented in Subsection 4.2.
The next proposition shows that the evenness assumption cannot be removed
from Theorem 13.
Proposition 14. For any 0 < a < 1, the quadratic form
p(x1 , x2 , x3 ) = (x1 + x2 + x3 )2 + a(x21 + x22 + x23 )

(14)

is positive definite but not r-sdsos for any r.
Proof. Positive definiteness is immediate from having a > 0. Let us show that
the quadratic form (x1 + x2 + x3 )2 is not r-sdsos for any r. Consider the polynomial
q(x) := (x1 + x2 + x3 )2 (x21 + x22 + x23 )r ,
which is a form of degree 2r + 2. Note that the coefficient of x2r+2
in q is 1. Observe
1
2r+1
also that q has the following two monomials: 2x2r+1
x
and
2x
x3 . Suppose for
2
1
1
the sake of contradiction that q was sdsos. Then it could be written as
X
X
X
+
+
−
−
q(x) =
αi m2i (x) +
(β̂ij
mi (x) + β̃ij
mj (x))2 +
(β̂ij
mi (x) − β̃ij
mj (x))2 ,
i

i,j

i,j

where all monomials mi (x) and mj (x) are of degree r+1. So to produce the monomials
2x12r+1 x2 and 2x2r+1
x3 in q, the right hand side above must include the terms
1
(xr1 x2 + xr+1
)2 + (xr1 x3 + xr+1
)2 .
1
1
But this also produces 2x2r+2
, which cannot be canceled. This contradicts the fact
1
that the coefficient of x2r+2
in q is 1.
1
To see that for any a ∈ (0, 1) the polynomial p in the statement of the proposition
cannot be r-sdsos either, one can repeat the exact same argument and observe that
the coefficient of x2r+2
would still not match.
1
We remark that the form in (14) can be proved to be positive the improvements
presented in Section 5, in particular with a “factor width 3” proof system from Subsection 5.3. Alternatively, one can work with the next theorem, which removes the
assumption of evenness from Theorem 13 at the cost of doubling the number of variables and the degree.8
8 Instead of Theorem 15, an earlier version of this paper was quoting a result by Habicht [33]
from [73, p. 1], [71, p. 8]. We thank Claus Scheiderer for informing us that the statement of
Habicht [33] in German was being misquoted by us and the above references.

16

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

Theorem 15 (see Section 4 of [2]). An n-variate form p : = p(x) of degree 2d is
positive definite if and only if there exists a positive integer r that makes the following
form r-dsos:
n
n
1 X 2
1 X 4
2 2 d
√
√
p(v − w ) −
( (v − wi ) ) +
( (v + wi4 ))d .
r i=1 i
2 r i=1 i
2

2

In [2, Section 4], this theorem is used to construct LP and SOCP-based converging
hierarchies of lower bounds for the general polynomial optimization problem (2) when
the feasible set is compact.
4. Numerical examples and applications. In this section, we consider several numerical examples that highlight the scalability of the approach presented in
this paper. Comparisons to existing methods on examples considered in the literature are provided whenever possible. A software package written using the Systems
Polynomial Optimization Toolbox (SPOT) [57] includes a complete implementation
of the presented methods and is available online9 . The toolbox features efficient polynomial algebra and allows us to setup the large-scale LPs and SOCPs arising from
our examples. The Appendix provides a brief introduction to the software package
and the associated functionality required for setting up the DSOS/SDSOS programs
considered in this paper. The LPs and SOCPs resulting from our programs are solved
with the MOSEK solver on a machine with 4 processors and a clock speed of 3.4
GHz and 16GB RAM. Running times for comparisons to SDP-based approaches are
presented for the MOSEK SDP solver10 .
4.1. Lower bounds on polynomial optimization problems. An important
application of sum of squares optimization is to obtain lower bounds on polynomial
optimization problems. In this subsection, we consider the particular problem of
minimizing a homogeneous polynomial p(x) of degree 2d on the unit sphere (i.e., the
set {x ∈ Rn | xT x = 1}). This well-studied problem is strongly NP-hard even when
d = 2. Its optimal value is easily seen to be equivalent to the optimal value of the
following problem:

(15)

maximize
γ

s.t.

γ
p(x) − γ(xT x)d ≥ 0, ∀x ∈ Rn .

By replacing the nonnegativity condition with a SOS/DSOS/SDSOS constraint,
the resulting problem becomes an SDP/LP/SOCP. Tables 1 and 2 compare optimal
values and running times respectively on problems of increasing size. We restrict
ourselves to quartic forms and vary the number of variables n between 5 and 70.
Table 1 compares the optimal values for the DSOS, SDSOS, 1-DSOS, 1-SDSOS and
SOS relaxations of problem (15) on quartic forms whose coefficients are fully dense
and drawn independently from the standard normal distribution (we consider a single
random instance for each n). We note that the optimal values presented here are
representative of numerical experiments we performed for a broader set of random
9 Link

to spotless isos software package:
https://github.com/anirudhamajumdar/spotless/tree/spotless isos
10 We note that the MOSEK solver is much faster on our problems than the more commonly used
SDP solver SeDuMi [81].

17

DSOS AND SDSOS OPTIMIZATION
n = 10 n = 15
n = 20 n = 25 n = 30 n = 40 n = 50 n = 60
DSOS
-5.31
-10.96
-18.012 -26.45
-36.85
-62.30 -94.26 -133.02
SDSOS
-5.05
-10.43
-17.33
-25.79
-36.04
-61.25 -93.22 -131.64
1-DSOS
-4.96
-9.22
-15.72
-23.58
NA
NA
NA
NA
1-SDSOS -4.21
-8.97
-15.29
-23.14
NA
NA
NA
NA
SOS
-1.92
-3.26
-3.58
-3.71
NA
NA
NA
NA
BARON -175.41 -1079.89 -5287.88
-28546.1
Table 1
Comparison of lower bounds on the minimum of a quartic form on the sphere
number of variables.

n = 70
-178.23
-176.89
NA
NA
NA
for varying

n = 10 n = 15 n = 20 n = 25 n = 30 n = 40 n = 50 n = 60 n = 70
0.30
0.38
0.74
15.51
7.88
10.68
25.99
58.10 342.76
0.27
0.53
1.06
8.72
5.65
18.66
47.90 109.54 295.30
0.92
6.26
37.98
369.08
∞
∞
∞
∞
∞
1.53
14.39
82.30 538.5389
∞
∞
∞
∞
∞
0.24
5.60
82.22
1068.66
∞
∞
∞
∞
∞
0.35
0.62
3.69
Table 2
Comparison of running times (in seconds) averaged over 10 instances for lower bounding a
quartic form on the sphere for varying number of variables.
DSOS
SDSOS
1-DSOS
1-SDSOS
SOS
BARON

instances (including coefficients drawn from a uniform distribution instead of the
normal distribution). Also, in all our experiments, we could verify that the SOS
bound actually coincides with the true global optimal value of the problem. This was
done by finding a point on the sphere that produced a matching upper bound.
As the tables illustrate, both DSOS and SDSOS scale up to problems of dimension
n = 70. This is well beyond the size of problems handled by SOS programming,
which is unable to deal with problems with n larger than 25 due to memory (RAM)
constraints. While the price we pay for this increased scalability is suboptimality, this
is tolerable in many applications (as we show in other examples in this section).
Table 1 also compares the approach presented in this paper to the lower bound
obtained at the root node of the branch-and-bound procedure in the popular global
optimization package BARON [78].11 This comparison was only done up to dimension
30 since beyond this we ran into difficulties passing the large polynomials to BARON’s
user interface.
Table 2 presents running times12 averaged over 10 random instances for each n.
There is a significant difference in running times between DSOS/SDSOS and SOS
beginning at n = 15. While we are unable to run SOS programs beyond n = 25,
DSOS and SDSOS programs for n = 70 take only a few minutes to execute. Memory
constraints prevent us from executing programs for n larger than 70. The results
indicate, however, that it may be possible to run larger programs within a reasonable
amount of time on a machine equipped with more RAM.
4.2. Copositive programming and combinatorial optimization. A symmetric matrix Q ∈ Rn×n is copositive if xT Qx ≥ 0 for all x ≥ 0 (i.e., for all x in the
nonnegative orthant). The problem of optimizing a linear function over affine sections
of the cone Cn of n × n copositive matrices has recently received a lot of attention
11 We

thank Aida Khajavirad for running the BARON experiments for us.
BARON instances were implemented on a 64-bit Intel Xeon X5650 2.66Ghz processor using
the CPLEX LP solver.
12 The

18

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

from the optimization community [27], since it can exactly model several problems of
combinatorial and nonconvex optimization [17], [27], [19]. For instance, the size α(G)
of the largest independent set 13 of a graph G on n nodes is equal to the optimal value
of the following copositive program [22]:

minimize
γ

(16)

s.t.

γ
γ(A + I) − J ∈ Cn ,

where A is the adjacency matrix of the graph, I is the identity matrix, and J is the
matrix of all ones.
It is easy to see that a matrix Q is copositive if and only if the quartic form
(x.2 )T Q(x.2 ), with x.2 := (x21 , x22 , . . . , x2n )T , is nonnegative. Hence, by requiring this
form to be r-sos/r-sdsos/r-dsos, one obtains inner approximations to the cone Cn that
can be optimized over using semidefinite, second order cone, and linear programming.
The sos version of this idea goes back to the PhD thesis of Parrilo [65].
In this section, we compare the quality of these bounds using another hierarchical
LP-based method in the literature [17], [22]. The r-th level of this LP hierarchy
(referred to as the Pólya LP from here onwards) requires that the coefficients of the
form (x.2 )T Q(x.2 )·(x21 +· · ·+x2n )r be nonnegative. This is motivated by a theorem by
Pólya, which states that this constraint will be satisfied for r large enough (assuming
xT Qx > 0 for all nonzero x ≥ 0). It is not difficult to see that the Pólya LP is always
dominated by the LP which requires (x.2 )T Q(x.2 ) to be r-dsos. This is because the
latter requires the underlying Gram matrix to be diagonally dominant, while the
former requires it to be diagonal, with nonnegative diagonal entries.
In [22], de Klerk and Pasechnik show that when the Pólya LP is applied to approximate the copositive program in (16), then the precise independent set number of
the graph is obtained within α(G)2 steps of the hierarchy. By the previous argument,
the same claim immediately holds for the r-dsos (and hence r-sdsos) hierarchies.
Table 3 revisits Example 5.2 of [17], where upper bounds on the independent set
number of the complement of the graph of an icosahedron are computed. (This is a
graph on 12 nodes with independent set number equal to 3; see [17] for its adjacency
matrix.) As the results illustrate, SOS programming provides an exact upper bound
(since the size of the largest stable set is necessarily an integer). The second levels
of the r-dsos and r-sdsos hierarchies also provide an exact upper bound. In contrast,
the LP based on Polya’s theorem in [17] gives the upper bound ∞ for the 0-th and
the 1-st level, and an upper bound of 6 at the second level. In contrast to the Pólya
LP, one can show that the 0-th level of the r-dsos LP always produces a finite upper
bound on the independent set number. In fact, this bound will always be smaller than
n − min di + 1, where di is the degree of node i [1, Theorem 5.1] .
4.3. Convex regression in statistics. Next, we consider the problem of fitting
a function to data subject to a constraint on the function’s convexity. This is an
important problem in statistics and has a wide domain of applications including value
function approximation in reinforcement learning, options pricing and modeling of
circuits for geometric programming based circuit design [34, 35]. Formally, we are
given N pairs of data points (xi , yi ) with xi ∈ Rn , yi ∈ R and our task is to find a
13 An

independent set of a graph is a subset of its nodes no two of which are connected. The
problem of finding upper bounds on α has many applications in scheduling and coding theory [52].

19

DSOS AND SDSOS OPTIMIZATION

r=0
r=1
r=2

Polya LP
∞
∞
6.000

r-DSOS
6.000
4.333
3.8049

r-SDSOS
6.000
4.333
3.6964

SOS
3.2362
NA
NA

Table 3
Comparison of upper bounds on the size of the largest independent set in the complement of
the graph of an icosahedron.

d=2
d=4

DSOS
35.11
14.86

SDSOS
33.92
12.94

SOS
21.28
NA

d=2
d=4

DSOS
0.49
145.85

SDSOS
0.66
240.18

SOS
0.53
∞

Table 4
Table 5
Comparison of fitting errors for convex reComparison of running times (in s) for congression.
vex regression.

convex function f from a family F that minimizes an appropriate notion of fitting
error (e.g., L1 error):
minimize
f ∈F

s.t.

N
X

|f (xi ) − yi |

i=1

f is convex.

The convexity constraint here can be imposed by requiring that the function wT H(x)w
in 2n variables (x, w) associated with the Hessian matrix H(x) of f (x) be nonnegative.
Restricting ourselves to polynomial functions f of bounded degree, we obtain the
following optimization problem:
minimize
f ∈Rd [x]

s.t.

N
X

|f (xi ) − yi |

i=1
T

w H(x)w ≥ 0, ∀x, w,

where H(x) is again the Hessian of f (x). As before, we can replace the nonnegativity
constraint with a dsos/sdsos/sos constraint. For our numerical experiment, we generated 300 random vectors xi in R20 drawn i.i.d. from the standard normal distribution.
The function values yi were computed as follows:
yi = exp(kxi k2 ) + η,
where η was chosen i.i.d. from the standard normal distribution. Tables 4 and 5
present the fitting errors and running times for the DSOS/SDSOS/SOS programs
resulting from restricting the class of functions F to polynomials of degree d = 2 and
d = 4. As the results illustrate, we are able to obtain significantly smaller errors with
polynomials of degree 4 using DSOS and SDSOS (compared to SOS with d = 2),
while the SOS program for d = 4 does not run due to memory constraints.
4.4. Options pricing. An important problem in financial economics is that of
determining the price of a derivative security such as an option given the price of the
underlying asset. Significant progress was made in tackling this problem with the
advent of the Black-Scholes formula, which operates under two assumptions: (i) the

20

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

underlying stock price is governed by geometric Brownian motion, and (ii) there is
no arbitrage. A natural question that has been considered in financial mathematics
is to see if one can provide bounds on the price of an option in the no-arbitrage
setting given much more minimal assumptions on the dynamics of the stock price.
One approach to this question—studied, e.g., in [13], [18], —is to assume that we are
only given the first k moments of the stock price and want to optimally bound the
price of the option.
More precisely, we consider that we are given m stocks, an associated option with
payoff function φ : Rm
+ → R (which will typically depend on the strike price of the
option), a vector of n moment functions, fi : R+
m → R, i = 0, 1, . . . , n (f0 is assumed
to be 1), and the corresponding vector of moments q = (q0 , . . . , qn ) (here q0 = 1). The
problem of upper bounding the price of the option can then be formulated as follows
[13]:
(17)

maximize

Eπ [φ(X)]

π

s.t.

Eπ [fi (X)] = qi , i = 0, . . . , n,
π(x) ≥ 0, ∀x ∈ Rm
+.

Here, the expectation is taken over all Martingale measures π, defined on Rm
+ . We
consider here the case where one is given the mean µ and covariance matrix σ of the
stock prices. Problem (17) can then be cast as the following problem (see [13, Section
6.2])
(18)

minimize
y,Y

y0 +

n
X
i=1

s.t.

yi µi +

n X
n
X

yij (σij + µi µj )

i=1 j=1

xT Y x + y T x + y0 ≥ φ(x), ∀x ≥ 0,

which for many common functions φ of interest (see e.g. below) gives rise to a copositive program. A well-known and obvious sufficient condition for a matrix M to be
copositive is for it to have a decomposition M = P + N , where P is psd and N is
element-wise nonnegative [65]. This allows one to apply SDP to problem (18) and obtain upper bounds on the option price. By replacing the psd condition on the matrix
P with a dd/sdd condition, we obtain a DDP/SDDP.
We first compare the DDP/SDDP/SDP approaches on an example from [18,
Section 4], which considers the problem of upper bounding the price of a European
call option with m = 3 underlying assets. The vector of means of the assets is
µ = [44.21, 44.21, 44.21]T and the covariance matrix is:


184.04 164.88 164.88
(19)
σ = 164.88 184.04 164.88 .
164.88 164.88 184.04
Further, we have the following payoff function, which depends on the strike price
K:
(20)

φ(x) = max(x1 − K, x2 − K, x3 − K, 0).

Table 6 compares the upper bounds for different strike prices using DDP/SDDP/SDP.
In each case, the upper bound obtained from the SDP is exact. This can be verified
by finding a distribution that achieves the upper bound (i.e., finding a matching

DSOS AND SDSOS OPTIMIZATION

K
K
K
K
K

= 30
= 35
= 40
= 45
= 50

Exact
21.51
17.17
13.20
9.84
7.30

SDP
21.51
17.17
13.20
9.84
7.30

SDDP
21.51
17.17
13.20
9.85
7.30

21

DDP
132.63
132.63
132.63
132.63
132.63

Table 6
Comparison of upper bounds on options prices for different strike prices obtained using SDP,
SDDP and DDP. Here, the number of underlying assets is m = 3.

Upper bound
Running time

SDP
18.76
2502.6 s

SDDP
19.42
24.36 s

DDP
252.24
11.85 s

Table 7
Comparison of upper bounds and running times for a large (m = 50) options pricing example
using SDP, SDDP and DDP.

lower bound). For example, when K = 30, the distribution supported on the set of
four points (5.971, 5.971, 5.971), (54.03, 46.02, 46.02, 46.02), (46.02, 54.03, 46.02) and
(46.02, 46.02, 54.03) in R3 with probability masses 0.105, 0.298, 0.298, 0.298 respectively does the job. The upper bound obtained using SDDP is almost identical to the
bound from SDP for each strike price. The DDP bound is loose, and interestingly
does not change with the strike price. Running times for the different methods on
this small example are negligible and hence not presented.
In order to demonstrate the scalability of DDP and SDDP, we consider a larger
scale randomized example with m = 50 underlying assets. The mean and covariance
are taken to be the sample mean and covariance of 100 instances of vectors belonging
to Rm with elements drawn uniformly and independently from the interval [0, 10].
The strike price is K = 5 and the payoff function is again
(21)

φ(x) = max(x1 − K, x2 − K, . . . , xm − K, 0).

Table 7 compares upper bounds and running times for the different methods.
SDDP provides us with a bound that is very close to the one obtained using SDP, with
a significant speedup in computation (approximately a factor of 100). The running
time for DDP is comparable to SDDP, but the bound is not as tight.
4.5. Sparse PCA. Next, we consider the problem of sparse principal component
analysis (sparse PCA). In contrast to standard PCA, where the principal components
(PCs) in general depend on all the observed variables, the goal of sparse PCA is to
identify principal components that only depend on small subsets of the variables [87].
While the statistical fidelity of the resulting representation of the data in terms of
sparse PCs will in general be lower than the standard PCs, sparse PCs can significantly
enhance interpretability of the results. This feature has proved to be very useful
in applications such as finance and analysis of gene expressions; see, e.g., [21] and
references therein.
Given a n × n covariance matrix A, the problem of finding sparse principal com-

22

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

ponents can be written as the following optimization problem:
(22)

maximize
x

s.t.

xT Ax
kxk2 = 1,
Card(x) ≤ k.

Here, Card(x) is the cardinality (number of non-zero entries) of x and k is a given
threshold. The problem can be reformulated as the following rank-constrained matrix
optimization problem [21]:
(23)

maximize

Tr(AX)

X

s.t.

Tr(X) = 1,
Card(X) ≤ k 2 ,
X  0, Rank(X) = 1.

In [21], the authors propose “DSPCA”, an SDP relaxation of this problem that is
obtained by dropping the rank constraint and replacing the cardinality constraint
with a constraint on the l1 norm:
(24)

maximize
X

s.t.

Tr(AX)
Tr(X) = 1,
1T |X|1 ≤ k,
X  0.

The optimal value of problem (24) is an upper bound on the optimal value of (22).
Further, when the solution X1 to (24) has rank equal to one, the SDP relaxation is
tight and the dominant eigenvector x1 of X1 is the optimal loading of the first sparse
PC [21]. When a rank one solution is not obtained, the dominant eigenvector can still
be retained as an approximate solution to the problem. Further sparse PCs can be
obtained by deflating A to obtain:
A2 = A − (xT1 Ax1 )x1 xT1
and re-solving the SDP with A2 in place of A (and iterating this procedure).
The framework presented in this paper can be used to obtain LP and SOCP
relaxations of (22) by replacing the constraint X  0 by the constraint X ∈ DDn∗
or X ∈ SDDn∗ respectively, where DDn∗ and SDDn∗ are the dual cones of DDn and
SDDn (see [3, Section 3.3] for a description of these dual cones). Since Sn+ ⊆ SDDn∗ ⊆
DDn∗ , we are guaranteed that the resulting optimal solutions will be upper bounds
on the SDP solution. Further, as in the SDP case, the relaxation is tight when a
rank-one solution is obtained, and once again, if a rank-one solution is not obtained,
we can still use the dominant eigenvector as an approximate solution.
We first consider Example 6.1 from [21]. In this example, there are three hidden
variables distributed normally:
V1 ∼ N (0, 290), V2 ∼ N (0, 300), V3 = −0.3V1 + 0.925V2 + ,  ∼ N (0, 1).
These hidden variables generate 10 observed variables:
Xi = Vj + ji , ji ∼ N (0, 1),

DSOS AND SDSOS OPTIMIZATION

23

X1
X2
X3
X4
X5
X6
X7
X8
X9
X10 expl. var.
PCA, PC1 0.116 0.116 0.116 0.116 -0.395 -0.395 -0.395 -0.395 -0.401 -0.401 60.0 %
PCA, PC2 -0.478 -0.478 -0.478 -0.478 -0.145 -0.145 -0.145 -0.145 0.010 0.010
39.6 %
Other, PC1
0
0
0
0
-0.5
-0.5
-0.5
-0.5
0
0
40.9 %
Other, PC2 0.5
0.5
0.5
0.5
0
0
0
0
0
0
39.5 %
Table 8
Comparison of loadings of principal components and explained variances for standard PCA and
sparse versions. Here, the label “Other” denotes sparse PCA based on SDP, DDP, and SDDP. Each
of these gives the optimal sparse solution.

DDP
SDDP
DSPCA
No. Time (s) NNZ Opt. Expl. Time (s) NNZ Opt. Expl. Time (s) NNZ Opt. Expl.
1
0.89
4 35.5 0.067 %
1.40
14 30.6 0.087 % 1296.9
5 30.5 0.086 %
2
1.18
4 36.6 0.071 %
1.42
14 34.4 0.099 % 1847.3
6 33.6 0.092 %
3
1.46
4 49.5 0.079 %
1.40
24 41.3 0.097 % 1633.0
16 40.2 0.096 %
4
1.14
4 44.1 0.072 %
1.80
17 38.7 0.10 %
1984.7
8 37.6 0.091 %
5
1.10
4 36.7 0.060 %
1.53
34 33.0 0.068 % 2179.6
10 31.7 0.105 %
Table 9
Comparison of running times, number of non-zero elements, optimal values, and explained
variances for five large-scale sparse PCA examples using covariance matrices of size 100 × 100.

with j = 1 for i = 1, . . . , 4, j = 2 for i = 5, . . . , 8 and j = 3 for i = 9, 10, and ji
independent for j = 1, 2, 3 and i = 1, . . . , 10. This knowledge of the distributions of the
hidden and observed variables allows us to compute the exact 10×10 covariance matrix
for the observed variables. The sparse PCA algorithm described above can then be
applied to this covariance matrix. Table 8 presents the first two principal components
computed using standard PCA, DDP, SDDP and DSPCA (corresponding to SDP).
As the table illustrates, the loadings corresponding to the first two PCs computed
using standard PCA are not sparse. All other methods (DDP, SDDP, and SDP) give
exactly the same answer and correspond to a rank-one solution, i.e., the optimal sparse
solution. As expected, the sparsity comes at the cost of a reduction in the variance
explained by the PCs (see [87] for computation of the explained variance).
Next, we consider a larger-scale example. We generate five random 100 × 100 covariance matrices of rank 4. Table 9 presents the running times, number of non-zero
entries (NNZ) in the top principal component14 , optimal value of the program (24)
(Opt.), and the explained variance (Expl.). The optimal values for DDP and SDDP
upper bound the SDP solution as expected, and are quite close in value. The number of non-zero entries and explained variances are comparable across the different
methods. The running times are between 1100 and 2000 times faster for DDP in comparison to SDP, and between 900 and 1400 times faster for SDDP. Hence the example
illustrates that one can obtain a very large speedup with the approach presented here,
with only a small sacrifice in quality of the approximate sparse PCs.
4.6. Applications in control theory. As a final application of the methods
presented in this paper, we consider two examples from control theory and robotics.
The applications of SOS programming in control theory are numerous and include
the computation of regions of attraction of polynomial systems [65], feedback control
synthesis [40], robustness analysis [83], and computation of the joint spectral radius
for uncertain linear systems [64]. A thorough treatment of the application of the
14 Note

that the entries of the PCs were thresholded lightly in order to remove spurious non-zero
elements arising from numerical inaccuracies in the solvers.

24

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

0.25

0.4
SOS
SDSOS
DSOS

0.2

SOS
SDSOS
DSOS

0.3

0.15
0.2
0.1
0.1
θ̇6 (rad/s)

θ̇1 (rad/s)

0.05

0

−0.05

0

−0.1

−0.1
−0.2
−0.15
−0.3

−0.2

−0.25
−0.05

Fig. 4. An illustration of the Nlink pendulum system (with N = 6).

−0.04

−0.03

−0.02

−0.01

0
0.01
θ1 (radians)

0.02

0.03

(a) θ1 -θ̇1 subspace.

0.04

0.05

−0.4
−0.05

−0.04

−0.03

−0.02

−0.01

0
0.01
θ6 (radians)

0.02

0.03

0.04

0.05

(b) θ6 -θ̇6 subspace.

Fig. 5. Figure reproduced from [53] comparing projections of the ROAs
computed for the 6-link pendulum system using DSOS, SDSOS and SOS
programming.

(S)DSOS approach to control problems is beyond the scope of this paper but can be
found in a different paper [53], which is joint work with Tedrake. Here we briefly
highlight two examples from [53] that we consider particularly representative.
Most of the applications in control theory mentioned above rely on SOS programming for checking Lyapunov inequalities that certify stability of a nonlinear system.
As discussed in Section 1.1, one can compute (inner approximations of) the region of
attraction (ROA) of a dynamical system ẋ = f (x) by finding a Lyapunov function
V : Rn → R that satisfies the following conditions:

(25)

V (x) > 0

∀x 6= 0,

and

V̇ (x) = h∇V (x), f (x)i < 0

∀x ∈ {x| V (x) ≤ β, x 6= 0}.

This guarantees that the set {x| V (x) ≤ β} is a subset of the ROA of the system,
i.e. initial conditions that begin in this sublevel set converge to the origin. For polynomial dynamical systems, one can specify sufficient conditions for (25) using SOS
programming [65]. By replacing the SOS constraints with (S)DSOS constraints, we
can compute inner approximations to the ROA more efficiently. While the approximations will be more conservative in general, the ability to tradeoff conservatism with
computation time and scalability is an important one in control applications.
4.6.1. Region of attraction for an inverted N -link pendulum. As an
illustration, we consider the problem of computing ROAs for the underactuated N -link
pendulum depicted in Figure 4. This system has 2N states x = [θ1 , . . . , θN , θ̇1 , . . . , θ̇N ]
composed of the joint angles (with the vertical line) and their derivatives. There
are N − 1 control inputs (the joint closest to the base is not actuated). We take
the unstable “upright” position of the system to be the origin of our state space
and design a Linear Quadratic Regulator (LQR) controller in order to stabilize this
equilibrium. A polynomial approximation of the dynamics of the closed loop system
is obtained by a Taylor expansion of order 3. Using Lyapunov functions of degree 2
results in the time derivative of the Lyapunov function being quartic and hence yields
dsos/sdsos/sos constraints on a polynomial of degree 4 in 2N variables.
Figures 5(a) and 5(b) compare projections of the ROAs computed for the system
with N = 6 onto two 2−dimensional subspaces of the state space. As the plots
suggest, the ROA computed using SDSOS is only slightly more conservative than the
ROA computed using SOS programming. In fact, comparing the volumes of the two
ROAs, we find:

DSOS AND SDSOS OPTIMIZATION

25

Fig. 6. A visualization of the model of the ATLAS humanoid robot, along with the hardware
platform (inset) on which the parameters of the model are based. (Picture of robot reproduced with
permission from Boston Dynamics.)

1/2N

V olROA-sdsos

(26)

1/2N

= 0.79.

V olROA-sos
Here, the volumes have been correctly rescaled in the standard manner by the
dimension of the ambient space. The ROA computed using DSOS is much smaller
(but still potentially useful in practice).
Comparing running times of the different approaches on this example, we find
that the LP corresponding to the DSOS program takes 9.67 seconds, and the SOCP
corresponding to the SDSOS program takes 25.9 seconds. The running time of the
SOS program is 1526.5 seconds using MOSEK and 23676.5 seconds using SeDuMi.
Hence, in particular DSOS is approximately 2500 times faster than SOS using SeDuMi
and 150 times faster than SOS using MOSEK, while SDSOS is 900 times faster in
comparison to SeDuMi and 60 times faster than MOSEK for SOS.
Of course the real benefit of our approach is that we can scale to problems where
SOS programming ceases to run due to memory/computation constraints. Table 10
illustrates this ability by comparing running times of the programs obtained using our
approach with SOS programming for different values of 2N (number of states). As
expected, for cases where the SOS programs do run, the DSOS and SDSOS programs
are significantly faster. Further, the SOS programs obtained for 2N > 12 are too large
to run (due to memory constraints). In contrast, our approach allows us to handle
almost twice as many states.
2N (# states)
DSOS
SDSOS
SOS (SeDuMi)
SOS (MOSEK)

4
<1
<1
<1
<1

6
0.44
0.72
3.97
0.84

8
2.04
6.72
156.9
16.2

10
12
14
16
18
20
22
3.08
9.67 25.1 74.2 200.5 492.0 823.2
7.78
25.9 92.4 189.0 424.74 846.9 1275.6
1697.5 23676.5 ∞
∞
∞
∞
∞
149.1 1526.5 ∞
∞
∞
∞
∞

Table 10
Table reproduced from [53] showing runtime comparisons (in seconds) for ROA computations
on N-link system.

4.6.2. Control synthesis for a humanoid robot. In our final example, we
highlight a robotics application considered in joint work with Tedrake in [53]. Control
of humanoid robots is an important problem in robotics and presents a significant
challenge due to the nonlinear dynamics of the system and high dimensionality of the
state space. Here we show how the approach described in this paper can be used to

26

(a) Nominal
pose
(fixed
point)

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

(b) Stabilized
pose 1

(c) Stabilized
pose 2

(d) Stabilized
pose 3

(e) Stabilized
pose 4

(f) Stabilized
pose 5

Fig. 7. Figure reproduced from [53] showing the nominal position of the robot, i.e., the fixed
point being stabilized (subplot (a)), and configurations of the robot that are stabilized by the controller
designed using SDSOS programing (subplots (b)-(f )). A video of simulations of the controller started
from different initial conditions is available online at http://youtu.be/lmAT556Ar5c.

design a balancing controller for a model of the ATLAS robot shown in Figure 6. This
robot was designed and built by Boston Dynamics Inc. and was used for the 2015
DARPA Robotics Challenge.
Our model of the robot is based on physical parameters of the hardware platform
and has 30 states and 14 inputs. The task considered here is to balance the robot on
its right toe. The balancing controller is constructed by searching for both a quadratic
Lyapunov function and a linear feedback control law in order to maximize the size of
the resulting region of attraction. The Lyapunov conditions are imposed using SDSOS
programming. We Taylor expand the dynamics about the equilibrium to degree 3 in
order to obtain polynomial dynamics. The total computation time is approximately
22.5 minutes. We note that SOS programming is unable to handle this system due to
memory (RAM) constraints.
Figure 7 demonstrates the performance of the resulting controller from SDSOS
programming by plotting initial configurations of the robot that are stabilized to the
fixed point. As the plot illustrates, the controller is able to stabilize a very wide range
of initial conditions. A video of simulations of the closed loop system started from
different initial conditions is available online at http://youtu.be/lmAT556Ar5c.

5. Improvements on DSOS and SDSOS optimization. While DSOS and
SDSOS techniques result in significant gains in terms of solving times and scalability,
they inevitably lead to some loss in solution accuracy when compared to the SOS
approach. In this section, we briefly outline three possible strategies to mitigate this
loss. The first two (Sections 5.1 and 5.2) generate sequences of linear and second order
cone programs, while the third (Section 5.3) works with “small” (fixed-size) semidefinite programs. The reader may recall that the r-DSOS and r-SDSOS hierarchies of
Section 3.2 could also be used to improve on DSOS and SDSOS techniques. Like the
first two methods that we present here, they did so while staying in the realm of LP
and SOCP. They differ from these methods on two crucial points however: (i) they are
more expensive to implement, and (ii) they approximate the cone of SOS polynomials
“blindly” irrespective of a particular objective function. This is in contrast to the
methods that we present in Section 5.1 and 5.2, which approximate the cone in the
direction of a specific linear objective function.
For brevity of exposition, we explain how all three strategies can be applied to

DSOS AND SDSOS OPTIMIZATION

27

approximate a generic semidefinite program:
(27)

minimize

Tr(CX)

X∈Sn

s.t.

Tr(Ai X) = bi ,

i = 1, . . . , m,

X  0.
A treatment tailored to the case of sum of squares programs can be found in the
references we provide.
5.1. Iterative change of basis. In [3], Ahmadi and Hall build on the notions
of diagonal and scaled diagonal dominance to provide a sequence of improving inner
approximations to the cone Pn of psd matrices in the direction of the objective function
of an SDP at hand. The idea is simple: define a family of cones15
DD(U ) := {M ∈ Sn | M = U T QU for some dd matrix Q},
parametrized by an n × n matrix U . Optimizing over the set DD(U ) is an LP
since U is fixed, and the defining constraints are linear in the coefficients of the two
unknown matrices M and Q. Furthermore, the matrices in DD(U ) are all psd; i.e.,
∀U, DD(U ) ⊆ Pn .
The proposal in [3] is to solve a sequence of LPs, indexed by k, by replacing the
condition X  0 by X ∈ DD(Uk ):
(28)

DSOSk := min

X∈Sn

s.t.

Tr(CX)
Tr(Ai X) = bi ,

i = 1, . . . , m,

X ∈ DD(Uk ).
The sequence of matrices {Uk } is defined as follows
(29)

U0 = I
Uk+1 = chol(Xk ),

where Xk is an optimal solution to the LP in (28).
Note that the first LP in the sequence optimizes over the set of diagonally dominant matrices. By defining Uk+1 as a Cholesky factor of Xk , improvement of the
T
IUk+1 , and the
optimal value is guaranteed in each iteration. Indeed, as Xk = Uk+1
identity matrix I is diagonally dominant, we see that Xk ∈ DD(Uk+1 ) and hence is
feasible for iteration k + 1. This implies that the optimal value at iteration k + 1 is at
least as good as the optimal value at the previous iteration; i.e., DSOSk+1 ≤ DSOSk
(in fact, the inequality is strict under mild assumptions; see [3]).
In an analogous fashion, one can construct a sequence of SOCPs that inner approximate Pn with increasing quality. This time, we define a family of cones
SDD(U ) := {M ∈ Sn | M = U T QU, for some sdd matrix Q},
parameterized again by an n × n matrix U . For any U , optimizing over the set
SDD(U ) is an SOCP and we have SDD(U ) ⊆ Pn . This leads us to the following
iterative SOCP sequence:
15 One

can think of DD(U ) as the set of matrices that are dd after a change of coordinates via
the matrix U .

28

(30)

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

Tr(CX)

SDSOSk := min

X∈Sn

s.t.

Tr(Ai X) = bi ,

i = 1, . . . , m,

X ∈ SDD(Uk ).
Assuming existence of an optimal solution Xk at each iteration, we can define
the sequence {Uk } iteratively in the same way as was done in (29). Using similar
reasoning, we have SDSOSk+1 ≤ SDSOSk . In practice, the sequence of upper
bounds {SDSOSk } approaches faster to the SDP optimal value than the sequence of
the LP upper bounds {DSOSk }. Figure 8 shows the improvement (in every direction)

(a) LP inner approximations

(b) SOCP inner approximations

Fig. 8. Figure reproduced from [3] showing improvement (in all directions) after one iteration
of the change of basis algorithm.

obtained just by a single iteration of this approach. The outer set in green in both
subfigures is the feasible set of a randomly generated semidefinite program. The sets in
black are the DD (left) and the SDD (right) inner approximations. What is shown in
dark blue in both cases is the boundary of the improved inner approximation after one
iteration. Note that the SOCP in particular fills up almost the entire spectrahedron
in a single iteration.
We refer the interested reader to [3] for more details, in particular for an explanation of how the same techniques can be used (via duality) to outer approximate
feasible sets of semidefinite programs.
5.2. Column generation. In [1], Ahmadi, Dash, and Hall design another iterative method for inner approximating the set of psd matrices using linear and second
order cone programming. Their approach combines DSOS/SDSOS techniques with
ideas from the theory of column generation in large-scale linear and integer programming. The high-level idea is to approximate the SDP in (27) by a sequence of LPs
(parameterized by t):
(31)

minimize
X∈Sn ,αi

s.t.

Tr(CX)
Tr(Ai X) = bi ,
X=

t
X

i = 1, . . . , m,

αi Bi ,

i=1

(32)

αi ≥ 0,

i = 1, . . . , t,

DSOS AND SDSOS OPTIMIZATION

29

where B1 , . . . , Bt are fixed psd matrices. These matrices are initialized to be the
extreme rays of DDn (recall Lemma 6), i.e., all rank one matrices vi viT , where the
vector vi has at most two nonzero components, each equal to ±1. Once this initial LP
is solved, then one adds one (or sometimes several) new psd matrices Bj to problem
(31) and resolves. This process then continues. In each step, the new matrices Bj are
picked carefully to bring the optimal value of the LP closer to that of the SDP. Usually,
the construction of Bj involves solving a “pricing subproblem” (in terminology of the
column generation literature), which adds appropriate cutting planes to the dual of
(31); see [1] for more details.
The SOCP analogue of this process is similar. The SDP in (27) is inner approximated by a sequence of SOCPs (parameterized by t):
(33)

minimize

X∈Sn ,Λi ∈S2

s.t.

Tr(CX)
Tr(Ai X) = bi ,
X=

t
X

i = 1, . . . , m,

Vi Λi ViT ,

i=1

Λi  0,

i = 1, . . . , t,

where V1 , . . . , Vt are fixed n × 2 matrices. They are initialized as the set of matrices
that have zeros everywhere, except for a 1 in the first column in position j and a 1
in the second column in position k 6= j. This gives exactly SDDn (via Lemma 9).
In subsequent steps, one (or sometimes several) appropriately-chosen matrices Vi are
added to problem (33). These matrices are obtained by solving pricing subproblems
and help bring the optimal value of the SOCP closer to that of the SDP in each
iteration.

(a) LP iterations.

(b) SOCP iterations.

Fig. 9. Figure reproduced from [1] showing the successive improvement on the dd (left) and sdd
(right) inner approximation of a spectrahedron via five iterations of the column generation method.

Figure 9 shows the improvement obtained by five iterations of this process on a
randomly generated SDP, where the objective is to maximize a linear function in the
northeast direction. 5
5.3. Factor-width k matrices with k > 2. A relevant generalization of sdd
matrices is through the notion of factor-width, defined by Boman et al. [16]. The
factor-width of a symmetric psd matrix Q is the smallest integer k for which Q can

30

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

be written as Q = V V T , where each column of V contains at most k nonzeros.
Equivalently, the factor-width of Q is the smallest k for which Q can be written as a
sum of psd matrices that are nonzero only on a single k × k principal submatrix.
If we denote the cone of n × n symmetric matrices of factor-width k by F Wnk ,
we have that F Wnk ⊆ Pn for all k ∈ {1, . . . , n}, F Wn2 = SDDn (cf. Lemma 9), and
F Wnn = Pn . For values of k ∈ {2, . . . , n − 1}, one gets an increasingly more accurate
inner approximation to the set of psd matrices, all outperforming the approximation
given by sdd matrices. We did not consider these cones in this paper as it is already
known that for k = 3, a representation based on second order cone programming is not
possible [28]. Nevertheless, working with F Wn3 can be useful as it involves semidefinite
constraints that are small (3×3) and hence efficiently solvable. Unfortunately though,
optimizing over F Wn3 requires O(n3 ) such semidefinite constraints which in many
applications can be prohibitive. We believe that a promising future direction would
be to use the column generation framework on Section 5.2 to work with a subset of
the extreme rays of F Wnk (for small k) and generate additional ones on the fly based
on problem structure.
Some initial experiments with F Wnk have been performed by Ding and Lim [24].
The authors also develop self-concordant barrier functions for these cones and show
that any second order cone program can be written as an optimization problem over
F Wn2 , i.e., a scaled diagonally dominant program (recall the definition from Section 3.1). In related work, Permenter and Parrilo [68] consider optimization problems
over F Wnk (and their dual cones) for facial reduction in semidefinite programming.
We believe the notion of factor-width is likely to receive more attention in upcoming
years in the theory of matrix optimization.
6. Conclusions. We have proposed more scalable alternatives to SOS programming by introducing inner approximations to the sum of squares cone in the form of
the cones of dsos and sdsos polynomials. These cones can be optimized over using
LP and SOCP respectively. The departure from SDP allows us to obtain significant
gains in terms of scalability. Our numerical examples from a diverse range of applications including polynomial optimization, combinatorial optimization, statistics and
machine learning, derivative pricing, control theory, and robotics demonstrate that
with reasonable tradeoffs in optimality, we can handle problem sizes that are well
beyond the current capabilities of SOS programming. In particular, we have shown
that our approach is able to tackle dense problems with as many as 70 polynomial
variables (with degree-4 polynomials). On the theoretical front, we have shown that
the (S)DSOS approach shares many of the theoretical asymptotic guarantees usually associated with SOS programming. Finally, in the last section of this paper, we
reviewed recent approaches that were developed with the idea of bridging the gap
between the (S)DSOS approach and the SOS one. One can obtain improved optimal
values using these methods at the cost of additional computation time. The appendix
provides a brief but complete tutorial on the toolbox that was used to generate the
numerical results in the paper. The accompanying code is freely available online.
We would like to emphasize that a key benefit of our approach is that it can
be applied in any application where SOS programming is used in order to obtain
potentially significant gains in scalability. Our hope is that the (S)DSOS approach
may open up application areas that have previously been beyond reach due to the
limitations in scalability of SOS programming. We also foresee possibilities for realtime applications [4], [55]. Indeed, technology for real-time linear and second order
cone programming is already coming to fruition [26], [56].

DSOS AND SDSOS OPTIMIZATION

31

7. Appendix. A complete implementation of the code used to generate the numerical results presented in this paper was written using the Systems Polynomial Optimization Toolbox (SPOT) [57] and is freely available online16 . The toolbox features
efficient polynomial algebra and allows us to setup the large-scale LPs and SOCPs
arising from our examples. In this appendix, we provide a brief introduction to the
software. This is not meant as a comprehensive tutorial on the SPOT toolbox (for
this one may refer to SPOT’s documentation17 ). Rather, the goal here is to provide
an introduction sufficient for setting up and solving the DSOS and SDSOS programs
in this paper. The code relevant for this purpose is included in a branch of SPOT
that we have named iSOS (for “inside sum of squares”).
7.1. Installing SPOT. Download the software package from Github available
here:
https://github.com/anirudhamajumdar/spotless/tree/spotless isos
Next, start MATLAB and run the spot install.m script. This script will setup the
MATLAB path and compile a few mex functions. The user may wish to save the new
MATLAB path for future use.
7.2. Variables and polynomials. Polynomials are defined and manipulated
by the @msspoly class of SPOT. In order to define a new variable, one can use the
msspoly.m function:
>> x = msspoly(‘x’)
which creates the polynomial p(x) = x. The argument to this function is the name
of the created variable and is restricted to four characters chosen from the alphabet
(lower case and upper case). A MATLAB vector of variables can be created by passing
a second argument to the function:
>> x = msspoly(‘x’,n).
This will create a n × 1 vector of variables that can be accessed using standard
array indexing. Multivariate polynomials can then be constructed as follows:

>> x = msspoly(‘x’,3);
>> p = 2*x(1)^2 - 5*x(2)^2 + x(3)^2
Variables can be manipulated and operated on using a variety of functions. These
include standard arithmetic operations (e.g. addition, multiplication, dot product)
and operations for manipulating vectors (e.g. concatenating, reshaping, transposing).
Other useful functions include:
deg.m: Returns the total degree of a polynomial. If a second argument in the
form of a msspoly is provided, the degree with respect to these variables is returned.
16 Link

to the software:
https://github.com/anirudhamajumdar/spotless/tree/spotless isos
17 Link to SPOT’s documentation:
https://github.com/spot-toolbox/spotless/blob/master/doc/manual.pdf

32

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

diff.m: Differentiates a polynomial (first argument) with respect to a set of
msspoly variables (second argument). The result is the matrix of partial derivatives.
subs.m: Substitutes the third argument in place of the second argument wherever
it appears in the first argument.
7.3. Programs. Programs involving DSOS/SDSOS/SOS constraints are constructed using the spotsosprog class. In this section, we demonstrate the workings
of spotsosprog with the help of an example. In particular, we take the problem of
minimizing a form on the unit sphere considered in this paper. The following example
is also available in doc/examples/sdsos_example.m.

% Construct polynomial which is to be minimized
x = msspoly(‘x’,6);
vx = monomials(x,4:4);
randn(‘state’,0)
cp = randn(1,length(vx));
p = cp*vx;
This block of code constructs the polynomial that is to be lower bounded. In order to
do this, we create a six dimensional vector x of variables using the msspoly command
introduced before. Next, we construct a vector of monomials using the monomials
function. The first input to this function is the msspoly variable over which the
monomials are defined. The second input to the function is the range of degrees the
monomials should have. In this case, since we are considering homogeneous quartics,
the range is simply 4 : 4 (i.e., just 4).

% Build program
prog = spotsosprog;
prog = prog.withIndeterminate(x);
[prog,gamma] = prog.newFree(1);
prog = prog.withDSOS(p - gamma*(x’*x)^2);
This block of code sets up the program and constraints. First, the program is
initialized in the form of the variable prog. This object will contain information
about constraints and decision variables. Next, we declare the variable x to be an
indeterminate or abstract variable (thus distinguishing it from a decision variable).
Decision variables are created using the newFree function in the class spotsosprog.
The input to this function is the number of new decision variables to be created.
The outputs are the updated program and a variable corresponding to the decision
variable. In our case, we need a single variable gamma to be declared. Finally, we
specify a DSOS constraint on p - gamma*(x’*x)^2 using the withDSOS command.
There are corresponding functions withSDSOS and withSOS that can be used to setup

DSOS AND SDSOS OPTIMIZATION

33

SDSOS or SOS constraints.
% Setup options and solve program
options = spot sdp default options();
% Use just the interior point algorithm to clean up
options.solveroptions.MSK IPAR BI CLEAN OPTIMIZER = ...
...

‘MSK OPTIMIZER INTPNT’;

% Don’t use basis identification
options.solveroptions.MSK IPAR INTPNT BASIS = ‘MSK BI NEVER’;
% Display solver output
options.verbose = 1;
sol = prog.minimize(-gamma, @spot mosek, options);
%Get value of gamma for optimal solution
gamma optimal = double(sol.eval(gamma))
Finally, in this block of code, we dictate an options structure for the program. In
particular, the field options.solveroptions contains options that are specific to the
solver to be used (MOSEK in our case). The minimize command is used to specify the
objective of the program (which must be linear in the decision variables), the function
handle corresponding to the solver to be used, and the structure of options. Currently,
MOSEK, Gurobi and SeDuMi are the supported solvers. So, for example, in order to
use the Gurobi solver for a (S)DSOS program, we would specify the function handle
@spot gurobi.
The output of the minimize command is a solution structure, which contains
diagnostic information and allows one to access the optimized decision variables. The
last line of our example code demonstrates how to obtain the optimized variable gamma
and convert it to a MATLAB double type.
7.4. Additional functionality. Next, we review additional functionality not
covered in the example above. While this section is not meant to be an exhaustive
list of all the functionality available in SPOT, it should be sufficient to reproduce
the numerical results presented in this paper. In the lines of code presented in the
following sections, it is assumed that prog is an object of the spotsosprog class, x is
a msspoly variable of size n and an indeterminate variable of prog. These variables
can be initialized as in the example above with the following lines of code:
prog = spotsosprog;
x = msspoly(‘x’,n);
prog = prog.withIndeterminate(x);
7.4.1. Creating decision variables. It is often useful to construct a polynomial whose coefficients are decision variables. This can be achieved with the newFree
and monomials functions introduced above. In particular, one can create a vector
of monomials in a msspoly variable x of a given degree d using the command v =
monomials(x,0:d). Then, the set of coefficients can be declared using [prog,c] =

34

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

prog.newFree(length(v)). Finally, one can obtain the desired polynomial by multiplying these two together: p = c’*v.
It is also possible to create other types of decision variables. For example, the
following functions can be used to create various types of matrix decision variables:
newSym: New symmetric matrix.
newDD: New symmetric matrix constrained to be diagonally dominant.
newSDD: New symmetric matrix constrained to be scaled diagonally dominant.
newPSD: New symmetric positive semidefinite matrix.
newDDdual: New symmetric matrix constrained to lie in the dual of the cone of
diagonally dominant matrices.
newSDDdual: New symmetric matrix constrained to lie in the dual of the cone of
scaled diagonally dominant matrices.
The input to these functions is the size of the desired (square) matrix. For
example, in order to create a 10 × 10 symmetric matrix Q constrained to be diagonally
dominant, one can use the following lines of code:
[prog,Q] = prog.newDD(10);
7.4.2. Specifying constraints. In addition to the withDSOS, withSDSOS and
withSOS functions introduced previously, constraints on existing decision variables
can be specified using the following functions:
withEqs: Sets up constraints of the form expr = 0, where expr is the input to
the function and is a matrix whose elements are to be constrained to be equal to 0.
withPos: Sets up constraints of the form expr ≥ 0, where expr is the input to
the function and is a matrix whose elements are to be constrained to be nonnegative.
Note that the functions withEqs and withPos allow one to specify constraints in
a vectorized manner, thus avoiding MATLAB’s potentially-slow for-loops.
withDD: Constrains a symmetric matrix (the input to the function) to be diagonally dominant.
withSDD: Constrains a symmetric matrix (the input to the function) to be scaled
diagonally dominant.
withPSD: Constrains a symmetric matrix (the input to the function) to be positive
semidefinite.
In each case above, the inputs must be affine functions of the decision variables of
the program. The output of these functions is an updated spotsosprog object that
contains the new constraints. We illustrate the use of the function withEqs with the
help of a simple example. The other functions can be used in a similar manner.
prog = spotsosprog;
[prog,tau1] = prog.newFree(10);
[prog,tau2] = prog.newFree(10);
(34)

prog = prog.withEqs(tau1 - tau2);

In this example, the elements of the 10 × 1 decision vector tau1 are constrained
to be equal to the elements of the vector tau2.
withPolyEqs: In order to constrain two polynomials p1 (x) and p2 (x) to be equal
to each other for all x, one may use the function withPolyEqs. This function will

DSOS AND SDSOS OPTIMIZATION

35

constrain the coefficients of the two polynomials to be equal. Here is a simple example:
prog = spotsosprog;
x = msspoly(‘x’,2);
prog = prog.withIndeterminate(x);
[prog,c] = prog.newFree(4);
p1 = (c(1) + c(2))*x(1)^2 + c(3)*x(2)^2;
p2 = 2*x(1)^2 + c(4)*x(2)^2;
prog = prog.withPolyEqs(p1 - p2);
This will constrain c(1) + c(2) = 2 and c(3) = c(4).
7.4.3. Checking if a polynomial is dsos. The three utility functions isDSOS,
isSDSOS and isSOS allow one to check if a given polynomial is dsos, sdsos, or sos
respectively. The only input to the functions is the polynomial to be checked. The
outputs are a boolean variable indicating whether the polynomial is in fact dsos/ sdsos/ sos, the Gram matrix, and the monomial basis corresponding to the Gram matrix
(these last two are non-empty only if the given polynomial is in fact dsos/sdsos/sos).
The following lines of code provide a simple example:
x = msspoly(‘x’,3);
p = x(1)^2 + 5*x(2)^2 + 3*x(3)^2;
(35)

[isdsos,Q,v] = isDSOS(p)

One can check that the polynomial p - v’*Q*v is the zero polynomial (up to
numerical tolerances).
Acknowledgments. We would like to thank Pablo Parrilo for acquainting us
with scaled diagonally dominant matrices, and Georgina Hall for many corrections and
simplifications on the first draft of this work. Our gratitude extends to Russ Tedrake
and the members of the Robot Locomotion Group at MIT for several insightful discussions, particularly around control applications. We thank Frank Permenter and
Mark Tobenkin for their help with the numerical implementation of DSOS/SDSOS
programs in SPOT, and Aida Khajavirad for running the experiments of Section 4.1
with BARON. Finally, this work has benefited from questions and comments from
several colleagues, among whom we gladly acknowledge Greg Blekherman, Stephen
Boyd, Sanjeeb Dash, Etienne de Klerk, Jesus de Loera, Lijun Ding, Didier Henrion,
Jean Bernard Lasserre, Lek-Heng Lim, Bruce Reznick, James Saunderson, and Bernd
Sturmfels.
REFERENCES
[1] A. A. Ahmadi, S. Dash, and G. Hall, Optimization over structured subsets of positive
semidefinite matrices via column generation, Discrete Optimization, (2016).
[2] A. A. Ahmadi and G. Hall, On the construction of converging hierarchies for polynomial
optimization based on certificates of global positivity. Available at https://arxiv.org/abs/
1709.09307, 2017.
[3] A. A. Ahmadi and G. Hall, Sum of squares basis pursuit with linear and second order cone
programming, Contemporary Mathematics, (2017).
[4] A. A. Ahmadi and A. Majumdar, Some applications of polynomial optimization in operations
research and real-time decision making, Optimization Letters, 10 (2016), pp. 709–729.

36

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

[5] F. Alizadeh and D. Goldfarb, Second-order cone programming, Mathematical programming,
95 (2003), pp. 3–51.
[6] E. Artin, Über die Zerlegung Definiter Funktionen in Quadrate, Hamb. Abh., 5 (1927),
pp. 100–115.
[7] C. Bachoc and F. Vallentin, New upper bounds for kissing numbers from semidefinite programming, Journal of the American Mathematical Society, 21 (2008), pp. 909–924.
[8] G. Barker and D. Carlson, Cones of diagonally dominant matrices, Pacific Journal of Mathematics, 57 (1975), pp. 15–32.
[9] A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications, vol. 2, Siam, 2001.
[10] D. P. Bertsekas and J. N. Tsitsiklis, Introduction to Probability, vol. 1, Athena Scientific
Belmont, MA, 2002.
[11] D. Bertsimas, R. M. Freund, and X. A. Sun, An accelerated first-order method for solving
SOS relaxations of unconstrained polynomial optimization problems, Optimization Methods and Software, 28 (2013), pp. 424–441.
[12] D. Bertsimas, D. A. Iancu, and P. A. Parrilo, A hierarchy of near-optimal policies for
multistage adaptive optimization, IEEE Transactions on Automatic Control, 56 (2011),
pp. 2809–2824.
[13] D. Bertsimas and I. Popescu, On the relation between option and stock prices: a convex
optimization approach, Operations Research, 50 (2002), pp. 358–374.
[14] D. Bertsimas and I. Popescu, Optimal inequalities in probability theory: a convex optimization approach, SIAM Journal on Optimization, 15 (2005), pp. 780–804.
[15] G. Blekherman, P. A. Parrilo, and R. Thomas, Semidefinite Optimization and Convex
Algebraic Geometry, SIAM Series on Optimization, 2013.
[16] E. G. Boman, D. Chen, O. Parekh, and S. Toledo, On factor width and symmetric Hmatrices, Linear Algebra and its Applications, 405 (2005), pp. 239–248.
[17] I. M. Bomze and E. De Klerk, Solving standard quadratic optimization problems via linear, semidefinite and copositive programming, Journal of Global Optimization, 24 (2002),
pp. 163–185.
[18] P. P. Boyle and X. S. Lin, Bounds on contingent claims based on several assets, Journal of
Financial Economics, 46 (1997), pp. 383–400.
[19] S. Burer, Copositive Programming, in Handbook on Semidefinite, Conic and Polynomial Optimization, Springer, 2012, pp. 201–218.
[20] CPLEX, V12. 2: Users manual for CPLEX, International Business Machines Corporation, 46
(2010), p. 157.
[21] A. d Aspremont, L. E. Ghaoui, M. I. Jordan, and G. R. Lanckriet, A direct formulation
for sparse pca using semidefinite programming, SIAM review, 49 (2007), pp. 434–448.
[22] E. de Klerk and D. V. Pasechnik, Approximation of the stability number of a graph via
copositive programming, SIAM Journal on Optimization, 12 (2002), pp. 875–892.
[23] E. De Klerk and R. Sotirov, Exploiting group symmetry in semidefinite programming relaxations of the quadratic assignment problem, Mathematical Programming, 122 (2010),
pp. 225–246.
[24] L. Ding and L. H. Lim, Higher-order cone programming. Available also as a Master’s thesis,
University of Chicago, 2016.
[25] A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri, Distinguishing separable and entangled
states, Physical Review Letters, 88 (2002).
[26] A. Domahidi, E. Chu, and S. Boyd, ECOS: An SOCP solver for embedded systems, in
Preecedings of the European Control Conference, IEEE, 2013, pp. 3071–3076.
[27] M. Dür, Copositive Programming–a Survey, in Recent advances in optimization and its applications in engineering, Springer, 2010, pp. 3–20.
[28] H. Fawzi, On representing the positive semidefinite cone using the second-order cone, arXiv
preprint arXiv:1610.04901, (2016).
[29] M. R. Garey and D. S. Johnson, Computers and Intractability, W. H. Freeman and Co., San
Francisco, Calif., 1979.
[30] K. Gatermann and P. A. Parrilo, Symmetry groups, semidefinite programs, and sums of
squares, Journal of Pure and Applied Algebra, 192 (2004), pp. 95–128.
[31] S. A. Gershgorin, Uber die Abgrenzung der Eigenwerte einer Matrix, Bulletin de l’Académie
des Sciences de l’URSS. Classe des sciences mathématiques et na, (1931), pp. 749–754.
[32] M. R. Gupta, A. Cotter, J. Pfeifer, K. Voevodski, K. Canini, A. Mangylov,
W. Moczydlowski, and A. Van Esbroeck, Monotonic calibrated interpolated look-up
tables, Journal of Machine Learning Research, 17 (2016), pp. 1–47.
[33] W. Habicht, Uber die Zerlegung strikte definiter Formen in Quadrate, Comment. Math. Helv.,

DSOS AND SDSOS OPTIMIZATION

37

12 (1940), pp. 317–322.
[34] L. Hannah and D. B. Dunson, Multivariate convex regression with adaptive partitioning.,
Journal of Machine Learning Research, 14 (2013), pp. 3261–3294.
[35] L. A. Hannah and D. B. Dunson, Ensemble methods for convex regression with applications
to geometric programming based circuit design, arXiv preprint arXiv:1206.4645, (2012).
[36] J. Harrison, Verifying nonlinear real formulas via sums of squares, in Theorem Proving in
Higher Order Logics, Springer, 2007, pp. 102–118.
[37] D. Hilbert, Über die Darstellung Definiter Formen als Summe von Formenquadraten, Math.
Ann., 32 (1888).
[38] C. Hillar and L.-H. Lim, Most tensor problems are NP-hard, arXiv preprint arXiv:0911.1393,
(2009).
[39] M. Huneault and F. Galiana, A survey of the optimal power flow literature, IEEE Transactions on Power Systems, 6 (1991), pp. 762–770.
[40] Z. Jarvis-Wloszek, R. Feeley, W. Tan, K. Sun, and A. Packard, Some controls applications of sum of squares programming, in Proceedings of the 42nd IEEE Conference on
Decision and Control, vol. 5, IEEE, 2003, pp. 4676–4681.
[41] H. Khalil, Nonlinear Systems, Prentice Hall, 2002. Third edition.
[42] J.-L. Krivine, Anneaux préordonnés, Journal d’analyse mathématique, 12 (1964), pp. 307–326.
[43] R. Laraki and J. B. Lasserre, Semidefinite programming for min–max problems and games,
Mathematical programming, 131 (2012), pp. 305–332.
[44] J. B. Lasserre, Global optimization with polynomials and the problem of moments, SIAM
Journal on Optimization, 11 (2001), pp. 796–817.
[45] J. B. Lasserre, Moments, Positive Polynomials and Their Applications, vol. 1, World Scientific, 2009.
[46] J. B. Lasserre, K.-C. Toh, and S. Yang, A bounded degree sos hierarchy for polynomial
optimization, EURO Journal on Computational Optimization, (2015), pp. 1–31.
[47] M. Laurent, Sums of squares, moment matrices and optimization over polynomials, in Emerging applications of algebraic geometry, Springer, 2009, pp. 157–270.
[48] M. Laurent, M. E. Nagy, and A. Varvitsiotis, Complexity of the positive semidefinite
matrix completion problem with a rank constraint, Discrete Geometry and Optimization,
69 (2013), pp. 105–120.
[49] X. Li, D. Sun, and K.-C. Toh, QSDPNAl: A two-phase proximal augmented Lagrangian
method for convex quadratic semidefinite programming, arXiv preprint arXiv:1512.08872,
(2015).
[50] L. Liberti, C. Lavor, N. Maculan, and A. Mucherino, Euclidean distance geometry and
applications, SIAM Review, 56 (2014), pp. 3–69.
[51] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, Applications of second-order cone
programming, Linear algebra and its applications, 284 (1998), pp. 193–228.
[52] L. Lovász, On the Shannon capacity of a graph, IEEE Transactions on Information theory, 25
(1979), pp. 1–7.
[53] A. Majumdar, A. Ahmadi, and R. Tedrake, Control and verification of high-dimensional
systems with dsos and sdsos programming, in 53rd IEEE Conference on Decision and
Control, IEEE, 2014, pp. 394–401.
[54] I. R. Manchester, M. M. Tobenkin, M. Levashov, and R. Tedrake, Regions of attraction
for hybrid limit cycles of walking robots, Proceedings of IFAC, 44 (2011), pp. 5801–5806.
[55] J. Mattingley and S. Boyd, Real-time convex optimization in signal processing, Signal Processing Magazine, IEEE, 27 (2010), pp. 50–61.
[56] J. Mattingley and S. Boyd, CVXGEN: a code generator for embedded convex optimization,
Optimization and Engineering, 13 (2012), pp. 1–27.
[57] A. Megretski, Systems polynomial optimization tools (spot), 2010.
[58] A. Mosek, The MOSEK optimization software, Online at http://www. mosek. com, 54 (2010),
pp. 2–1.
[59] T. S. Motzkin, The arithmetic-geometric inequality, in Inequalities (Proc. Sympos. WrightPatterson Air Force Base, Ohio, 1965), Academic Press, New York, 1967, pp. 205–224.
[60] K. G. Murty and S. N. Kabadi, Some NP-complete problems in quadratic and nonlinear
programming, Mathematical Programming, 39 (1987), pp. 117–129.
[61] Y. Nesterov, Squared functional systems and optimization problems, in High performance
optimization, vol. 33 of Appl. Optim., Kluwer Acad. Publ., Dordrecht, 2000, pp. 405–440.
[62] J. Nie and M. Schweighofer, On the complexity of Putinar’s Positivstellensatz, Journal of
Complexity, 23 (2007), pp. 135–150.
[63] J. Nie and L. Wang, Regularization methods for SDP relaxations in large-scale polynomial
optimization, SIAM Journal on Optimization, 22 (2012), pp. 408–428.

38

AMIR ALI AHMADI AND ANIRUDHA MAJUMDAR

[64] P. Parrilo and A. Jadbabaie, Approximation of the joint spectral radius using sum of squares,
Linear Algebra and its Applications, 428 (2008), pp. 2385–2402.
[65] P. A. Parrilo, Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization, PhD thesis, California Institute of Technology, May 2000.
[66] P. A. Parrilo, Semidefinite programming relaxations for semialgebraic problems, Mathematical Programming, 96 (2003), pp. 293–320.
[67] P. A. Parrilo, Polynomial games and sum of squares optimization, in Proceedings of the 45th
IEEE Conference on Decision and Control, 2006.
[68] F. Permenter and P. Parrilo, Partial facial reduction: simplified, equivalent SDPs via
approximations of the PSD cone, arXiv preprint arXiv:1408.4685, (2014).
[69] G. Pólya, Über Positive Darstellung von Polynomen, Vierteljschr. Naturforsch. Ges. Zürich,
73 (1928), pp. 141–145.
[70] I. Popescu, A semidefinite programming approach to optimal-moment bounds for convex
classes of distributions, Mathematics of Operations Research, 30 (2005), pp. 632–657.
[71] V. Powers, Positive polynomials and sums of squares: theory and practice, Real Algebraic
Geometry, (2011).
[72] M. Putinar, Positive polynomials on compact semi-algebraic sets, Indiana University Mathematics Journal, 42 (1993), pp. 969–984.
[73] B. Reznick, Uniform denominators in Hilbert’s 17th problem, Math Z., 220 (1995), pp. 75–97.
[74] B. Reznick, Some concrete aspects of Hilbert’s 17th problem, in Contemporary Mathematics,
vol. 253, American Mathematical Society, 2000, pp. 251–272.
[75] B. Reznick, Blenders, in Notions of Positivity and the Geometry of Polynomials, Springer,
2011, pp. 345–373.
[76] C. Riener, T. Theobald, L. J. Andrén, and J. B. Lasserre, Exploiting symmetries in SDPrelaxations for polynomial optimization, Mathematics of Operations Research, 38 (2013),
pp. 122–141.
[77] M. Roozbehani, Optimization of Lyapunov invariants in analysis and implementation of
safety-critical software systems, PhD thesis, Massachusetts Institute of Technology, 2008.
[78] N. V. Sahinidis, BARON: A general purpose global optimization software package, Journal of
global optimization, 8 (1996), pp. 201–205.
[79] N. Z. Shor, Class of global minimum bounds of polynomial functions, Cybernetics, 23 (1987),
pp. 731–734. (Russian orig.: Kibernetika, No. 6, (1987), 9–11).
[80] G. Stengle, A Nullstellensatz and a Positivstellensatz in semialgebraic geometry, Mathematische Annalen, 207 (1974), pp. 87–97.
[81] J. F. Sturm, Using SeDuMi 1.02, a Matlab toolbox for optimization over symmetric cones,
Optimization Methods and Software, 11 (1999), pp. 625 – 653.
[82] R. Tae, B. Dumitrescu, and L. Vandenberghe, Multidimensional FIR filter design via
trigonometric sum-of-squares optimization, IEEE Journal of Selected Topics in Signal Processing, 1 (2007), pp. 641–650.
[83] U. Topcu, A. Packard, P. Seiler, and G. Balas, Robust region-of-attraction estimation,
IEEE Transactions on Automatic Control, 55 (2010), pp. 137 –142.
[84] F. Vallentin, Symmetry in semidefinite programs, Linear Algebra and its Applications, 430
(2009), pp. 360–369.
[85] L. Vandenberghe and S. Boyd, Semidefinite programming, SIAM Review, 38 (1996), pp. 49–
95.
[86] X.-Y. Zhao, D. Sun, and K.-C. Toh, A Newton-CG augmented Lagrangian method for
semidefinite programming, SIAM Journal on Optimization, 20 (2010), pp. 1737–1765.
[87] H. Zou, T. Hastie, and R. Tibshirani, Sparse principal component analysis, Journal of
computational and graphical statistics, 15 (2006), pp. 265–286.

