COMPLEXITY OF FIRST-ORDER METHODS FOR DIFFERENTIABLE CONVEX OPTIMIZATION

This is a short tutorial on complexity studies for differentiable convex optimization. A complexity study is made for a class of problems, an "oracle" that obtains information about the problem at a given point, and a stopping rule for algorithms. These three items compose a scheme, for which we study the performance of algorithms and problem complexity. Our problem classes will be quadratic minimization and convex minimization in ℝn. The oracle will always be first order. We study the performance of steepest descent and Krylov spacemethods for quadratic function minimization and Nesterov’s approach to the minimization of differentiable convex functions.


INTRODUCTION
Due to the huge increase in the size of problems tractable with modern computers, the study of problem complexity and algorithm performance became essential. This was recognized very early by computer scientists and mathematicians working on combinatorial problems, and has recently become a central issue in continuous optimization. Complexity studies for these problems started in the former Soviet Union, and the main results are described in the book by Nemirovski & Yudin [14].
The special case of Linear Programming, which will not be tackled in this paper, initiated with Khachiyan [10], also in Russia in 1978, and had an explosive expansion in the West with the creation of interior point methods in the 80's and 90's. This paper starts with a brief introduction to the main concepts in the study of algorithm performance and complexity, following Nemirovski and Yudin, and then apply to the study of the convex optimization problem: minimize where f : R n → R is a continuously differentiable function.
In Section 2 we introduce the general framework for the study of algorithm performance and problem complexity and present a simple example.
We dedicate Section 3 to study the special case of convex quadratic functions, because they are the simplest non-linear functions: if a method is inefficient for quadratic problems, it will certainly be inefficient for more general problems; if it is efficient, it has a good chance of being adaptable to general differentiable convex problems, because near an optimal solution the quadratic approximation of the function uses to be precise. We study the performance of steepest descent and of Krylov space methods.
Section 4 will describe and analyze a basic method for unconstrained convex optimization devised by Nesterov [15], with "accelerated steepest descent" iterations. This method has become very popular, and the presentation and complexity proofs will be based on our paper [7].
Finally, we comment in Section 5 on improvements of this basic algorithm, presenting without proofs its extension to problems restricted to "simple sets" (sets onto which projecting a vector is easy).

SCHEMES, PERFORMANCE AND COMPLEXITY
A complexity study is associated with a scheme ( , O, τ ε ) as follows.
(i) is a class of problems. Examples: linear programming problems, unconstrained minimization of convex functions.
(ii) O is an oracle associated to .
The oracle is responsible for accessing the available information O(x) about a given problem in at a given point x.
• The numerical performance is a bound on the number of arithmetical operations needed in the worst case. The numerical performance is usually proportional to the iteration performance for each given algorithm. In this paper we only study the iteration performance of algorithms.
• The complexity of the scheme ( , O, τ ε ) is the performance of the best possible algorithm for the scheme. It is frequently unknown, and finding it for different schemes is the main purpose of complexity studies.
The performance of any algorithm for a scheme gives an upper bound to its complexity. A lower bound to the complexity may sometimes be found by constructing an example of a (difficult) problem in and finding a lower bound for any algorithm based on the same oracle and stopping rule. This will be the case in the end of Section 3.
If the performance of an algorithm for a scheme matches its complexity or a fixed multiple of it, it is called an optimal algorithm for the scheme.

First-order algorithms:
In most of this paper we study first-order algorithms for solving the problem minimize where f : R n → R is a differentiable function. The problem classes will be the special cases of quadratic and convex functions.
A first-order algorithm starts from a given point x 0 and constructs a sequence (x k ) using the where span(S) stands for the subspace generates by S.
In particular, the most well-known minimization algorithm is the steepest descent method, in which where λ k is a steplength. Each different choice of steplength (the rules of the method) defines a different steepest descent algorithm. This will be studied ahead in this paper.
Remark: In our algorithm model we used a single oracle, but there may be more than one.
and O 0 may be called more than once in each iteration. This is the case when line searches are used. The performance evaluation must then be adapted.

