Performance of first-order methods for smooth convex minimization: a novel approach

Introduction

First-order convex optimization methods have recently gained in popularity both in theoretical optimization and in many scientific applications, such as signal and image processing, communications, machine learning, and many more. These problems are very large scale, and first-order methods, which in general involve very cheap and simple computational iterations, are often the best option to tackle such problems in a reasonable time, when moderate accuracy solutions are sufficient. For convex optimization problems, there exists an extensive literature on the development and analysis of first-order methods, and in recent years, this has been revitalized at a quick pace due to the emergence of many fundamental new applications alluded above, see e.g., the recent collections and references therein.

where $f^{\prime}(\cdot)$ stands for the gradient of $f(\cdot)$. Note that the algorithm has another implicit knowledge, i.e., that the distance from its initial point $x_{0}$ to a minimizer $x_{\ast}\in X_{\ast}(f)$ of $f$ is bounded by some constant $R>0$, see more precise definitions in the next section.

Given a desired accuracy $\varepsilon>0$, applying the given algorithm on the function $f$ in the class $\mathcal{F}$, the algorithm stops when it produces an approximate solution $x_{\varepsilon}$ which is $\varepsilon$-optimal, that is such that

The performance (or complexity) of a first order black-box optimization algorithm is then measured by the number of oracle calls the algorithm needs to find such an approximate solution. Equivalently, we can measure the performance of an algorithm by looking at the absolute inaccuracy

where $x_{N}$ is the result of the algorithm after making $N$ calls to the oracle. Throughout this paper we will use the latter form to measure the performance of a given algorithm.

Building on this model, in this work we introduce a novel approach for analyzing the performance of a given first order scheme. Our approach relies on the observation that by definition, the worst case behavior of a first-order black-box optimization algorithm is by itself an optimization problem which consists of finding the maximal absolute inaccuracy over all possible inputs to the algorithm. Thus, with $x_{N}$ being the output of the algorithm after making $N$ calls to the oracle, we look at the solution of the following Performance Estimation Problem (PEP):

At first glance this problem seems very hard or impossible to solve. We overcome this difficulty through an analysis that relies on various types of relaxations, including duality and semi-definite relaxation techniques. The problem and setting, and an outline of the underlying idea of the proposed approach for analyzing (P) are described in Section 2. In order to develop the basic idea and tools underlying our proposed approach, we first focus on the fundamental gradient method (GM) for smooth convex minimization, and then extend it to a broader class of first order black box minimization methods. Obviously, the gradient method is a particular case of this broader class that will be analyzed below. However, it is quite important to start with the gradient method for two reasons. First, it allows to acquaint the reader in a more transparent way with the techniques and methodolgy we need to develop in order to analyze (PEP), thus paving the way to tackle more general schemes. Secondly, for the gradient method, we are able to prove a new and tight bound on its performance which is given analytically, see Section 3. Capitalizing on the methodology and tools developed in the past section, in Section 4, we consider a broader class of first-order black-box methods, which among others, is shown to include the so-called heavy-ball and fast gradient schemes . Although an analytical solution is not available for this general case, we show that for this broader class of methods, it is always possible to compute numerical bounds for an adequate relaxation of the corresponding PEP, allowing to derive new bounds on the performance of these methods. We then derive in Section 5 an efficient procedure for finding optimal step sizes which results in a first-order method that achieves best performance. Our approach and analysis give rise to some interesting problems leading us to suggest some conjectures. Finally, an appendix includes the proof of a technical result.

The Problem and the Main Approach

Let $\mathcal{A}$ be a first-order algorithm for solving the optimization problem

Throughout the paper we make the following assumptions:

We assume that (M) is solvable, i.e., the optimal set $X_{*}(f):=\operatorname*{arg\,min}f$ is nonempty, and for $x_{*}\in X_{*}(f)$ we set $f^{*}:=f(x_{*}).$

There exists $R>0$, such that the distance from $x_{0}$ to an optimal solution $x_{\ast}\in X_{*}(f)$ is bounded by $R$.In general, the terms $L$ and $R$ are unknown or difficult to compute, in which case some upper bound estimates can be used in place. Note that all currently known complexity results for first order methods depend on $L$ and $R$.

