Exploiting Sparsity in the Coefficient Matching Conditions in Sum-of-Squares Programming using ADMM
aa r X i v : . [ m a t h . O C ] M a y Exploiting Sparsity in the Coefficient MatchingConditions in Sum-of-Squares Programmingusing ADMM
Yang Zheng, Giovanni Fantuzzi and Antonis Papachristodoulou,
Senior Member, IEEE
Abstract —This paper introduces an efficient first-order methodbased on the alternating direction method of multipliers (ADMM)to solve semidefinite programs (SDPs) arising from sum-of-squares (SOS) programming. We exploit the sparsity of the coefficient matching conditions when SOS programs are formu-lated in the usual monomial basis to reduce the computationalcost of the ADMM algorithm. Each iteration of our algorithmrequires one projection onto the positive semidefinite cone andthe solution of multiple quadratic programs with closed-formsolutions free of any matrix inversion. Our techniques areimplemented in the open-source MATLAB solver SOSADMM.Numerical experiments on SOS problems arising from uncon-strained polynomial minimization and from Lyapunov stabilityanalysis for polynomial systems show speed-ups compared to theinterior-point solver SeDuMi, and the first-order solver CDCS.
Index Terms —Sum-of-squares, ADMM, large-scale problems.
I. I
NTRODUCTION C HECKING whether a given polynomial is nonnegativehas applications in many areas (see [1], [2] and the ref-erences therein). For example, the unconstrained polynomialoptimization problem min x ∈ R n p ( x ) is equivalent to max γ subject to p ( x ) − γ ≥ . (1)Moreover, the stability of an equilibrium x ∗ of a polynomialdynamical system ˙ x ( t ) = f ( x ( t )) , x ( t ) ∈ R n in a neighbour-hood D of x ∗ (we assume x ∗ = 0 without loss of generality)—a fundamental problem in control theory—can be establishedby constructing a polynomial V ( x ) (called Lyapunov function )that satisfies the polynomial inequalities [3] ( V ( x ) > , ∀ x ∈ D\{ } , − ˙ V ( x ) = −h∇ V ( x ) , f ( x ) i ≥ , ∀ x ∈ D . (2)Throughout this work, h· , ·i denotes the inner product in theappropriate Hilbert space.A powerful way to test polynomial inequalities is to employa sum-of-squares (SOS) relaxation (we refer the reader to [4],[5] for details on SOS relaxations in polynomial optimization, Y. Zheng is supported by the Clarendon Scholarship and the Jason HuScholarship. A. Papachristodoulou is supported by project EP/M002454/1.Y. Zheng and A. Papachristodoulou are with the Department of EngineeringScience, University of Oxford, Parks Road, Oxford, OX1 3PJ, U.K. (e-mail: { yang.zheng, antonis } @eng.ox.ac.uk)G. Fantuzzi is with the Department of Aeronautics, Imperial CollegeLondon, South Kensington Campus, London, SW7 2AZ, United Kingdom(e-mail: [email protected]) and to [3] for a tutorial on SOS techniques for systems analy-sis). In fact, while testing the non-negativity of a polynomialis NP-hard in general, the existence of an SOS decompositioncan be checked in polynomial time by solving a semidefiniteprogram (SDP) [4]. Unfortunately, however, the size of theSDP for the SOS relaxation of a degree- d polynomial in n variables is (cid:0) n + dd (cid:1) . Consequently, SOS relaxations are limitedto small problem instances; with the current technology, forexample, Lyapunov-based analysis is impractical for generalsystems with ten or more states.In order to mitigate scalability issues, one can act at the modeling level, i.e. one can try to replace the SDP obtainedfrom an SOS relaxation with an optimization problem that ischeaper to solve still using second-order interior-point methods(IPMs), implemented in efficient solvers such as SeDuMi [6].One approach is to exploit structural properties of the polyno-mial whose positivity is being tested [7]–[11]. For example,computing the Newton polytope [7] or checking for diagonalinconsistency [8] can restrict the monomial basis required inthe SOS decomposition by eliminating redundant monomials.Further improvements are possible by group-theoretic sym-metry reduction techniques [9] and graph-theoretic correlativesparsity [10]. Facial reduction has also been applied to select areduced monomial basis for SOS programs in [11]. A secondapproach is to approximate the positive semidefinite (PSD)cone using diagonally dominant or scaled diagonally dominantmatrices [12], [13]. These relaxations can be solved with lin-ear programs (LPs) or second-order-cone programs (SOCPs),rather than SDPs, and the conservativeness introduced byapproximating the PSD cone can be reduced with a recentlyproposed basis pursuit algorithm [14].Further improvements are available on the computational level if IPMs are replaced by more scalable first-order methods(FOMs) at the cost of reduced accuracy. The design of efficientfirst-order algorithms for large-scale SDPs has received partic-ular attention in recent years. For instance, Wen et al. proposedan alternating direction augmented Lagrangian method forlarge-scale dual SDPs [15]. O’Donoghue et al. developedan operator-splitting method to solve the homogeneous self-dual embedding of conic programs [16], which has recentlybeen extended by the authors to exploit aggregate sparsityvia chordal decomposition [17], [18]. In the context of SOSprogramming, Bertsimas et al. proposed an accelerated FOMfor unconstrained polynomial optimization [19], while Henrion& Malick introduced a projection-based method for SOSrelaxations [20]. However, both approaches are only applicableo a small subset of SOS programs as they rely on theconstraint of the corresponding SDP being orthogonal, whichis not the case for SOS problems with free variables.In this paper, we propose a first-order algorithm based onthe alternating direction method of multipliers (ADMM) tosolve the SDPs arising from SOS optimization. In contrastto [19] and [20], we exploit the sparsity in the coefficientmatching conditions, making our approach suitable for a largerclass of SOS programs. While the aggregate sparsity patternof these SDPs is dense (so that the methods of [17], [18] arenot very advantageous), each equality constraint in the SDPonly involves a small subset of decision variables when anSOS program is formulated in the usual monomial basis. Thissparsity can be exploited to formulate an efficient ADMMalgorithm, the iterations of which consist of conic projectionsand optimization problems with closed-form solutions that—crucially—are free of any matrix inversion. We implementour techniques in SOSADMM, an open-source MATLABsolver. The efficiency of our methods compared to the IPMsolver SeDuMi [6] and the first-order solver CDCS [21] isdemonstrated on SOS problems arising from unconstrainedpolynomial optimization and from Lyapunov stability analysisof polynomial systems.The rest of this paper is organized as follows. Section IIreviews SOS polynomials and the ADMM algorithm. Sparsityfor SDPs arising in SOS programs is discussed in Section III,and we show how to exploit it to build an efficient ADMMalgorithm in Section IV. Numerical experiments are reportedin Section V. Section VI concludes the paper.II. P RELIMINARIES
A. SOS polynomials and SDPs
The sets of real and natural numbers (including zero) aredenoted by R and N , respectively. Let x ∈ R n , α ∈ N n , andlet x α = x α x α · · · x α n n denote a monomial in x of degree | α | = P ni =1 α i . Given an integer d ∈ N , we denote N nd = { α ∈ N n : | α | ≤ d } , and the vector of all monomials ofdegree no greater than d by v d ( x ) = { x α | α ∈ N nd } = [1 , x , x , . . . , x n , x , x x , . . . , x dn ] T . (3)The length of v d ( x ) is | N nd | = (cid:0) n + dd (cid:1) . A real polynomial p ( x ) is a finite, real linear combination of monomials of xp ( x ) = X α ∈ N n p α x α , p α ∈ R . The degree of p ( x ) is the maximum of the degrees of allmonomials with nonzero coefficients. We denote the set ofreal polynomials in x by R [ x ] . Definition 1:
A polynomial p ( x ) ∈ R [ x ] of degree d is asum-of-squares (SOS) if there exist polynomials f i ( x ) ∈ R [ x ] , i = 1 , . . . , m of degree no greater than d such that p ( x ) = m X i =1 [ f i ( x )] . Clearly, the existence of an SOS representation guaranteesthat p ( x ) ≥ . The following theorem gives an equivalentcharacterization of SOS polynomials. Proposition 1 ([4]):
A polynomial p ( x ) ∈ R [ x ] of degree d is an SOS polynomial if and only if there exists a (cid:0) n + dd (cid:1) × (cid:0) n + dd (cid:1) symmetric PSD matrix X (cid:23) such that p ( x ) = v d ( x ) T Xv d ( x ) . (4)The equality in (4) gives a set of affine equalities on theelements of X to match the coefficients of p ( x ) . Togetherwith X (cid:23) , this makes the problem of finding an SOSrepresentation for p ( x ) an SDP. The formulation of such SDPscan be assisted by software packages, such as SOSTOOLS[22] and GloptiPoly [23]. Remark 1:
The size of the PSD matrix X in (4) is (cid:0) n + dd (cid:1) × (cid:0) n + dd (cid:1) because we have used the full set of monomialsof degree no greater than d in our representation. This numbermight be reduced by inspecting the structural properties of p ( x ) to identify and eliminate redundant monomials in v d ( x ) ;well-known techniques include Newton polytope [7], diagonalinconsistency [8], symmetry property [9], and facial reduc-tion [11]. B. ADMM algorithm
The ADMM algorithm solves the optimization problem min y,z f ( y ) + g ( z ) subject to Ay + Bz = c, (5)where y ∈ R n and z ∈ R m are the decision variables, f : R n → R and g : R m → R are convex functions, and A ∈ R l × n , B ∈ R l × m and c ∈ R l are the constraint data. Given apenalty parameter ρ > and a multiplier λ ∈ R l (known asthe dual variable ), the ADMM algorithm solves (5) by findinga saddle point [24, Chapter 5.4] of the augmented Lagrangian L ρ ( y, z, λ ) = f ( y ) + g ( z ) + λ T ( Ay + Bz − c )+ ρ k Ay + Bz − c k (6)with the following steps: y k +1 = arg min y L ρ ( y, z k , λ k ) , (7a) z k +1 = arg min z L ρ ( y k +1 , z, λ k ) , (7b) λ k +1 = λ k + ρ ( Ay k +1 + Bz k +1 − c ) . (7c)In these equations, the superscript k denotes the value of avariable at the k -th iteration of the algorithm, and k · k denotesthe standard Euclidean norm, i.e. , k x k = √ x T x for x ∈ R n .Then, from a computational perspective, steps (7a) and (7b)are equivalent to the minimizations of ˜ L ρ ( y, z, λ ) = f ( y ) + g ( z ) + ρ (cid:13)(cid:13)(cid:13)(cid:13) Ay + Bz − c + 1 ρ λ (cid:13)(cid:13)(cid:13)(cid:13) (8)over y and z , respectively, with λ fixed. Under very mildconditions, ADMM converges to a solution with a rate O ( k ) ,which is independent of ρ , although its value can affectconvergence in practice; see [25, Section 3.2] for details. ABLE ID
ENSITY OF NONZERO ELEMENTS IN THE EQUALITY CONSTRAINTS OF
SDP (12) n d = 4 1 . × − . × − . × − . × − . × − . × − . × − d = 6 4 . × − . × − . × − . × − . × − . × − . × − d = 8 2 . × − . × − . × − . × − . × − . × − . × − III. R OW S PARSITY IN
SDP
S FROM
SOS
PROGRAMS
A. SDP formulations of SOS relaxations
Let A α be the indicator matrix for the monomials x α inthe rank-one matrix v d ( x ) v d ( x ) T ; in other words, the entry of A α with row index β and column index γ (where the naturalordering for multi-indices β, γ ∈ N nd is used) satisfies ( A α ) β,γ = ( if β + γ = α otherwise . (9)The SOS constraint (4) can then be reformulated as p ( x ) = h v d ( x ) v d ( x ) T , X i = X α ∈ N n d h A α , X i x α . (10)Matching the coefficients of the left- and right-hand sides givesthe equality constraints h A α , X i = p α ∀ α ∈ N n d . (11)We refer to these equalities as coefficient matching conditions .The existence of an SOS decomposition for p ( x ) (or lackthereof) can then be checked with the feasibility SDPfind X subject to h A α , X i = p α , α ∈ N n d ,X (cid:23) . (12)When the full monomial basis is used, as in this case, thedimension of X and the number of constraints in (12) are,respectively, N = | N nd | = (cid:18) n + dd (cid:19) , m = | N n d | = (cid:18) n + 2 d d (cid:19) . (13) B. Properties of the coefficient matching conditions
In this section, for simplicity, we re-index the constraintmatching conditions (11) using integer indices i = 1 , . . . , m instead of the multi-indices α .The conditions (11) inherit two important properties fromthe data matrices A i , i = 1 , . . . , m . The first one follows fromthe fact that the matrices A i are orthogonal. If n i denotes thenumber of nonzero entries in A i we have h A i , A j i = ( n i if i = j, otherwise . (14)After letting vec : S N → R N be the usual operator mappinga matrix to the stack of its columns, and defining A = (cid:2) vec ( A ) · · · vec ( A m ) (cid:3) T , (15)the equality constraints in (12) can be rewritten as the matrix-vector product A · vec ( X ) = b , where b ∈ R m is a vectorcollecting the coefficients p i , i = 1 , . . . , m . Property (14) Fig. 1. Sparsity pattern of AA T for the example (16) directly implies the following lemma, which forms the basisof the FOMs of [19], [20]. Lemma 1 (Orthogonality of constraints): AA T is an m × m diagonal matrix with ( AA T ) ii = n i .The second property of the coefficient matching conditionsis that they are sparse, in the sense that each equality constraintin (12) only involves a small subset of entries of X , becauseonly a small subset of entries of the product v d ( x ) v d ( x ) T areequal to a given monomial x α . Thus, the vectorized matrix A is row sparse , meaning that each row is a sparse vector. Inparticular, the following result holds. Lemma 2 (Sparsity of constraints):
Let A be the vectorizedmatrix for (12), and let N and m be as in (13). The numberof nonzero elements in A is N , and the density of nonzeroelements in A is equal to m − = O ( n − d ) . Proof:
Since the matrix v d ( x ) v d ( x ) T contains all mono-mials x α , α ∈ N n d , all entries of the PSD matrix X enter atleast one of the equality constraints in (12). Moreover, (14)implies that each entry of X enters at most one constraint.Therefore, A must contain N nonzero elements. Its densityis then given by N N × m = 1 m = (cid:20)(cid:18) n + 2 d d (cid:19)(cid:21) − = O ( n − d ) . Remark 2:
While the constraint matrix A is sparse (seetypical values in Table I), the aggregate sparsity pattern of theSDP (12) is dense because all entries of the matrix variable X appear in the equality constraints. This implies that X isgenerally a dense variable, so the first-order algorithms of [17],[18] are not particularly suitable. Remark 3:
The property of orthogonality in Lemma 1 holdsfor standard SOS feasibility problems. However, this propertyfails for the following example:find a, b subject to ax + bx + x + 1 is SOS ,bx + ax + x + 1 is SOS . (16)Fig. 1 shows the sparsity pattern of AA T for (16) obtainedusing SOSTOOLS, demonstrating that the constraints are notorthogonal. The reason is that (16) involves the free parameters a, b as well as the PSD matrices for the SOS representation,nd this destroys the orthogonality. This issue is common incontrol applications; see, e.g. , the condition (2) when findingLyapunov functions. Consequently, the first-order algorithmsin [19], [20] cannot be applied to many problems with SOSconstraints because they rely on constraint orthogonality.IV. E XPLOITING R OW S PARSITY IN
SDP S As we have seen, the algorithms in [17]–[20] are not usefulfor generic SOS problems because their aggregate sparsitypattern is dense, and the orthogonality property only holds forsimple SOS feasibility problems. However, the data matrix A is always row-sparse due to the coefficient matching conditionsin SOS programs. This property can be exploited to constructan efficient ADMM algorithm. In the following, we considera generic SDP in the vectorized form min x c T x subject to Ax = b, x ∈ K , (17)where x is the optimization variable, A ∈ R m × ˆ n , b ∈ R m and c ∈ R ˆ n are the problem data, and K is a product of cones, atleast one of which is the PSD cone. Throughout this section, δ S ( x ) denotes the indicator function of a set S , δ S ( x ) = ( , if x ∈ S , + ∞ , if x / ∈ S . A. Reformulation considering individual row sparsity
Let us represent A = [ a , a , . . . , a m ] T , so each vector a i is a row of A , and let H i ∈ R | a i |× ˆ n , i = 1 , . . . , m be “entry-selector” matrices of 1’s and 0’s selecting the nonzero elementsof a i , where | a i | denotes the number of nonzero elements of a i . Note that the rows of H i are orthonormal, since each selectsa different entry of a i . Then, Ax = b ⇔ ( ( H i a i ) T z i = b i , i = 1 , . . . , m,z i = H i x, i = 1 , . . . , m. (18)In (18), z i is a copy of the elements of x which enter the i -th affine constraint. It is also convenient to introduce anadditional slack variable u = x , so the affine constraints in (18)and conic constraint in (17) are decoupled when applying theADMM algorithm. We can then reformulate (17) as min z i ,u,x c T x subject to ( H i a i ) T z i = b i i = 1 , . . . , m,z i = H i x, i = 1 , . . . , m,u = x, u ∈ K . (19) B. ADMM steps
To apply ADMM, we move the affine constraints ( H i a i ) T z i = b i and the conic constraint u ∈ K in (19) tothe objective using the indicator functions δ ( · ) and δ K ( · ) ,respectively: min z i ,u,x c T x + δ K ( u ) + m X i =1 δ (cid:0) ( H i a i ) T z i − b i (cid:1) subject to z i = H i x, i = 1 , . . . , m,u = x. (20) The augmented Lagrangian of (20) is L = c T x + δ K ( u ) + m X i =1 δ (cid:2) ( H i a i ) T z i − b i (cid:3) + m X i =1 µ Ti ( z i − H i x ) + ρ k z i − H i x k + ξ T ( u − x ) + ρ k u − x k , (21)and we group the variables as Y = { x } , Z = { u, z , . . . , z m } , D = { µ , . . . , µ m , ξ } . Then, the ADMM steps (7a)–(7c) become the following:
1) Minimization over Y : The minimization of (21) overthe variables in Y is equivalent to an unconstrained quadraticprogram, min x c T x + ρ m X i =1 (cid:13)(cid:13)(cid:13)(cid:13) z ki − H i x + µ ki ρ (cid:13)(cid:13)(cid:13)(cid:13) + ρ (cid:13)(cid:13)(cid:13)(cid:13) u k − x + ξ k ρ (cid:13)(cid:13)(cid:13)(cid:13) . (22)The updated variable x k +1 is then simply given by x k +1 = D − " m X i =1 H Ti (cid:18) z ki + µ ki ρ (cid:19) + (cid:18) u k + ξ k ρ (cid:19) − ρ c , (23)where the matrix D = I + P mi =1 H Ti H i is diagonal because therows of each matrix H i are orthonormal. This means that (23)is cheap to calculate.
2) Minimization over Z : Minimizing (21) over the vari-ables in Z amounts to a conic projection, min u (cid:13)(cid:13) u − x k +1 + ρ − ξ k (cid:13)(cid:13) subject to u ∈ K , (24)plus m independent quadratic programs min z i (cid:13)(cid:13) z i − H i x k +1 + ρ − µ ki (cid:13)(cid:13) subject to ( H i a i ) T z i = b i . (25)The projection (24) is easy to compute when K is a productof R n , the non-negative orthant, second-order cones, and PSDcones; for example, a projection onto the PSD cone onlyrequires one eigen-decomposition. As for problem (25), itsKKT conditions are z i − H i x k +1 + ρ − µ ki + ( H i a i ) ω i = 0 , (26a) ( H i a i ) T z i = b i , (26b)where ω i is the Lagrangian multiplier for the equality con-straint in (25). Simple algebra shows that ω i = 1 k H i a i k (cid:18) − b i + ( H i a i ) T H i x k +1 − ρ ( H i a i ) T µ ki (cid:19) , so the solution z k +1 i to (25) can be calculated easily with (26a).Note that this step is free of any matrix inversion.
3) Update multipliers D : According to (7c), the multipliersin D are updated with inexpensive and parallelizable gradientascent steps: µ k +1 i = µ ki + ρ ( z k +1 i − H i x k +1 ) , i = 1 , . . . , m,ξ k +1 = ξ k + ρ ( u k +1 − x k +1 ) . (27) . Summary of the computations in the ADMM algorithm In the proposed ADMM algorithm, subproblems (7a) and(7b) have explicit closed-form solutions. Each iteration re-quires solving1) one unconstrained quadratic program, given by (22);2) one conic projection, given by (24);3) m independent quadratic programs, given by (25).Note that only the nonzero elements of a i appear in (25).Since we have assumed that a i is sparse, only operations onvectors of small size are required. Besides, our algorithm isfree of matrix inversion (with the exception of the m × m diagonal matrix D , requiring O ( m ) flops), which results fromintroducing the local variables z i so each affine constraintcan be considered individually. The cost is that our algorithmneeds to maintain multiple local variables z i , which may haveadverse effects on the convergence speed of ADMM.In contrast, the FOMs in [19], [20] fail to deal with generalSOS programs since they rely on orthogonality of constraints,and those in [17], [18] require factorizing the m × m matrix AA T (in general, O ( m ) flops). In SDPs arising from genericSOS relaxations, this step is computationally demanding be-cause AA T is not diagonal (as seen in Remark 3, this is oftenthe case) and the number m is usually large ( m = 18564 if n = 12 and d = 6 in (12)). Note that all these algorithmsconverge at rate O ( k ) because they are based on ADMM orits variants. V. N UMERICAL E XPERIMENTS
We implemented our techniques in SOSADMM, an open-source first-order MATLAB solver for conic programs withrow sparsity. Currently, SOSADMM supports cartesian prod-ucts of the following cones: R n , non-negative orthant, second-order cone, and the positive definite cone. SOSADMM, thenumerical examples presented in this section, and additionalexamples are available from https://github.com/oxfordcontrol/SOSADMM We tested SOSADMM on random unconstrained polyno-mial optimization problems and Lyapunov stability analysisof polynomial systems. To assess the suboptimality of thesolution returned by SOSADMM, we compared it to the accu-rate one computed with the interior-point solver SeDuMi [6].CPU times were compared to the first-order solver CDCS [21],which exploits aggregate sparsity in SDPs; in particular, theprimal method in CDCS was used [17]. In our experiments,the termination tolerance for SOSADMM and CDCS was setto − , and the maximum number of iterations was .To improve convergence, SOSADMM employs an adaptivepenalty parameter update rule [25], with an initial value ρ = 1 .All tests were run on a PC with a 2.8 GHz Intel R (cid:13) Core TM i7CPU and 8GB of RAM. A. Unconstrained polynomial optimization
Consider the global polynomial minimization problem min x ∈ R n p ( x ) , (28)where p ( x ) is a given polynomial. This problem is equivalentto (1), and we can obtain an SDP relaxation by replacing the TABLE IICPU
TIME ( S ) TO SOLVE THE
SDP
RELAXATIONS OF (28). N IS THE SIZEOF THE
PSD
CONE , m IS THE NUMBER OF CONSTRAINTS .Dimensions CPU time (s) n N m
SeDuMi CDCS(primal) SOS-ADMM
15 69 0.13 0.11 0.06
28 209 0.24 0.16 0.14
45 494 1.16 0.18 0.18
66 1000 3.17 0.25 0.39
91 1819 13.89 0.46 0.55
120 3059 54.63 0.79 0.84
153 4844 187.0 0.92 0.82
190 7314 610.2 2.91 1.92
231 10625 1739 4.93 2.32TABLE IIIL
YAPUNOV FUNCTIONS FOR THE SYSTEM (29)Solver Time (s) Lyapunov function V ( x ) SeDuMi 0.054 . x + 4 . x + 2 . x CDCS-primal 0.21 . x + 1 . x + 2 . x SOSADMM 0.58 . x + 1 . x + 2 . x non-negativity constraint with an SOS condition on p ( x ) − γ .Motivated by [20], we generated p ( x ) according to p ( x ) = p ( x ) + n X i =1 x di , where p ( x ) is a random polynomial of degree strictly lessthan d . We used GloptiPoly [23] to generate the examples.Table II compares the CPU time (in seconds) requiredto solve the SOS relaxation as the number n of variableswas increased with d = 2 . Both SOSADMM and CDCS-primal were faster than SeDuMi on these examples (note thatSeDuMi’s runtime reduces if a weaker termination toleranceis set, but not significantly). Also, the optimal value returnedby SOSADMM was within 0.05% of the high-accuracy valuereturned by SeDuMi. For all examples in Table II, the cone size N is moderate (less than 300), while the number of constraints m is large. SeDuMi assembles and solves an m × m linearsystem at each iteration, which is computationally expensive. B. Finding Lyapunov functions
Next, we consider the problem of constructing Lyapunovfunctions to check local stability of polynomial/rational sys-tems when (2) is replaced by SOS conditions. We usedSOSTOOLS [22] to generate the corresponding SDPs.The first system we study is ˙ x = − x − x x , ˙ x = − x − x x , ˙ x = − x − x x + 1 + 3 x x , (29)which is demo 2 in SOSTOOLS. The system has an equilib-rium at the origin, and we search for a homogeneous quadraticpolynomial Lyapunov function V ( x ) = ax + bx + cx toprove its global stability. The results given by SeDuMi, CDCS-primal and SOSADMM are listed in Table III. For such a small ABLE IVCPU
TIME ( S ) TO CONSTRUCT A QUADRATIC L YAPUNOV FUNCTION FORRANDOMLY GENERATED POLYNOMIAL SYSTEMS .Statistics CPU time (s) n Size of A nonzerodensity SeDuMi CDCS(primal) SOS-ADMM10 × . × − × . × − × . × − × . × − × . × − × . × − × . × − × . × − * 729.1 104.730 × . × − * 3509.2 259.0 * SeDuMi fails due to memory requirements. system, SeDuMi was slightly faster than CDCS-primal andSOSADMM, which is expected since IPMs are well-suited forsmall-scale SDPs. Note that since the problem of constructinga Lyapunov functions is a feasibility problem, the solutionsreturned by SeDuMi, CDCS-primal and SOSADMM need notbe the same (see Table III).As the last example, we consider randomly generatedpolynomial dynamical systems ˙ x = f ( x ) of degree threewith a locally asymptotically stable equilibrium at the ori-gin, and checked for local nonlinear stability in the ball D = { x ∈ R n | . − k x k ≥ } using a complete quadraticpolynomial as the candidate Lyapunov function. Table IVsummarizes the average CPU times required to search for sucha Lyapunov function, when successful (note that we cannotdetect infeasible problems because we only solve the primalform (17)). The results clearly show that SOSADMM is fasterthan both SeDuMi and CDCS-primal for the largest probleminstances ( n ≥ ). Also, FOMs have much lower memoryrequirements, and SOSADMM can solve problems that arenot accessible with IPM: SeDuMi failed due to memory issueswhen n > . Finally, note that for the problem of findingLyapunov functions the m × m linear system solved in SeDuMiand CDCS is not diagonal, and solving it is expensive: when n = 30 it took over 150 s for CDCS just to factorize AA T ,which is over 50% of the total time taken by SOSADMM toreturn a solution. VI. C ONCLUSION
In this paper, we proposed an efficient ADMM algorithmto exploit the row-sparsity of SDPs that arise from SOSprogramming, which are implemented in SOSADMM. Thesubproblems of our algorithm consist of one conic projectionand multiple quadratic programs with closed-form solutions,which can be computed efficiently and—most importantly—do not require any matrix inversion.Our numerical experiments on random unconstrained poly-nomial optimization and on Lyapunov stability analysis ofpolynomial/rational systems demonstrate that our method canprovide speed-ups compared to the interior-point solver Se-DuMi and the first-order solver CDCS. One major drawback ofour method is the inability to detect infeasibility; future workwill try to exploit the sparsity of SDPs from SOS relaxations in a homogeneous self-dual embedding formulation similar tothat of [16], [18]. R
EFERENCES[1] J. B. Lasserre,
Moments, positive polynomials and their applications .World Scientific, 2009, vol. 1.[2] G. Chesi and D. Henrion, “Guest editorial: Special issue on positivepolynomials in control,”
IEEE Trans. Autom. Control , vol. 54, no. 5,pp. 935–936, 2009.[3] A. Papachristodoulou and S. Prajna, “A tutorial on sum of squarestechniques for systems analysis,” in
Am. Control Conf. (ACC) . IEEE,2005, pp. 2686–2700.[4] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraicproblems,”
Math. Program. , vol. 96, no. 2, pp. 293–320, 2003.[5] J. B. Lasserre, “Global optimization with polynomials and the problemof moments,”
SIAM J. Optim. , vol. 11, no. 3, pp. 796–817, 2001.[6] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimizationover symmetric cones,”
Optim. Methods Softw. , vol. 11, no. 1-4, pp.625–653, 1999.[7] B. Reznick et al. , “Extremal PSD forms with few terms,”
Duke Math.J. , vol. 45, no. 2, pp. 363–374, 1978.[8] J. L¨ofberg, “Pre-and post-processing sum-of-squares programs in prac-tice,”
IEEE Trans. Autom. Control , vol. 54, no. 5, pp. 1007–1011, 2009.[9] K. Gatermann and P. A. Parrilo, “Symmetry groups, semidefinite pro-grams, and sums of squares,”
J. Pure Appl. Algebra , vol. 192, no. 1, pp.95–128, 2004.[10] H. Waki, S. Kim, M. Kojima, and M. Muramatsu, “Sums of squares andsemidefinite program relaxations for polynomial optimization problemswith structured sparsity,”
SIAM J. Control Optim. , vol. 17, no. 1, pp.218–242, 2006.[11] F. Permenter and P. A. Parrilo, “Basis selection for SOS programs viafacial reduction and polyhedral approximations,” in , 2014, pp. 6615–6620.[12] A. Majumdar, A. A. Ahmadi, and R. Tedrake, “Control and verificationof high-dimensional systems with DSOS and SDSOS programming,” in
IEEE 53rd Conf. Decis. Control (CDC) , 2014, pp. 394–401.[13] A. A. Ahmadi and A. Majumdar, “DSOS and SDSOS optimization: LPand SOCP-based alternatives to sum of squares optimization,” in . IEEE, 2014, pp. 1–5.[14] A. A. Ahmadi and G. Hall, “Sum of squares basis pursuit with linearand second order cone programming,” arXiv preprint arXiv:1510.01597 ,2015.[15] Z. Wen, D. Goldfarb, and W. Yin, “Alternating direction augmentedLagrangian methods for semidefinite programming,”
Math. Program.Comput. , vol. 2, no. 3-4, pp. 203–230, 2010.[16] B. O’Donoghue, E. Chu, and S. Parikh, Nealand Boyd, “Conic opti-mization via operator splitting and homogeneous self-dual embedding,”
J. Optim. Theory Appl. , vol. 169, no. 3, pp. 1042–1068, 2016.[17] Y. Zheng, G. Fantuzzi, A. Papachristodoulou, P. Goulart, and A. Wynn,“Fast ADMM for semidefinite programs with chordal sparsity,” arXivpreprint arXiv:1609.06068 , 2016.[18] ——, “Fast ADMM for homogeneous self-dual embedding of sparseSDPs,” arXiv preprint arXiv:1611.01828 , 2016.[19] D. Bertsimas, R. M. Freund, and X. A. Sun, “An accelerated first-order method for solving SOS relaxations of unconstrained polynomialoptimization problems,”
Optim. Methods Softw. , vol. 28, no. 3, pp. 424–441, 2013.[20] D. Henrion and J. Malick, “Projection methods in conic optimization,”in
Handbook on Semidefinite, Conic and Polynomial Optimization .Springer, 2012, pp. 565–600.[21] Y. Zheng, G. Fantuzzi, A. Papachristodoulou, P. Goulart, andA. Wynn, “CDCS: Cone decomposition conic solver, version 1.1,”https://github.com/oxfordcontrol/CDCS, Sep. 2016.[22] A. Papachristodoulou, J. Anderson, G. Valmorbida, S. Prajna, P. Seiler,and P. Parrilo, “SOSTOOLS version 3.00 sum of squares optimizationtoolbox for MATLAB,” arXiv preprint arXiv:1310.4716 , 2013.[23] D. Henrion and J.-B. Lasserre, “GloptiPoly: Global optimization overpolynomials with MATLAB and SeDuMi,”
ACM Tans. Math. Softw.(TOMS) , vol. 29, no. 2, pp. 165–194, 2003.[24] S. Boyd and L. Vandenberghe,
Convex optimization . CambridgeUniversity Press, 2004.[25] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributedoptimization and statistical learning via the alternating direction methodof multipliers,”