The notation O(·).
Given two real positive functions g(·) and h(·), we say that g = O(h) if there exists some constant K > 0 such that g(·) ≤ K h(·). This notation is very useful in complexity studies. For example, we shall prove that a certain steepest descent algorithm stops for k ≤ C 4 log 1 ε , where C is a parameter that identifies the problem in and ε is the precision. We may write k = C O(log(1/ε)), ignoring the coefficient 1/4.

Example: root of a continuous function
Here we present a simple example to illustrate how a complexity analysis works. Consider the following example of ( , O, τ ε ) given by: : Given a continuous function f : Remark: The stopping rule above is obviously not computable. We shall use a practical rule that implies τ ε .

Algorithm 2. Bisection
Performance: the following facts are straightforward at all iterations: If the algorithm does not stop at iteration k, then 2 −k > ε, and then k < log 2 (1/ε). We conclude that the stopping rule will be satisfied for k = log 2 (1/ε) , where r is the smallest integer above r. Thus that the performance of the scheme above is k = log 2 (1/ε) = O(log(1/ε)).
It is possible to prove that this is the best possible algorithm for this scheme, and hence the complexity of the scheme is this performance.

Remarks:
(i) Note that the only assumption here was the continuity of f . With stronger assumptions (Lipschitz constants for instance), better algorithms are described in numerical calculus textbooks.
(ii) The rules of the method are in the bisection calculation. Only the present oracle information O(m) is used at step k.

MINIMIZATION OF A QUADRATIC FUNCTION: FIRST-ORDER METHODS
Quadratic functions are the simplest nonlinear functions, and so an efficient algorithm for minimizing nonlinear functions must also be efficient in the quadratic case. On the other hand, near a minimizer, a twice differentiable function is usually well approximated by a quadratic function. A quadratic function is defined by where c ∈ R n and H is an n × n symmetric matrix. Then for x ∈ R n , If x * is a minimizer or a maximizer of f , then and hence finding an extremal of f is equivalent to solving the linear system H x * = −c, one of the most important problems in Mathematics.
The behavior of a quadratic function depends on the eigenvalues of its Hessian H . Since H is symmetric, it is known that H has n real eigenvalues which may be associated with n orthonormal (mutually orthogonal with unit norm) eigenvectors v 1 , v 2 , . . . , v n . There are four cases, represented in Figure 1:  (iv) If c T v i < 0 for some i = 1, . . . , k, then f is unbounded below.
In this section, we study the following scheme: : the class of quadratic functions that are bounded below (cases (i) and (iii) above). A function in has at least one minimizer x * . Without loss of generality, the study of algorithmic properties (not the implementation) may assume that x * = 0, and so the function becomes τ ε : Given an initial point x 0 ∈ R n , two rules will be used in the analysis: These rules are not implementable, because they require the knowledge of x * , but are very useful in the performance analysis.
Simplification: As we explained above, we assume that f (·) has a minimizer x * = 0. After performing the analysis with this simplification, we substitute x − x * for x. A further simplification may be done by diagonalizing H , also without loss of generality.

Steepest descent algorithms
In the first half of the 19 th century, Cauchy found that the gradient ∇ f (x) of a function f is the direction of maximum ascent of f from x, and stated the gradient method. It is the most basic of all optimization algorithms, and its performance is still an active research topic.

Algorithm 3. Steepest descent algorithm (model)
Data: Steplengths: Each different method for choosing the steplengths λ k defines a different steepest descent algorithm. Let us describe the two best known choices for the steplengths: • The Cauchy step, or exact step, the unique minimizer of f along the direction −g with g = ∇ f (x k ). The steplength is computed by setting ∇ f (x k − λ k g)⊥g and simplifying, which results in • The short step: λ k < 2/μ n , a fixed steplength.