2 Basic Idea and Main Approach

We are interested in measuring the worst-case behavior of a given algorithm $\mathcal{A}$ in terms of the absolute inaccuracy $f(x_{N})-f(x_{\ast})$, by solving problem (P) defined in the introduction, namely

To tackle this problem, we suggest to perform a series of relaxations thereby reaching a tractable optimization problem.

An un-formal description of the underlying idea consists of two main steps as follows:

To define constraints that relate the new variables, we use relevant/useful properties characterizing the family of convex functions in $C^{1,1}_{L}$, as well as the rule(s) describing the given algorithm ${\cal A}$.

This approach can, in principle, be applied to any optimization algorithm. Note that any relaxation performed on the maximization problem (P) may increase its optimal value, however, the optimal value of the relaxed problem still remains a valid upper bound on $f(x_{N})-f^{\ast}$.

A formal description on how this approach can be applied to the gradient method is described in the next section, which as we shall see, allows us to derive a new tight bound on the performance of the gradient method.

An Analytical Bound for the Gradient Method

To develop the basic idea and tools underlying the proposed approach for analyzing the performance of iterative optimization algorithms, in this section we focus on the simplest fundamental method for smooth convex minimization, the Gradient Method (GM). It will also pave the way to tackle more general first-order schemes as developed in the forthcoming sections.

Consider the gradient algorithm with constant step size, as applied to problem $(M)$, which generates a sequence of points as follows:

Here $h>0$ is fixed. At this point, we recall that for $h=1$, the convergence rate of the GM algorithm can be shown to be (see for example ):

To begin our analysis, we first recall a fundamental well-known property for the class of convex $C^{1,1}_{L}$ functions, see e.g., [14, Theorem 2.1.5].

and note that we always have $\delta_{\ast}=0$ and $g_{\ast}=0$.

In terms of $\delta_{i}$, $g_{i}$, condition (3.3) becomes

The obtained problem remains nontrivial to tackle. We will now perform some simplifications on this problem that will be useful for the forthcoming analysis.

Secondly, we consider (3.4) for the four cases $i=\ast$, $j=\ast$, $i<j$ and $j<i$, and use the equality constraints

to eliminate the variables $x_{1},\dots,x_{N}$. After some algebra, we reach the following form for the PEP:

where $i<j=0,\dots,N$ is a shorthand notation for $i=0,\dots,N-1,j=i+1,\dots,N$.

Finally, we note that the optimal solution for this problem is attained when $\|x_{\ast}-x_{0}\|=R$, and hence we can also eliminate the variables $x_{0}$ and $x_{\ast}$. This produces the following PEP for the gradient method, a nonconvex quadratic minimization problem:

Therefore, by defining the following $(N+1)\times(N+1)$ symmetric matrices

Problem (G) is a nonhomogeneous Quadratic Matrix Program, a class of problems recently introduced and studied by Beck .

2 A Tight Performance Estimate for the Gradient Method

We now proceed to establish the two main results of this section. First, we derive an upper bound on the performance of the gradient method, this is accomplished via using duality arguments. Then, we show that this bound can actually be attained by applying the gradient method on a specific convex function in the class $C^{1,1}_{L}$.

In order to simplify the following analysis, we will remove some constraints from (G) and consider the bound produced by the following relaxed problem:

As we shall show below, it turns out that this additional relaxation has no damaging effects and produces the desired performance bound when $0<h\leq 1$.

To establish our dual result, the next lemma shows that a dimension reduction is possible when minimizing a quadratic matrix function sharing the special form as the one that appears in problem (G′).

and using (3.7) it follows $\inf_{\xi}f(\xi b^{T})=f({\bar{\xi}}b^{T})=\|b\|^{2}a^{T}{\bar{\xi}}=\operatorname{tr}(ba^{T}{\bar{\xi}}b^{T})$. Now, using (3.6)-(3.7), one obtains ${\bar{X}}^{T}Q=-ba^{T}$ and $Q({\bar{X}}-{\bar{\xi}}b^{T})=0$, and hence it follows that

Equipped with Lemma 3.1, we now derive a Lagrangian dual for problem (G′).

For convenience, we recast (G′) as a minimization problem, and we also omit the fixed term $LR^{2}$ from the objective. That is, we consider the equivalent problem (G′′) defined by

