{'fraction_non_alphanumeric': 0.05091315989289853, 'fraction_numerical': 0.012708308356312193, 'mean_word_length': 4.290486257928118, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Design of First-Order Optimization Algorithms via Sum-of-Squares Programming

In this paper, we propose a framework based on sum-of-squares programming to design iterative first-order optimization algorithms for smooth and strongly convex problems. Our starting point is to develop a polynomial matrix inequality as a sufficient condition for exponential convergence of the algorithm. The entries of this matrix are polynomial functions of the unknown parameters (exponential decay rate, stepsize, momentum coefficient, etc.). We then formulate a polynomial optimization, in which the objective is to optimize the exponential decay rate over the parameters of the algorithm. Finally, we use sum-of-squares programming as a tractable relaxation of the proposed polynomial optimization problem. We illustrate the utility of the proposed framework by designing a first-order algorithm that shares the same structure as Nesterov's accelerated gradient method.


I. INTRODUCTION
Many applications in science and engineering involve solving convex optimization problems using iterative methods. In the development of iterative optimization algorithms, there are several objectives that algorithm designers should consider [1]: • Robust performance guarantees: Algorithms should perform well for a wide range of optimization problems in their class. • Time and space efficiency: Algorithms should be efficient in terms of both time and storage, although these two may conflict. • Accuracy: Algorithms should be able to provide arbitrarily accurate solutions to the problem at a reasonable computational cost. In general, these goals may conflict. For example, a rapidly convergent method (e.g., Newton's method) may require too much computer storage and/or computation. In contrast, a robust method, resilient to noise and uncertainties, may also be too slow in finding an optimal solution. Tradeoffs between, for example, convergence rate and storage requirements, or between robustness and speed, are central issues in numerical optimization [1].
In recent years, there has been an increasing interest in using tools from control theory to study the convergence properties of iterative optimization algorithms. The connection between control and optimization is made by interpreting these iterative algorithms as discrete-time dynamical systems with feedback. This interpretation provides many insights and new directions of research. In particular, we can use control tools to study disturbance rejection properties of optimization algorithms [2]- [4], study robustness to parameter and model uncertainty [5], and analyze tracking and adaptation capabilities [6]. This interpretation also opens the door to the use of control tools for algorithm design [5], [7], [8].
The main aim of this paper is to develop an optimization framework, based on tools from robust control theory and polynomial optimization, to design first-order optimization algorithms for the class of smooth and strongly convex problems, in which the convergence rate is exponential (O(ρ k ) for some 0 ≤ ρ < 1). To this end, we start with a result in [8], in which the authors derive a matrix inequality as a sufficient condition for exponential stability of the algorithm. The entries of this matrix are polynomial functions of (i) the unknown parameters of the algorithm (e.g., stepsize, momentum coefficient) and (ii) the unknown exponential decay rate. We then formulate a polynomial optimization problem in which the cost function is the exponential decay rate ρ, the constraint is the polynomial matrix inequality described above, and the decision variables are the tunable parameters of the algorithm. Finally, we use sum-of-squares programming as a tractable relaxation of the polynomial optimization problem to tune the parameters of the algorithm. We illustrate our approach by designing a first-order method sharing the same structure as Nesterov's accelerated method. To the best of our knowledge, this is the first work on principled and computationally efficient numerical algorithm design.
The rest of the paper is organized as follows. In Section II, we use a Lyapunov analysis framework, in which we cast the problem of finding the worst-case exponential rate bound of a first-order optimization algorithm as a semidefinite program. In Section III, we turn to algorithm synthesis and use the results of Section II to formulate the algorithm design problem as a polynomial optimization and use sum-of-squares machinery to solve the algorithm design problem. Finally, we illustrate the performance of our design framework via numerical simulations.

II. ALGORITHM ANALYSIS A. Preliminaries
We denote the set of real numbers by R, the set of real n-dimensional vectors by R n , the set of m × n-dimensional matrices by R m×n , and the n-dimensional identity matrix by I n . We denote by S n , S n + , and S n ++ the sets of nby-n symmetric, positive semidefinite, and positive definite arXiv:1803.10928v2 [math.OC] 23 Sep 2018 matrices, respectively. We denote by vec(·) a linear transformation which converts the matrix into a column vector. For a function f : R d → R, we denote by dom f = {x ∈ R n : f (x) < ∞} the effective domain of f . The p-norm (p ≥ 1) is displayed by · p : R d → R + . For two matrices A ∈ R m×n and B ∈ R p×q of arbitrary dimensions, we denote their Kronecker product by A ⊗ B. We define R [x] as the polynomial ring in n variables and R [x] n,d as the polynomial in n variables of degree at most d.
holds for some L f > 0 and all x, y ∈ S. (1a) implies for all x, y ∈ S.
holds for some m f > 0 and all x, y ∈ S. An equivalent definition is that for all x, y ∈ S.
We denote the class of functions satisfying (1a) and (2a) by F(m f , L f ). A differentiable function f belongs to the class F(m f , L f ) on S if and only if the inequality holds for all x, y ∈ S [5], [9], where the indefinite matrix Q f ∈ S 2d is given by The

B. Algorithm Representation
Consider the following optimization problem where f ∈ F(m f , L f ) is smooth and strongly convex. Under this assumption, the well-known optimality condition of a point y is that Note that y is unique as f is strongly convex. Consider a first-order algorithm that generates a sequence of points {y k } that solves (4). In other words, we have lim where ∇f (y ) = 0. We can represent the algorithm in the following state-space form [5], [8]: where ξ k ∈ R n (n ≥ d) is the state of the algorithm and y k ∈ R d is the output at which the gradient is evaluated. Furthermore, we assume that x k is another output whose fixed point is optimal, i.e., lim k→∞ f (x k ) = f (x ) where x = y . Therefore, the fixed points of the algorithm must naturally satisfy In particular, one of the eigenvalues of A(θ) is equal to one. Note that the matrices A(θ) ∈ R n×n , B(θ) ∈ R n×d , C(θ) ∈ R d×n , and E(θ) ∈ R d×n in (5) are parameterized by the vector θ ∈ R p , which is the concatenation of the parameters of the algorithm (e.g., stepsize, momentum coefficient, etc.). We give two examples below.

Example 1
The gradient method is given by the recursion where h ≥ 0 is the stepsize. For this algorithm, θ = h is the only parameter, and the matrices A(θ), B(θ), C(θ) and E(θ) are given by Example 2 Consider the following recursion defined on the two sequences {x k } and {y k }, where h, β and γ are nonnegative scalars. By defining the state vector ξ k = [x k−1 x k ] ∈ R 2d and the parameter vector θ = [h β γ] , we can represent (9) in the canonical form (5), as follows, Therefore, the matrices A(θ), B(θ), C(θ) and E(θ) are given by Notice that depending on the selection of β and γ, (9) describes various existing algorithms. For example, the gradient method corresponds to the case β = γ = 0 (see Example 1). In Nesterov's accelerated method, we have β = γ [10]. Finally, we recover Heavy-ball method by setting γ = 0 [11]. In Table I, we provide various parameter selections and convergence rates for the gradient method, the heavyball method, and Nesterov's accelerated method [5], [10], [11].

C. Exponential Convergence via SDPs
To measure the progress of the algorithm in (5) towards optimality, we make use of the following Lyapunov function [8]: where x is the unique minimizer of f and P ∈ S n + is an unknown positive semidefinite matrix that does not depend on k. The first term in (12) is the suboptimality of x k and the second term is a weighted "distance" between the state ξ k and the fixed point ξ . Notice that by this definition, V k is positive everywhere and zero at optimality. Suppose we select P such that the Lyapunov function satisfies the inequality for some 0 ≤ ρ < 1. By iterating down (13) to k = 0, we can conclude that V k ≤ ρ 2k V 0 for all k. This implies that In other words, the algorithm exhibits an O(ρ 2k ) convergence rate-in terms of objective values-if the Lyapunov function satisfies the decrease condition (13). In the following result, developed in [8], we present a matrix inequality whose feasibility implies (13) and, hence, the exponential convergence of the algorithm. For notational convenience, we drop the argument θ wherever the dependence of the matrices A(θ), B(θ), C(θ), and E(θ) on θ is clear from the context.
. Consider the dynamical system in (5), whose fixed point satisfies (6). Define the matrices where Suppose there exists a positive semidefinite P ∈ S n + , and nonnegative scalars 0 < ρ ≤ 1, λ ≥ 0 that satisfy the following matrix inequality: Then the sequence {x k } satisfies Proof: See [8]. According to Theorem 1, any triple (ρ, P, λ) ∈ [0, 1] × S n + × R + that satisfies the matrix inequality in (16) certifies an O(ρ 2k ) convergence rate for the algorithm. In particular, the fastest convergence rate can be found by solving the following optimization problem: where θ is given (the algorithm parameters), and the decision variables are λ, ρ, P . Note that (18) is a quasiconvex program since the constraint is affine in λ and P for a fixed ρ. We can therefore use a bisectioning search to find the smallest possible value of the decay rate ρ. This approach has been pursued in [5] using a different Lyapunov function. Furthermore, the corresponding matrix inequality for proximal variants of (5) has been developed in [12]. Note that for finding an -accurate (0 < ≤ 1) solution for the decay rate ρ, the computational complexity is at most O(q 3 log 2 ( −1 )) where q, the total number of unknowns, is independent of the dimension d of the optimization problem (4) as we remark below.

Remark 1
We can often exploit some special structure in the matrices A, B, C, and E to reduce the dimension of (16). For many algorithms, these matrices are in the form where nowĀ,B,C, andĒ are lower dimensional matrices independent of d [5, §4.2]. By selecting P =P ⊗I d , whereP is a lower dimensional matrix, we can factor out all the Kronecker products ⊗I d from the matrices M 0 , M 1 , M 2 , M 3 and make the dimension of the corresponding matrix inequality in (16) independent of d. In particular, a multi-step method with r ≥ 1 steps yields an (r + 1) × (r + 1) matrix inequality. For instance, the gradient method (r = 1) and Nesterov's accelerated method (r = 2) yield 2 × 2 and 3 × 3 matrix inequalities, respectively.

Remark 2 (Non-monotone Algorithms)
We emphasize that in the development of Theorem 1 although we require the Lyapunov function to decrease geometrically at each iteration (see (13)), the resulting bound in (14) does not imply that the sequence of objective values is monotone, i.e., (17) does not imply f (x k+1 ) ≤ f (x k ). This allows us to analyze the convergence properties of nonmonotone algorithms such as Nesterov's accelerated method.
Remark 3 (Feasibility of (18)) The matrix inequality (16) provides an upper bound on the true decay rate and, therefore, (16) is sufficient for the exponential convergence result in (17). In other words, there might be an exponentially convergence algorithm for which (16) is not feasible. Nevertheless, it has been shown in [8] that the bounds are not conservative. For instance, for Nesterov's accelerated method, the rate bound obtained by solving (18) turns out to be even better than the theoretical rate bound proved by Nesterov [8].

III. ALGORITHM SYNTHESIS
We saw in the previous section that the exponential stability of a given first-order algorithm can be certified by solving an SDP feasibility problem. More precisely, given an algorithm in (5) with a prespecified value of θ (the tuning parameters of the algorithm), we can search for a suitable Lyapunov function and establish a rate bound for the algorithm by solving a quasiconvex program. A natural question to ask is whether we can leverage the same framework to do algorithm design. We formalize this problem as follows.
Given a parameterized family of first-order methods given by (5), tune the parameters θ of the algorithm, within a compact set Θ ⊂ R p , such that the resulting algorithm converges at an O(ρ 2k ) rate to x with a minimal ρ.  Using the result of Theorem 1, Problem 1 can be formally written as the following nonconvex optimization problem: subject to M (θ, ρ, λ, P ) 0 where the decision variables are now θ, P, λ, and ρ. We recall from (15) that the parameter θ appears in (19) through the matrices A, B, C, E, M 1 , M 2 , M 3 . Intuitively, (19) searches for the parameters θ of the algorithm for optimal performance (minimal ρ) while respecting the stability condition (13), which is imposed by the matrix inequality constraint in (19). Note that if we fix θ = θ 0 , which means the algorithm parameters are given, then (19) reduces to the quasiconvex program in (18). This suggests a natural but inefficient way to solve (19): We could do an exhaustive search over the parameter space Θ, and solve (18) to find the optimal ρ 2 for each value of θ. We, therefore, need to solve a sequence of quasiconvex programs in order to find the optimal tuning. As an illustration, we implement this approach to tune the parameters of the Nesterov's accelerated method (the algorithm in (9) with γ = β), where the tuning parameters are the stepsize h and the momentum coefficient β. In Figure  1 Note that, in general, the exhaustive search approach described above becomes prohibitively costly as the dimension and/or the granularity of the search space Θ increase. We therefore need an efficient way to solve (19). Although this problem is nonconvex, the special structure of the constraint set makes the problem tractable. To see this, we note that all the entries of the constraint matrix in (19) are polynomial functions of the decision variables. This matrix inequality constraint can be alternatively "scalarized" and rewritten in terms of scalar polynomial inequalities (e.g., by considering minors, or coefficients of the characteristic polynomial). The resulting optimization problem is of the form It is here that we can draw a direct connection from the algorithm design problem in (19) to polynomial optimization, which are tractable problems in many cases [13]. We briefly introduce polynomial optimization next.

A. Sum-of-Squares Programs
The main difficulty in solving problems involving polynomial constraints, such as the one in (19), is the lack of efficient numerical methods able to handle multivariate nonnegativity conditions. A computationally efficient approach is to use sum-of-squares (SOS) relaxations [14], [15]. In what follows, we introduce the basic results used in our derivations.
Definition 3 A multivariate polynomial of degree 2d in n variables with real coefficients, p(x 1 , x 2 , ..., x n ) = p(x) ∈ R[x] n,2d , is a sum-of-squares (SOS) if there exist polynomials q 1 (x), ..., q m (x) ∈ R[x] n,d such that We will denote the set of SOS polynomials in n variables of degree at most 2d by Σ n,2d [x]. A polynomial being an SOS is sufficient to certify its global nonnegativity, since any p ∈ Σ n,2d [x] satisfies p(x) ≥ 0 for all x ∈ R n . Hence, Σ n,2d ⊆ P n,2d , where P n,2d [x] is the set of nonnegative polynomials in R [x] n,2d . Given a polynomial p (x) ∈ R [x] n,2d , the existence of an SOS decomposition of the form (20) is equivalent to the existence of a positive semidefinite matrix , called the Gram matrix, such that where [x] d = 1, x 1 , . . . , x n , x 2 1 , x 1 x 2 , . . . , x d n is the vector of all n+d d monomials in x 1 , . . . , x n of degree at most d. Notice that the equality constraint in (21) is affine in the matrix Q, since the expansion of the right-hand side results in a polynomial whose coefficients depend affinely on the entries of Q and must be equal to the corresponding coefficients of the given polynomial p (x). Hence, finding an SOS decomposition is computationally equivalent to finding a positive semidefinite matrix Q subject to the affine constraint in (21), which is a semidefinite program [14], [16].
Using the notion of sos polynomials, we can now define the class of sum-of-squares programs (SOSP). An SOSP is an optimization program in which we maximize a linear function over a feasible set given by the intersection of an affine family of polynomials and the set Σ [x] of SOS polynomials in R [x], as described below [17]: where y = [y 1 · · · y m ] is the vector of decision variables, p i (x; y) = c i (x) + a i1 (x) y 1 + · · · + a im (x) y m , and c i , a ij are given multivariate polynomials in R [x]. Note that in the above optimization problem, x is the vector of indeterminates and not the decision variables. Despite their generality, it can be proved that SOSPs are equivalent to SDPs; hence, they are convex programs and can be solved in polynomial time [16]. In recent years, SOSPs have been used as convex relaxations for various computationally hard optimization and control problems (see, for example, [14], [15], [18]- [20] and the volume [21]).
The notion of positive definiteness and sum-of-squares of scalar-valued polynomials can be extended to polynomial matrices, i.e., matrices with entries in R[x]. The definition of an sos matrix is as follows [22].
Since an m × m matrix is simply a representation of an mvariate quadratic form, we can always interpret an sos matrix in terms of a polynomial with m additional variables. The following lemma makes this precise.
Lemma 1 Consider a symmetric matrix with polynomial entries P(x) ∈ R[x] m×m , and let z = [z 1 , . . . , z m ] T be a vector of indeterminates. Then P(x) is a sum-of-squares matrix (SOSM) if z T P(x)z is an SOS polynomial in R[x, z].
Obviously, a polynomial matrix P (x) being SOSM provides an explicit certificate for P (x) being positive semidefinite for all x ∈ R n .

B. Polynomial Optimization Problems
One application of SOSP is the global optimization of a polynomial p(x). To this end, rather than directly computing a minimizer x of p(x), we instead focus on obtaining the best possible lower bound on its optimal value p(x ). This viewpoint is based on the observation that a real-valued number γ is a global lower bound of p(x) if and only if the polynomial p(x) − γ is nonnegative for all x. The best lower bound on p(x) is thus obtained by solving the optimization problem where the decision variable is now γ and the original decision variable x acts as an indeterminate. Observe that (22) is a convex optimization problem with infinitely many constraints. By replacing the nonnegativity condition with an sos constraint, we obtain the following optimization problem Note that for multivariate polynomials, nonnegativity and sos are not equivalent. More precisely, we have the inclusion Σ n,2d ⊆ P n,2d . Therefore, since the feasible set of the second problem is a subset of the feasible set of the first problem, we have the inequality p sos ≤ p . However, for relatively small problems, we often have p sos = p , i.e., there is no loss of optimality by using sos relaxations. Nevertheless, even in those situations where p sos < p , we can improve the lower bound p sos by producing stronger sos conditions [13]. In our particular application, we need to optimize a multivariate polynomial over a set described by polynomial inequalities: where p (x), g i (x) ∈ R [x] for all i = 1, . . . , m. Similar to the unconstrained case, we can instead find the best lower bound on p(x) on the constraint set, as follows: The corresponding sos relaxation requires p(x) to be sos on S. Recalling the formal similarity with weak duality and Lagrange multipliers, it is natural to consider the following decomposition for p(x): where s 0 (x) and s i (x) are sos polynomials that are determined by matching the coefficients of the left-and righthand side. This particular decomposition implies that for any x 0 ∈ S, we have g i (x 0 ) ≥ 0 for all i and therefore, the condition p(x 0 ) ≥ 0 is automatically satisfied. Considering this decomposition, we obtain the SOCP problem Following the discussion in §III-A, the above problem is an SOCP and can be converted to an SDP.

IV. NUMERICAL SIMULATIONS
In this section, we validate the proposed approach on two examples: the gradient method and Nesterov's accelerated method. There are several software packages that convert polynomial optimization problems into a hierarchy of SDP relaxations and then call an external solver to solve them. In this paper, we use the software Gloptipoly3 [23], which is oriented toward global polynomial optimization problems of the form (23).

A. The Gradient Method
Consider the gradient method of Example 1. By choosing P = pI d (p ≥ 0), we can apply the dimensionality reduction outlined in Remark 1 and reduce the dimension of the LMI. After dimensionality reduction, the matrices M 0 , M 1 , M 2 , and M 3 in the LMI (16) read as By substituting these matrices back in (16), we obtain the following matrix inequality constraint: where the entries M ij , 1 ≤ i, j ≤ 2 of M are all polynomials functions of ρ 2 , h, λ, p and are given by Recalling that the condition M 0 is equivalent to Tr(M ) ≤ 0 and det(M ) ≥ 0, the design problem in (19) for the gradient method is equivalent to the following polynomial optimization problem: with decision variables ρ 2 , h, λ, and p. In Figure 2, we plot the optimal rate bound ρ sos , obtained by solving an sos relaxation of (25), for various values of the condition number κ f . We also plot the analytical rate outlined in Table I. We observe that the rate obtained from the sos formulation coincides with the analytical rate, which is known to be tight.

B. Nesterov's Accelerated Method
Consider the algorithm of Example 2 with γ = β, which corresponds to Nesterov's accelerated method. The statespace matrices of this algorithm are given in (11). By applying the dimensionality reduction of Remark 1, we arrive at the following matrices that appear in the LMI (15): where Q f is given in (3b) and P is now a 2 × 2 positive semidefinite matrix. By defining x = [ρ 2 h β λ vec(P ) ] ∈ R 7 , p(x) = ρ 2 , and M (x) = M 0 + ρ 2 M 1 + (1 − ρ 2 )M 2 + λM 3 ∈ S 3 , the corresponding optimization problem can be written as minimize p(x) subject to − M (x) 0.
Using Lemma 1 we can convert the polynomial matrix inequality −M (x) 0 to a single polynomial inequality in a higher-dimensional space (see Lemma 1). Here we scalarize the positive semidefiniteness constraint by considering the principal minors (Sylvester's criterion for positive semidefinite matrices). Furthermore, we fix the stepsize h to the value 1/L f for solving the sos relaxation. In other words, we optimize ρ 2 over (β, λ, P ).
In Figure 3, we plot the optimal decay rate ρ sos , obtained by solving an sos relaxation of (27), for various values of κ f . We observe that the obtained rate is slightly better than Nesterov's rate and worse than the theoretical lower bound, which is achieved by the Heavy-ball method on quadratic objective functions-see Table I.  The blue curve is obtained by tuning the parameters using the developed sos framework. The red curve is for the parameter tuning proposed by Nesterov [10] (see also Table  I). Finally, the yellow curve is the theoretical lower bound for all first-order methods, which is achieved by the Heavy-ball method on quadratic objective functions-see Table I.

V. CONCLUSIONS
We have considered the problem of designing first-order iterative optimization algorithms for solving smooth and strongly convex optimization problems. By using a family of parameterized nonquadratic Lyapunov functions, we presented a polynomial matrix inequality as a sufficient condition for exponential stability of the algorithm. All the entries of this matrix have a polynomial dependence on the unknown parameters, which are the parameters of the algorithm, the parameters of the Lyapunov function, and the exponential decay rate. We then formulated a polynomial optimization problem to search for the optimal convergence rate subject to the polynomial matrix inequality. Finally, we proposed a sum-of-squares relaxation to solve the resulting design problem. We illustrated the proposed approach via numerical simulations.