Complexity results
Now we study the iteration performance of the steepest descent methods with these two steplength choices for minimizing a strictly convex quadratic function (case (i)). Given ε > 0 and In both cases the algorithm stops in O(C log(1/ε)) iterations, where C = μ n /μ 1 . At this moment the following question is open: find a steepest descent algorithm (by a different choice of λ k ) with performance O( √ C log(1/ε)). This performance is achieved in practice for "normal problems" (but not for particular worst case problems) by Barzilai-Borwein and spectral methods, described in [3]. Theorem 1. Let C = μ n /μ 1 ≥ 1 be the condition number of H . The iteration performance of the steepest descent method with Cauchy steplength for minimizing f starting at x 0 ∈ R n and with stopping criterion (5) is given by Proof. We begin by stating a classical result for the steepest descent step, which is based on the Kantorovich inequality and is proved for instance in [12, p. 238], Using this recursively, we obtain is an increasing function and that for t > 1, If τ ε is not satisfied at an iteration k, then by Example. In this example we show that the bound obtained in Theorem 1 is sharp, i.e., it cannot be improved. Take the following problem in R 2 : Assume that the initial point of some iteration has the shape x = (C, 1) z, for some z ∈ R. Then Computing the steplength λ by (4), we obtain λ = 2/(C + 1), and then the next iterate will be It follows that and this will be repeated on all iterations, with the worst possible performance as in Theorem 1.

Theorem 2.
Let C = μ n /μ 1 ≥ 1 be the condition number of H . The iteration performance of the steepest descent method with short steps λ k = 1/μ n , for minimizing f starting at x 0 ∈ R n and with stopping criterion (5) is given by

Proof.
A simplification in the analysis can be made by diagonalizing the matrix H by using the orthonormal matrix whose columns are the eigenvectors of H . Thus, we can consider, without loss of generality, that H = diag(μ 1 , . . . , μ n ). Given x 0 ∈ R n , by the steepest descent algorithm with short steps, Thus, for all i = 1, . . . , n, Consequently, by the definition of f , Proceeding like in the proof of Theorem 1, we obtain So, if τ ε is not satisfied at an iteration k, then k < C 2 log 1 ε , completing the proof.

Remarks:
• When the short steps 1/μ n are used, the result is that for • The diagonalization of H can be made without loss of generality for the performance analysis, as we did in the proof of Theorem 2. This leads to an interesting observation about the constant C, for the case in which there are null eigenvalues (case(iii)). Assuming that μ 1 = μ 2 = . . . = μ p = 0, we see that for i = 1, . . . , p, (∇ f (x)) i = 0, and the variables x i remain constant forever having no influence on the performance. The bounds in Theorems 1 and 2 remain valid for C = μ n /μ p+1 .

Krylov methods
Krylov space methods are the best possible algorithms for minimizing a quadratic function using only first-order information. Let us describe the geometry of a Krylov space method for the quadratic (2).
• Starting at a point x 0 , define the line This is actually the Cauchy step. We may write and the next iterate will be This is a two-dimensional problem.
• k-th step: adding ∇ f (x k−1 ) to the set of gradients, we construct the set and the next point will be It is then clear that x n is an optimal solution because V n = R n . This gives us a first performance bound k ≤ n for the Krylov space method. This bound is bad for high-dimensional spaces.
Main question: how to solve (P k ). Without proof (see for instance [15,20]), it is known that the directions (x k − x k−1 ) are conjugate, and any conjugate direction algorithm like Fletcher-Reeves [4,18] solves (P k ) at each iteration with about the same work as in the steepest descent method.
From now on we do a performance analysis of the Krylov space method, with stopping criterion A result for the relative error bound will also be discussed in the end of the section. The analysis is quite technical and will result in (16). Definition 1. Given x 0 ∈ R n and k ∈ N, define the k-th Krylov space by Consider V k = x 0 + K k and define the sequence (x k ) by Let P k be the set of polynomials p : R → R of degree k such that p(0) = 1, i.e, P k = 1 + a 1 t + a 2 t 2 + · · · + a k t k | a i ∈ R, i = 1, . . . , k .
From now on we deal with matrix polynomials, setting t = H .

Lemma 1.
A point x ∈ V k if, and only if, x = p(H )x 0 for some polynomial p ∈ P k . Furthermore,

Proof.
A point x ∈ V k if, and only if, where p ∈ P k . Furthermore, As H is symmetric, p(H ) T H = H p(H ), completing the proof.