The dual objective function is then defined by

and the dual problem of (G′′) is then given by

Since $L_{1}(\cdot,\lambda,\tau)$ is linear in $\delta$, we have $\min_{\delta}L_{1}(\delta,\lambda,\tau)=0$ whenever

and $-\infty$ otherwise. Invoking Lemma 3.1, we get

Therefore for any $(\lambda,\tau)$ satisfying (3.10), we have obtained that the dual objective reduces to

Now, recalling the definition of the matrices $A_{i-1,i},D_{i}$ (see (3.5)), we obtain

Finally, using the relations (3.10) to eliminate $\tau_{i}$, and recalling that $\operatorname{val}\!\eqref{P:agpr}$ was defined as $-LR^{2}\operatorname{val}\!\eqref{P:gpr}$, the desired form of the stated dual problem follows. ∎

The next lemma will be crucial in invoking duality in the forthcoming theorem. The proof for this lemma is quite technical and appears in the appendix.

We are now ready to establish a new upper bound on the complexity of the gradient method for values of $h$ between 0 and 1. To the best of our knowledge, the tightest bound thus far is given by (3.1).

First note that both (G) and (G′) are clearly feasible and $\operatorname{val}(G)\leq\operatorname{val}\eqref{P:gpr}$. Invoking Lemma 3.2, by weak duality for the pair of primal-dual problems (G′) and (DG′), we thus obtain that $\operatorname{val}\eqref{P:gpr}\leq\operatorname{val}\eqref{P:dualgpr}$ and hence:

Now consider the following point $(\lambda,t)$ for the dual problem (DG′):

Assuming that this point is (DG′)-feasible, it follows from (3.12) that

which proves the desired result. Thus, it remains to show that the above given choice $(\lambda,t)$ is feasible for (DG′). First, it is easy to see that all the linear constraints of (DG′) on the variables $\lambda_{i},i=1,\ldots,N$ described through the set $\Lambda$ hold true. Now we prove that the matrix $S\equiv S(\lambda,t)$ is positive semidefinite. From Lemma 3.3, with $h\in$, we get that $(1-h)S_{0}+hS_{1}$ is positive definite, as a convex combination of positive definite matrices. Next, we argue that the determinant of $S$ is zero. Indeed, take $u:=(1,\dots,1,-(2Nh+1))^{T}$, then from the definition of $S$ and the choice of $\lambda_{i}$ and $t$ it follows by elementary algebra that $Su=0$. To complete the argument, we note that the determinant of $S$ can also be found via the identity (see, e.g., [7, Section A.5.5]):

Since we have just shown that $(1-h)S_{0}+hS_{1}\succ 0$, then $\det((1-h)S_{0}+hS_{1})>0$ and we get from the above identity that the value of $t-q^{T}((1-h)S_{0}+hS_{1})^{-1}q$, which is the Schur complement of the matrix $S$, is equal to 0. By a well known lemma on the Schur complement [6, Lemma 4.2.1], we conclude that $S$ is positive semidefinite. ∎

The next theorem gives a lower bound on the complexity of Algorithm GM for all values of $h$. In particular, it shows that the bound (3.11) is tight and that it can be attained by a specific convex function in $C^{1,1}_{L}$.

We will describe two functions that attain the two parts of the $\max$ expression in the above claimed statement. For the sake of simplicity we will assume that $L=1$ and $R=\|x_{\ast}-x_{0}\|=1$. Generalizing this proof to general values of $L$ and $R$ can be done by an appropriate scaling.

To show the first part of the $\max$ expression, consider the function

To show the second part of the $\max$ expression, we apply the gradient method on

with $x_{0}=\nu$. We then get that for $i=0,\dots,N$:

Numerical experiments we have performed on problem (G) suggest that the lower complexity bound given by Theorem 3.2 is in fact the exact complexity bound on the gradient method with $0<h<2$.

Suppose the sequence $x_{0},\dots,x_{N}$ is generated by the gradient method GM with $0<h<2$, then

We now turn our attention to the problem of choosing the step size, $h$. Assuming the above conjecture holds true, the optimal step size (i.e., the step size that produces the best complexity bound) for the gradient method with constant step size is given by the unique positive solution to the equation

The solution of this equation approaches $2$ as $N$ grows. Therefore, assuming the conjecture holds true and $N$ is large enough, the complexity of the gradient method with the optimal step size approaches $\frac{LR^{2}}{8N+2}$. This represents a significant improvement, by the factor of 4, upon the best known bound on the gradient method (3.1), and also supports the observation that the gradient method often performs better in practice than in theory.

A Class of First-Order Methods: Numerical Bounds

The interest in the analysis of first-order algorithms of this type is motivated by the fact that it covers some fundamental first-order schemes beyond the gradient method. In particular, to motivate (FO), let us consider the following two algorithms which are of particular interest, and as we shall see below can be seen as special cases of (FO).

We start with the so-called Heavy Ball Method (HBM). For earlier work on this method see Polyak , and for some interesting modern developments and applications, we refer the reader to Attouch et al. and references therein.

Here the step sizes $\alpha$ and $\beta$ are chosen such that $0\leq\beta<1$ and $0<\alpha<2(1+\beta)$, see .

By recursively eliminating the term $x_{i}-x_{i-1}$ in step 2 of HBM, we can rewrite this step as

Therefore, the heavy ball method is a special case of Algorithm FO with the choice

The next algorithm is Nesterov’s celebrated Fast Gradient Method .

A major breakthrough was achieved by Nesterov in , where he proved that the FGM, which requires almost no increase in computational effort when compared to the basic gradient scheme, achieves the improved rate of convergence $O(1/N^{2})$ for function values. More precisely, one hasThe bound given here, which improves the original bound derived in , was recently obtained in Beck-Teboulle .

This fundamental algorithm discovered about 30 years ago by Nesterov has been recently revived and is currently subject of intensive research activities. For some of its extensions and many applications, see e.g., the recent survey paper Beck-Teboulle and references therein.

At first glance, Algorithm FGM seems different than the scheme FO defined above. Here two sequences are defined: the main sequence $x_{0},\dots,x_{N}$ and an auxiliary sequence $y_{1},\dots,y_{N}$. Observing that the gradient of the function is only evaluated on the auxiliary sequence of points $\{y_{i}\}$, we show in the next proposition that FGM fits in the scheme FO through the following algorithm:

The points $y_{1},\dots,y_{N},x_{N}$ generated by Algorithm FGM′ are identical to the respective points generated by Algorithm FGM.

We will show by induction that the sequence $y_{i}$ generated by Algorithm FGM′ is identical to the sequence $y_{i}$ generated by Algorithm FGM, and that the value of $x_{N}$ generated by Algorithm FGM′ is equal to the value of $x_{N}$ generated by Algorithm FGM.

First note that the sequence $t_{i}$ is defined by the two algorithms in the same way. Now let $\{x_{i}$, $y_{i}\}$ be the sequences generated by FGM and denote by $\{y_{i}^{\prime}\}$, $x_{N}^{\prime}$ the sequence generated by FGM′. Obviously, $y^{\prime}_{1}=y_{1}$ and since $t_{1}=1$ we get using the relations 4.2:

Assuming $y_{i}^{\prime}=y_{i}$ for $i=1,\dots,n$, we then have

2 A Numerical Bound for FO

To build the performance estimation problem for Algorithm FO, from which a complexity bound can be derived, we follow the approach used to derive problem (G) for the gradient method. The only difference being that here, of course, the relation between the variables $x_{i}$ is derived from the main iteration of Algorithm FO. After some algebra, the resulting PEP for the class of algorithms (FO) reads

In view of the difficulties in the analysis required to find the solution of (G), an analytical solution to this more general case seems unlikely. However, as we now proceed to show, we can derive a numerical bound on this problem that can be efficiently computed.

Following the analysis of the gradient method, (cf.(G′) in §3.2) we consider the following simpler relaxed problem:

With the same proof as given in Lemma 3.2 for problem (Q′), we obtain that a dual problem for (Q′) is given by the following convex semidefinite optimization problem:

Note that the data matrices of both primal-dual problems (Q′) and (DQ′) depend on the step-sizes $h^{(i)}_{k}$. To avoid a trivial bound for problem $\eqref{P:fopr}$, here we need the following assumption on the dual problem (DQ′):