Lemma 2. For any polynomial p ∈
Proof. Consider an arbitrary polynomial p ∈ P k . From Lemma 1, the point x = p(H )x 0 belongs to V k . As Using (8) we complete the proof.
Proof. As A is a symmetric matrix, there exists an orthogonal matrix P such that A = P D P T , with D = diag(λ 1 , λ 2 , . . . , λ n ). If q(t ) = a 0 + a 1 t + · · · + a k t k , then q(H ) = a 0 I + a 1 P D P T + · · · + a k (P D P T ) k = P a 0 I + a 1 D + · · · + a k D k P T .

Chebyshev Polynomials
The Chebyshev polynomials will be needed in the performance analysis of Krylov methods.
The next lemma shows that T k is, in fact, a polynomial (even though it does not look like one).

Lemma 5.
If T k (t ) = a k t k + · · · + a 2 t 2 + a 1 t + a 0 , then a k = 2 k−1 . Furthermore, Proof. We prove by induction. The results are trivial for k = 0 and k = 1. Suppose that the results hold for all natural number less than or equal to k. Using the induction hypothesis, consider leading to the first statement. Suppose that (k + 1) is even. Then k is odd and (k − 1) is even. Thus, by induction hypothesis, T k has only odd powers of t and T k−1 has only even powers. In this way, by (9), T k+1 has only even powers of t . Furthermore, its independent term is

COMPLEXITY OF FIRST-ORDER METHODS FOR DIFFERENTIABLE CONVEX OPTIMIZATION
On the other hand, if (k + 1) is odd, then k is even and (k − 1) is odd. Again by the induction hypothesis, T k only has even powers of t and T k−1 has only odd powers. Thus by (9), T k+1 has only odd powers of t . Furthermore, its linear term is completing the proof.
The next lemma discusses a relationship between a Chebyshev polynomial of odd degree and polynomials of the set P k , defined in (7). Lemma 6. Consider L > 0 and k ∈ N. Then there exists p ∈ P k such that, for all t ∈ [0, L], Proof. By Lemma 5, for all t ∈ [−1, 1], we have where the polynomial in parentheses has only even powers of t . So, for all t ∈ [0, L], Defining we complete the proof.

Complexity results
Now we present the main result about the performance of Krylov methods for minimizing a convex quadratic function. This result is based on [19, Thm. 3, p. 170].
We use the matrix norm defined by Theorem 3. Let μ n be the largest eigenvalue of H and consider the sequence (x k ) defined by (6).
and f (x k ) − f * ≤ ε is satisfied for Proof. Without loss of generality, assume that x * = 0. By Lemma 2, for all polynomial p ∈ P k , But, from Lemma 3 and (10), Considering the polynomial p ∈ P k given in Lemma 6 and using the fact that all eigenvalues of H belong to (0, μ n ], we have proving (11). If τ ε is not satisfied at an iteration k, then f (x k ) > ε and consequently ε < μ n x 0 2 2(2k + 1) 2 < μ n x 0 2 8k 2 which implies (12) and completes the proof.
Performance of the method for the relative error bound: A similar analysis for τ ε given by (5), also using Chebyshev polynomials, can be done using the condition number C. This is done in [14,23], and the result is clearly better than the best performance of the steepest descent algorithm for the steplength rules studied above, and for reasonable values of μ 1 , better than (12).
Complexity bound. The Krylov space methods uses at each iteration all the information gathered in the previous steps, and hence it seems to be the best possible algorithm based on first order information. In fact, Nemirovskii & Yudin [14] prove that no algorithm using a first order oracle can have a performance more than twice as good as the Krylov space method.
For methods based on accumulated first order information there is a negative result described by Nesterov [15, p. 59]: he constructs a quadratic problem (which he calls "the worst problem in the world") for which such methods need at least iterations to reach τ ε .
We conclude that the best performance for a first order method must be between the bounds (12) and (15). So the complexity of the scheme is

CONVEX DIFFERENTIABLE FUNCTIONS: THE BASIC ALGORITHM
In this section we study the performance of algorithms for the unconstrained minimization of differentiable convex functions. Quadratic functions are a particular case, and hence the performance bounds for first order algorithms will not be better than those found in the former section.
The role played by μ n in quadratic functions will be played by a Lipschitz constant L for the gradient of f (indeed, for a quadratic function the largest eigenvalue is a Lipschitz constant for the gradient), and we shall see that there are optimal algorithms, i.e., algorithms with the performance given by (16) with μ n replaced by L. These algorithms were developed by Nesterov [15], and are also studied in our papers [7,8].
Consider the scheme ( , O, τ ε ) where : the class of minimization problems of a convex continuously differentiable function f , with a Lipschitz constant L > 0 for the gradient. It means that for all x, y ∈ R n , Simple quadratic functions. The following definition will be useful in our development: we shall call "simple" a quadratic function φ : R n → R with ∇ 2 φ(x) = γ I , γ ∈ R, γ > 0. The following facts are easily proved for such functions: • φ(·) has a unique minimizer v ∈ R n (which we shall refer as the center of the quadratic), and the function can be written as • Given and

The algorithm
We now state the main algorithm and then study its properties. We include in the statement of the algorithm the definitions of the relevant functions (approximations of f (·) and the simple quadratic defined below).
We begin by summarizing the geometrical construction at an iteration k, represented in Figure 2. The iteration starts with two points x k , v k ∈ R n and a simple quadratic function whose global minimizer is v k . A point y k = x k + α(v k − x k ) is chosen between x k and v k . The choice of α is a central issue, and will be discussed later. All the action is centered on y k , with the following construction: • Take a gradient step from y k , generating x k+1 .

Analysis of the algorithm
The most important procedure in the algorithm is the choice of the parameter α k , which then determines y k at each iteration. The choice of α k is the one devised by Nesterov in [15, Scheme (2.2.6)]. Instead of "discovering" the values for these parameters, we shall simply adopt them and show their properties.
Once y k is determined, two independent actions are taken: (i) A steepest descent step from y k computes x k+1 .
(ii) A new simple quadratic is constructed by combining φ k (·) and the linear approximation Our scope will be to prove two facts: From these facts we shall conclude that f (x k ) → f * with the same speed as γ k → 0, which easily leads to the desired performance result.
The first lemma shows our main finding about the geometry of these points. All the action happens in the two-dimensional space defined by x k , v k , v k+1 . Note the beautiful similarity of the triangles in Figure 3.

Lemma 7.
Consider the sequences generated by Algorithm 4. Then for k ∈ N, Proof. By the algorithm, we know that Lα 2 k = γ k+1 , and completing the proof.

Lemma 8. Consider the sequences generated by Algorithm 4. Then for k ∈ N,
Proof. By the definition of φ α k , Using this, the definition of and the fact that α k ∈ (0, 1), we have By the definition of y k in the algorithm, v Substituting this in (21), we complete the proof.

Lemma 9. Consider the sequences generated by Algorithm 4. Then for k ∈ N,
Proof. The first inequality follows trivially from the convexity of f and the definition of u.
Since x k+1 and v k+1 are respectively global minimizers of u(·) and φ α k (·), we have from (18) that, for all x ∈ R n , As f (y k ) = u(y k ) and, from the last lemma, u(y k ) ≤ φ α k (v k ), we only need to show that The construction is shown in Fig. 3: since, by Lemma 7, Using this, (24) and the fact that by construction, proving the second inequality of (22).

Lemma 10.
For any x ∈ R n and k ∈ N, Proof. By the definition of φ α k and the fact that Subtracting f (x) in both sides, using (23) and the definition of γ k+1 , we have Using this recursively, we get the result and complete the proof.

Complexity
The following lemma was proved by Nesterov [15, p. 77] with a different notation.
Theorem 4. Consider the sequences generated by Algorithm 4 and assume that x * is an optimal solution. Then for k ∈ N, Proof. From Lemma 10,(25) holds in particular at x * , Using the fact that f (x k ) = φ k (v k ) ≤ φ k (x * ) and the definition of φ 0 , we get, Since x * is a minimizer of the convex function f , Applying this and the result of Lemma 11 in (28), As γ 0 = L, we have (26). If τ ε is not satisfied at an iteration k, then f (x k ) − f * > ε and consequently which implies (27) and completes the proof.
So, the iteration performance of Algorithm 4 is which corresponds to the complexity (16) for quadratic performance. Then, the algorithm is optimal.