Assumption 1 Problem (DQ′) is solvable, i.e., the minimum is finite and attained for the given step-sizes $h^{(i)}_{k}$.

Actually, the attainment requirement can be avoided if we can exhibit a feasible point $(\lambda,\tau,t)$ for the problem (DQ′). As noted earlier, given the difficulties already encountered for the simpler gradient method, finding explicitly such a point for the general class of algorithms (FO) is unlikely. However, the structure of problem (DQ′) will be very helpful in the analysis of the next section which further addresses the role of the step-sizes.

The promised numerical bound now easily follows showing that a complexity bound for (FO) is determined by the optimal value of the dual problem (DQ′) which can be computed efficiently by any numerical algorithm for SDP (see e.g., ).

Follows from weak duality for the pair of primal-dual problems (Q′)-(DQ′) ∎

3 Numerical Illustrations

We apply Proposition 4.2 to find bounds on the complexity of the heavy ball method (HBM) withAccording to our simulations, this choice for the values of $\alpha,\beta$ produce results that are typical of the behavior of the algorithm. $\alpha=1$ and $\beta=\frac{1}{2}$ and on the fast gradient method (FGM) with $h_{k}^{(i)}$ as given in (4.2), which as shown earlier, can both be viewed as particular realizations of (FO).

The resulting SDP programs were solved for different values of $N$ using CVX . These results, together with the classical bound on the convergence rate of the main sequence of the fast gradient method (4.1), are summarized in Figures 1 and 2.

Note that as far as the authors are aware, there is no known convergence rate result for the HBM on the class of convex functions in $C^{1,1}_{L}$. As can be seen from the above results, the numerical bound for HBM behaves slightly better than the gradient method (compare with the explicit bound given in Theorem 3.1), but remains much slower than the fast gradient scheme (FGM).

Considering the results on the FGM, note that the numerical bounds for the main sequence of point $x_{i}$ and the corresponding values at the auxiliary sequence $y_{i}$ of the fast gradient method are very similar and perform slightly better than predicted by the classical bound (4.1). To the best of our knowledge, the complexity of the auxiliary sequence is yet unknown, thus these results encourage us to raise the following conjecture.

Let $x_{0},x_{1},\dots$ and $y_{1},y_{2},\dots$ be the main and auxiliary sequences defined by FGM (respectively), then $\{f(x_{i})\}$ and $\{f(y_{i})\}$ converge to the optimal value of the problem with the same rate of convergence.

A Best Performing Algorithm: Optimal Step Sizes for The Algorithm Class FO

We now consider the problem of finding the “best” performing algorithm of the form FO with respect to the new bounds. Namely, we consider the problem of minimizing $\operatorname{val}\eqref{P:fopr}$, the optimal value of (Q′), with respect to the step sizes $h:=(h^{(i)}_{k})_{0\leq k<i\leq N}$ defining the algorithm FO, and which are now considered as unknown variables in FO.

Note that the feasibility of (BIL) follows from the proof of Theorem 3.2, where an explicit feasible point is given to (DG′), which is a special instance of (BIL) when the steps $(h^{(i)}_{k})$ are chosen as in the gradient method.

and denoting $r=(r_{i,k})_{0\leq k<i\leq N}$, we obtain the following linear SDP relaxation of (BIL):

This convex SDP can now be efficiently solved by numerical methods. As the following theorem shows, its solution can be used to construct a solution for (BIL) with optimal step sizes $h$.

Suppose $(r^{\ast},\lambda^{\ast},\tau^{\ast},t^{\ast})$ is an optimal solution for (LIN), then $(h,\lambda^{\ast},\tau^{\ast},t^{\ast})$ is an optimal solution for (BIL), where $h=(h_{k}^{(i)})_{0\leq k<i\leq N}$ is defined by the following recursive rule

As (LIN) is a relaxation of (BIL), it is enough to show that (BIL) can achieve the same objective value. Let $(r^{\ast},\lambda^{\ast},\tau^{\ast},t^{\ast})$ be an optimal solution for (LIN). If $\lambda_{i}^{\ast}\neq 0$ for all $1\leq i\leq N$, then (5.2) satisfies all the equations in (5.1) and therefore $(h,\lambda^{\ast},\tau^{\ast},t^{\ast})$ is feasible for (BIL).