CONVEX DIFFERENTIABLE FUNCTIONS: ENHANCED ALGORITHMS
The algorithm presented in the former section may be extended in several ways: the need for a previous knowledge of a Lipschitz constant L may be eliminated, a strong convexity constant akin to the smallest eigenvalue in the quadratic case may be used, and the algorithm may be extended to problems constrained to so-called simple sets. These extensions are treated in our paper [8] and in references therein.
In this section we state the extension of the basic algorithm to problems with simple constraints, without proofs. Consider the scheme ( , O, τ ε ) where : the class of problems where ⊂ R n is a closed convex set and f : R n → R is convex and continuously differentiable, with a Lipschitz constant L > 0 for the gradient. We assume that is a "simple" set, in the following sense: given an arbitrary point x ∈ R n , an oracle is available to compute P (x) = argmin y∈ x − y , the orthogonal projection onto the set . We now state the basic algorithm for constrained problems, without proofs. We keep in the statement the definition of the functions used in the analysis made in [8].

Algorithm 5.
Data: Consider the sequences generated by Algorithm 5. Then, as proved in [8,Thm. 2.6], at any iteration k before stopping, and hence This expression is similar to (27). In fact, if x * is a global minimizer, then may be introduced in the last expression, retrieving (27).
Estimations of the Lipschitz constant. Both Algorithms 4 and 5 and the algorithms discussed by Nesterov in [15, Chapter 2] make explicit use of a Lipschitz constant L for the function gradient. In [16], Nesterov describes a method for a more general problem, easily applied to the situations studied in this paper. This method includes a scheme for estimating the Lipschitz constant. In [7,8], the authors eliminate the use of L at the cost of an extra imprecise line search, and obtain an algorithm which keeps the optimal complexity properties and also inherits the global convergence properties of the steepest descent method for general continuously differentiable optimization. Besides this, the algorithm takes advantage of the knowledge of the strong convexity constant for the function and develop in [7] an adaptive procedure for estimating it. In another context -constrained minimization of non-smooth homogeneous convex functions -Richtárik [21] uses an adaptive scheme for guessing bounds for the distance between a point x 0 and an optimal solution x * . This bound determines the number of subgradient steps needed to obtain a desired precision.
Extensions. Nesterov's approach is applied to penalty methods by Lan, Lu & Monteiro [11], and interior descent methods based on Bregman distances are described by Auslender & Teboulle [1]. This method has been generalized to a non-interior method using projections by Rossetto [22]. Results for higher order methods are discussed in Nesterov & Polyak [17]. Accelerated versions of first-order algorithms following Nesterov's approach were developed by Monteiro, Ortiz & Svaiter [13] and by Beck & Teboulle [2], with improved performance for benchmark problems.

CONCLUSIONS
In this paper we described what we believe to be the basic results in the study of algorithm performance and problem complexity for the minimization of convex functions, both unconstrained and with simple constraints.
Algorithms with proved low worst-case performance are not necessarily efficient for practical problems. Khachiyan's algorithm [10] for linear programming is very inefficient, but had a great impact on the development of both continuous and discrete optimization. Karmarkar's algorithm [9] for linear programming improved Khachiyan's performance bound, and his bound was again improved later (see [6]). The effort to improve complexity led to better algorithms, which are nowadays used for solving large scale linear and quadratic in many domains. In fact, the largest linear programming problem treated up to now had 1.1 billion variables and 380 million constraints, solved by Gondzio & Grothey [5] using an interior point algorithm.
The conjugate gradient algorithm (Krylov space method) has optimal performance for quadratic problems, but its extension to more general problems is not straightforward. It was superseded by quasi-Newton methods, which are more efficient for non-quadratic problems, but coincide with it in the convex quadratic case. Note that the conjugate gradient method was not motivated by the complexity study, which came later.
Accelerated gradient methods are now in the phase of development, and we are not aware of any extensive comparison with classical algorithms. Research in this field is presently very active, and it is not clear to what classes of problems this approach will be applied and which methods will be the winners in practical applications to large-scale problems.