Suppose $\lambda^{\ast}_{m}=0$ for some $1\leq m\leq N$ and that $m$ is the maximal index with this property. Then by the equality and non-negativity constraints in (LIN), we get that $\lambda^{\ast}_{1}=\lambda^{\ast}_{2}=\dots=\lambda^{\ast}_{m}=0$ and $\tau^{\ast}_{0}=\tau^{\ast}_{1}=\dots=\tau^{\ast}_{m-1}=0$. Let $S:=S(r,\lambda^{\ast},\tau^{\ast})$, then by the positive semidefinite constraint in (LIN), we have $S\succeq 0$. From the linear equalities connecting $\lambda$ and $\tau$ it follows that

and we get that $S_{1,1}=\dots=S_{m,m}=0$. By the properties of positive semidefinite matrices we now get that $r_{i,k}^{\ast}=0$ for $i=1,\dots,m$ and $k=0,\dots,i-1$, hence the set of equations (5.1) with the chosen values of $h_{k}^{(i)}$ is consistent. ∎

The optimal value of (LIN) for various values of $N$ is summarized in Figure 3. The resulting new algorithm with the computed optimal step sizes $h_{k}^{(i)}$ is illustrated for $N=5$ and given in Figure 4. As can be seen from these results, (compare with Figure 2) the performance of the new algorithm is almost exactly two times better than the performance of the fast gradient method.

Acknowledgements

This work was initiated during our participation to the “Modern Trends in Optimization and Its Application” program at IPAM, (UCLA), September-December 2010. We would like to thank IPAM for their support and for the very pleasant and stimulating environment provided to us during our stay. We would also like to thank Simi Haber, Ido Ben-Eliezer and Rani Hod for their help in the proof of Lemma 3.3.

Appendix A Proof of Lemma 3.3

We now establish the positive definiteness of the matrices $S_{0}$ and $S_{1}$ given in (3.8) and (3.9), respectively.

We begin by showing that $S_{0}$ is positive definite. Recall that

Let us look at $x^{T}S_{0}x$ for any $x=(x_{0},\dots,x_{N})^{T}$:

which is always positive for $x\neq 0$. We conclude that $S_{0}$ is positive definite.

We will show that $S_{1}$ is positive definite using Sylvester’s criterionDespite the interesting structure of the matrix $S_{1}$, this proof is quite involved. A simpler proof would be most welcome!.

We begin by deriving a recursion rule for the determinant of matrices of the following form:

To find the determinant of $M_{k}$, subtract the one before last row multiplied by $\frac{a_{k}}{a_{k-1}}$ from the last row: the last row becomes

Expanding the determinant along the last row we get

where $(M_{k})_{k,k-1}$ denotes the $k,k-1$ minor:

If we multiply the last column of $(M_{k})_{k,k-1}$ by $\frac{a_{k-1}}{a_{k}}$ we get a matrix that is different from $M_{k-1}$ by only the corner element. Thus by basic determinant properties we get that

Combining these two results, we have found the following recursion rule for $\det M_{k}$, $k\geq 2$:

Obviously, the recursion base cases are given by

Closed form expressions for the determinants

Going back to our matrix, $S_{1}$, by choosing

we get that $M_{k}$ is the $k+1$’th leading principal minor of the matrix $S_{1}$. The recursion rule (A.1) can now be solved for this choice of $a_{i}$ and $d_{i}$. The solution is given by:

Verification

We now proceed to verify the expressions (A.2) and (A.3) given above. We will show that these expressions satisfy the recursion rule (A.1) and the base cases of the problem. We begin by verifying the base cases:

then the recursion rule (A.1) can be written as

Substituting (A.2) in the RHS of (A.1) we get that for $k=2,\dots,N$

It is straightforward (although somewhat involved) to verify that for $k<N$

and thus (A.2) satisfies (A.1). It is also possible to show that

and the expression (A.3) is also verified.

To complete the proof, note that the closed form expressions for $\det M_{k}$ consist of sums and products of positive values, hence $\det M_{k}$ is positive, and thus follows from Sylvester’s criterion that $S_{1}$ is positive definite.

References