Adaptive Restart for Accelerated Gradient Schemes

Introduction

Accelerated gradient schemes were first proposed by Yurii Nesterov in 1983, . He demonstrated a simple modification to gradient descent that could obtain provably optimal performance for the complexity class of first-order algorithms applied to minimize smooth convex functions. The method, and its successors, are often referred to as ‘accelerated methods’. In recent years there has been a resurgence of interest in first-order optimization methods , driven primarily by the need to solve very large problem instances unsuited to second-order methods.

Accelerated gradient schemes can be thought of as momentum methods, in that the step taken at the current iteration depends on the previous iterations, and where the momentum grows from one iteration to the next. When we refer to restarting the algorithm we mean starting the algorithm again, taking the current iteration as the new starting point. This erases the memory of previous iterations and resets the momentum back to zero.

Unlike gradient descent, accelerated methods are not guaranteed to be monotone in the objective value. A common observation when running an accelerated method is the appearance of ripples or bumps in the trace of the objective value; these are seemingly regular increases in the objective, see Figure (1) for an example. In this paper we demonstrate that this behavior occurs when the momentum has exceeded a critical value (the optimal momentum value derived by Nesterov in ) and that the period of these ripples is proportional to the square-root of the (local) condition number of the function. Separately, we show that the optimal restart interval is also proportional to the square root of the condition number. Combining these results we show that restarting when we observe an increase in the function value allows us to recover the optimal linear convergence rate in many cases. Indeed if the function is locally well-conditioned we can use restarting to obtain a linear convergence rate inside the well-conditioned region.

We wish to minimize a smooth convex function of a variable $x\in{\mbox{\bf R}}^{n}$ ,

where $f:{\mbox{\bf R}}^{n}\rightarrow{\mbox{\bf R}}$ has a Lipschitz continuous gradient with constant $L$, i.e.,

We shall denote by $f^{\star}$ the optimal value of the above optimization problem, if the minimizer exists and is unique then we shall write it as $x^{\star}$. Further, a function is said to be strongly convex if there exists a $\mu>0$ such that

where $\mu$ is referred to as the strong convexity parameter. The condition number of a smooth, strongly convex function is $L/\mu$.

Accelerated methods

Accelerated first-order methods to solve (1) were first developed by Nesterov , this scheme is from :

There are many variants of the above scheme, see, ${\it e.g.}$, . Note that by setting $q=1$ in the above scheme we recover gradient descent. For a smooth convex function the above scheme converges for any $t_{k}\leq 1/L$; setting $t_{k}=1/L$ and $q=0$ obtains a guaranteed convergence rate of

If the function is also strongly convex with parameter $\mu$, then a choice of $q=\mu/L$ (the reciprocal of the condition number) will achieve

This is often referred to as linear convergence. With this convergence rate we can achieve an accuracy of $\epsilon$ in

In the case of a strongly convex function the following simpler scheme obtains the same guaranteed rate of convergence :

Note that in Algorithm 1, using the optimal choice $q=\mu/L$, we have that $\beta_{k}\uparrow\beta^{\star}$. Taking $\beta_{k}$ to be a momentum parameter, then for a strongly convex function $\beta^{\star}$ is the maximum amount of momentum we should apply; when we have a value of $\beta$ higher than $\beta^{\star}$ we refer to it as ‘high momentum’. We shall return to this point later.

The convergence of these schemes is optimal in the sense of the lower complexity bounds derived by Nemirovski and Yudin in . However, this convergence is only guaranteed when the function parameters $\mu$ and $L$ are known in advance.

A natural question to ask is how robust are accelerated methods to errors in the estimates of the Lipschitz constant $L$ and strong convexity parameter $\mu$? For the case of an unknown Lipschitz constant we can estimate the optimal step-size by the use of backtracking; see, e.g., . Estimating the strong convexity parameter is much more challenging.

In Nesterov demonstrated a method to bound $\mu$, similar to the backtracking scheme for $L$ described above. His scheme achieves a convergence rate quite a bit slower than Algorithm 1 with a known value of $\mu$. In practice, we often assume or guess that $\mu$ is zero, which corresponds to setting $q=0$ in Algorithm 1. Indeed many discussions of accelerated algorithms do not even include a $q$ term; the original algorithm in did not use a $q$. However, this can dramatically slow down the convergence of the algorithm. Figure 1 shows Algorithm 1 applied to minimize a positive definite quadratic function in $n=200$ dimensions, with optimal choice of $q$ being $q^{\star}=\mu/L=4.1\times 10^{-5}$ (a condition number of about $2.4\times 10^{4}$), and step size $t=1/L$. Each trace is the progress of the algorithm with a different choice of $q$ (hence a different estimate of $\mu$).

We observe that slightly over or underestimating the optimal value of $q$ for the function can have a severe detrimental effect on the rate of convergence of the algorithm. We also note the clear difference in behavior between the cases where we underestimate and where we overestimate $q^{\star}$; in the latter we observe monotonic convergence but in the former we notice the appearance of regular ripples or bumps in the traces.

Interpretation.

The optimal momentum depends on the condition number of the function; specifically, higher momentum is required when the function has a higher condition number. Underestimating the amount of momentum required leads to slower convergence. However we are more often in the other regime, that of overestimated momentum, because generally $q=0$, in which case $\beta_{k}\uparrow 1$; this corresponds to high momentum and rippling behavior, as we see in Figure 1. This can be visually understood in Figure (2), which shows the trajectories of sequences generated by Algorithm 1 minimizing a positive definite quadratic in two dimensions, under $q=q^{\star}$, the optimal choice of $q$, and $q=0$. The high momentum causes the trajectory to overshoot the minimum and oscillate around it. This causes a rippling in the function values along the trajectory. Later we shall demonstrate that the period of these ripples is proportional to the square root of the (local) condition number of the function.

Lastly we mention that the condition number is a global parameter; the sequence generated by an accelerated scheme may enter regions that are locally better conditioned, say, near the optimum. In these cases the choice of $q=q^{\star}$ is appropriate outside of this region, but once we enter it we expect the rippling behavior associated with high momentum to emerge, despite the optimal choice of $q$.

Restarting

For strongly convex functions an alternative to choosing the optimal value of $q$ is to use restarting, . One example of a fixed restart scheme is as follows:

We restart the algorithm every $k$ iterations, taking as our starting point the last point produced by the algorithm, where $k$ is a fixed restart interval. In other words we ‘forget’ all previous iterations and reset the momentum back to zero.

We can obtain an upper bound on the optimal restart interval. If we restart every $k$ iterations we have, at outer iteration $j$, inner loop iteration $k$ (just before a restart),

where the first inequality is the convergence guarantee of Algorithm 1, and the second comes from the strong convexity of $f$. So after $jk$ steps we have

If we assume we have $jk=c$ total iterations and we wish to minimize $(8L/\mu k^{2})^{j}$ over $j$ and $k$ jointly, we obtain

Using this as our restart interval we obtain an accuracy of $\epsilon$ in less than $\mathcal{O}(\sqrt{L/\mu}\log(1/\epsilon))$ iterations, i.e., the optimal linear convergence rate as in equation (4).

The drawbacks in using fixed restarts are that firstly it depends on unknown parameters $L$ and, more importantly, $\mu$, and secondly it is a global parameter that may be inappropriate in better conditioned regions.

2 Adaptive restart

The above analysis suggests that an adaptive restart technique may be useful. In particular we want a scheme that makes some computationally cheap observation and decides whether or not to restart based on that observation. In this paper we suggest two schemes that perform well in practice and provide some analysis to show accelerated convergence when these schemes are used.

Empirically we observe that these two schemes perform similarly well. The gradient scheme has two advantages over the function scheme. Firstly near to the optimum the gradient scheme may be more numerically stable. Secondly all quantities involved in the gradient scheme are already calculated in accelerated schemes, so no extra computation is required.

We can give rough justifications for each scheme. The function scheme restarts at the bottom of the troughs as in Figure 1, thereby avoiding the wasted iterations where we are moving away from the optimum. The gradient scheme restarts whenever the momentum term and the negative gradient are making an obtuse angle. In other words we restart when the momentum seems to be taking us in a bad direction, as measured by the negative gradient at that point.

Figure 3 shows the effect of different restart intervals on minimizing a positive definite quadratic function in $n=500$ dimensions. In this particular case the upper bound on the optimal restart interval is every $700$ iterations. We note that when this interval is used the convergence is better than when no restart is used, however not as good as using the optimal choice of $q$. We also note that restarting every $400$ iterations performs about as well as restarting every $700$ iterations, suggesting that the optimal restart interval is somewhat lower than $700$. We have also plotted the performance of the two adaptive restart schemes. The performance is on the same order as the algorithm with the optimal $q$ and much better than using the fixed restart interval. (Conjugate gradient methods, , will generally outperform an accelerated gradient scheme when minimizing a quadratic; we use quadratics here simply for illustrative purposes.)

Figure 4 demonstrates the function restart scheme trajectories in the two dimensional example, restarting resets the momentum and prevents the characteristic spiralling behavior.

Analysis

In this section we consider applying an accelerated scheme to minimizing a positive definite quadratic. We shall see that once the momentum is larger than a critical value we observe periodicity in the iterates. We use this to prove linear convergence when using adaptive restarting. The analysis presented in this section is similar in spirit to the analysis of the heavy ball method in [18, §3.2].

Consider minimizing a convex quadratic. Without loss of generality we can assume that $f$ has the following form:

2 The algorithm as a linear dynamical system

We shall assume a fixed step-size $t=1/L$ for simplicity. Given quantities $x^{0}$ and $y^{0}=x^{0}$, Algorithm 1 is carried out as follows,

For the rest of the analysis we shall take $\beta_{k}$ to be constant and equal to some $\beta$ for all $k$. This is a somewhat crude approximation, but by making it we can show that there are two regimes of behavior of the system, depending on the value of $\beta$. Consider the eigenvector decomposition of $A=V\Lambda V^{T}$. Denote by $w^{k}=V^{T}x^{k}$, $v^{k}=V^{T}y^{k}$. In this basis the update equations can be written

These are $n$ independently evolving dynamical systems. The $i$th system evolves according to

where $\lambda_{i}$ is the $i$th eigenvalue of $A$. Eliminating the sequence $v_{i}^{(k)}$ from the above we obtain the following recurrence relation for the evolution of $w_{i}$:

where $w_{i}^{0}$ is known and $w_{i}^{1}=w_{i}^{0}(1-\lambda_{i}/L)$, ${\it i.e.}$, a gradient step from $w_{i}^{0}$.

The update equation for $v_{i}$ is identical, differing only in the initial conditions,

where $v_{i}^{0}=w_{i}^{0}$ and $v_{i}^{1}=((1+\beta)(1-\lambda_{i}/L)-\beta)v_{i}^{0}$.

3 Convergence properties

The behavior of this system is determined by the characteristic polynomial of the recurrence relation,

Let $\beta^{\star}_{i}$ be the critical value of $\beta$ for which this polynomial has repeated roots, i.e.,

If $\beta\leq\beta^{\star}_{i}$ then the polynomial (7) has two real roots, $r_{1}$ and $r_{2}$, and the system evolves according to

When $\beta=\beta^{\star}_{i}$ the roots coincide at the point $r^{\star}=(1+\beta)(1-\lambda_{i}/L)/2=(1-\sqrt{\lambda_{i}/L})$; this corresponds to critical damping. We have the fastest monotone convergence at rate $\propto(1-\sqrt{\lambda_{i}/L})^{k}$. Note that if $\lambda_{i}=\mu$ then $\beta^{\star}_{i}$ is the optimal choice of $\beta$ as given by equation (5) and the convergence rate is the optimal rate, as given by equation (3). This is the case because, as we shall see, the smallest eigenvalue will come to dominate the convergence of the entire system.

If $\beta<\beta^{\star}_{i}$ we are in the low momentum regime, and we say the system is over-damped. The convergence rate is dominated by the larger root, which is greater than $r^{\star}$,i.e., the system exhibits slow monotone convergence.

If $\beta>\beta^{\star}_{i}$ then the roots of the polynomial (7) are complex; we are in the high momentum regime and the system is under-damped and exhibits periodicity. In that case the characteristic solution is given by

and $\delta_{i}$ and $c_{i}$ are constants that depend on the initial conditions; in particular for $\beta\approx 1$ we have $\delta_{i}\approx 0$ and we will ignore it. Similarly,

where $\hat{\delta}_{i}$ and $\hat{c}_{i}$ are constants, and again $\hat{\delta}_{i}\approx 0$. For small $\theta$ we know that $\cos^{-1}(\sqrt{1-\theta})\approx\sqrt{\theta}$, and therefore if $\lambda_{i}\ll L$, then

In particular the frequency of oscillation for the mode corresponding to the smallest eigenvalue $\mu$ is approximately given by $\psi_{\mu}\approx\sqrt{\mu/L}$.

To summarize, based on the value of $\beta$ we observe the following behaviors:

$\beta>\beta_{i}^{\star}$: high momentum, under-damped

$\beta<\beta_{i}^{\star}$: low momentum, over-damped

$\beta=\beta_{i}^{\star}$: optimal momentum, critically damped.

4 Observable quantities

We don’t observe the evolution of the modes, but we can observe the evolution of the function value; which is given by

and if $\beta>\beta^{\star}=(1-\sqrt{\mu/L})/(1+\sqrt{\mu/L})$ we are in the high momentum regime for all modes and thus

The function value will quickly be dominated by the smallest eigenvalue and we have that

In other words observing the quantities in (9) or (10) we expect to see oscillations at a frequency proportional to $\sqrt{\mu/L}$, i.e., the frequency of oscillation is telling us something about the condition number of the function.

5 Convergence with adaptive restart

If we apply Algorithm 1 with $q=0$ to minimize a quadratic we start with $\beta_{0}=0$, i.e., the system is in the low momentum, monotonic regime. Eventually $\beta_{k}$ becomes larger than $\beta^{\star}$ and we enter the high momentum, oscillatory regime. It takes about $(3/2)\sqrt{L/\mu}$ iterations for $\beta_{k}$ to exceed $\beta^{\star}$. After that the system is under-damped and the iterates obey equations (9) and (10). Under either adaptive restart scheme, equations (9) and (10) indicate that we shall observe the restart condition after a further $(\pi/2)\sqrt{L/\mu}$ iterations. We restart and the process begins again, with $\beta_{k}$ set back to zero. Thus under either scheme we restart approximately every

iterations (cf., the upper bound on optimal fixed restart interval (6)). Following a similar derivation to §3.1, this restart interval guarantees us an accuracy of $\epsilon$ within $\mathcal{O}(\sqrt{L/\mu}\log(1/\epsilon))$ iterations, i.e., we have recovered the optimal linear convergence rate of equation (4) via adaptive restarting, with no prior knowledge of $\mu$.

6 Extension to smooth convex minimization

In many cases the function we are minimizing is well approximated by a quadratic near the optimum, i.e., there is a region inside of which

The effect of these restart schemes outside of the quadratic region is unclear. In practice we observe that restarting based on one of the criteria described above is almost always helpful, even far away from the optimum. However, we have observed cases where restarting far from the optimum can slow down the early convergence slightly, until the quadratic region is reached and the algorithm enters the rapid linear convergence phase.

Numerical examples

In this section we describe three further numerical examples that demonstrate the improvement of accelerated algorithms under an adaptive restarting technique.

Here we minimize a smooth convex function that is not strongly convex. Consider the following optimization problem

where $x\in{\mbox{\bf R}}^{n}$. The objective function is smooth, but not strongly convex, it grows linearly asymptotically. Thus, the optimal value of $q$ in Algorithm 1 is zero. The quantity $\rho$ controls the smoothness of the function, as $\rho\rightarrow 0$, $f(x)\rightarrow\max_{i=1,\ldots,m}(a_{i}^{T}x-b_{i})$. As it is smooth, we expect the region around the optimum to be well approximated by a quadratic (we consider only examples where the optimal value is finite), and thus we expect to eventually enter a region where our restart method will obtain linear convergence without any knowledge of where this region is, the size of the region or the local function parameters within this region. For smaller values of $\rho$ the smoothness of the objective function decreases and thus we expect to take more iterations before we enter the region of linear convergence.

As a particular example we took $n=20$ and $m=100$; we generated the $a_{i}$ and $b_{i}$ randomly. Figure 5 demonstrates the performance of four different schemes for four different values of $\rho$. We selected the step size for each case using backtracking. We note that both restart schemes perform well, eventually beating both gradient descent and the accelerated scheme. Both the function and gradient schemes eventually enter a region of fast linear convergence. For large $\rho$ we see that even gradient descent performs well, as, similar to the restarted method, it is able to automatically exploit the local strong convexity of the quadratic region around the optimum. Notice also the appearance of the periodic behavior.

2 Sparse linear regression

Consider the following optimization problem:

over $x\in{\mbox{\bf R}}^{n}$, where $A\in{\mbox{\bf R}}^{m\times n}$ and in general $n\gg m$. This is a widely studied problem in the field of compressed sensing, see e.g., . Loosely speaking problem (11) seeks a sparse vector with a small measurement error. The quantity $\rho$ trades off these two competing objectives. The iterative soft-threshold algorithm (ISTA) can be used to solve (11) . ISTA relies on the soft-thresholding operator:

where all the operations are applied elementwise. The ISTA algorithm, with constant step-size $t$, is given by

The convergence rate of ISTA is guaranteed to be at least $\mathcal{O}(1/k)$, making it analogous to gradient descent.

The fast iterative soft thresholding algorithm (FISTA) was developed in ; a similar algorithm was also developed by Nesterov in . FISTA essentially applies acceleration to the ISTA algorithm; it is carried out as follows,

It is easy to show that the function adaptive restart scheme can be performed without an extra application of the matrix $A$, which is the costly operation in the algorithm.

In performing FISTA we do not evaluate a gradient, however FISTA can be thought of as a generalized gradient scheme, in which we take

to be a generalized gradient step, where $G(y^{k})$ is the generalized gradient at $y^{k}$. In this case the gradient restart scheme amounts to restarting whenever

3 Quadratic programming

Consider the following quadratic program,

over $x\in{\mbox{\bf R}}^{n}$, where $Q\in{\mbox{\bf R}}^{n\times n}$ is positive definite and $a,b\in{\mbox{\bf R}}^{n}$ are fixed vectors. The constraint inequalities are to be interpreted element-wise, and we assume that $a<b$. We denote by $\Pi_{\mathcal{C}}(z)$ the projection of a point $z$ onto the constraint set, which amounts to thresholding the entries in $z$.

Projected gradient descent can solve (13); it is carried out as follows,

Projected gradient descent obtains a guaranteed convergence rate of $\mathcal{O}(1/k)$. Acceleration has been successfully applied to the projected gradient method, .

The presence of constraints make this a non-smooth optimization problem, however once the constraints that are active have been identified the problem reduces to minimizing a quadratic on a subset of the variables, and we expect adaptive restarting to increase the rate of convergence. As in the sparse regression example we can use the generalized gradient in our gradient based restart scheme, i.e., we restart based on condition (12).

Summary

In this paper we have demonstrated a simple heuristic adaptive restart technique that can improve the convergence performance of accelerated gradient schemes for smooth convex optimization. We restart the algorithm whenever we observe a certain condition on the objective function value or gradient value. We provided some qualitative analysis to show that we can recover the optimal linear rate of convergence in many cases; in particular near the optimum of a smooth function we can potentially dramatically accelerate the rate of convergence, even if the function is not globally strongly convex. We demonstrated the performance of the scheme on some simple numerical examples.

Acknowledgments

We are very grateful to Stephen Boyd for his help and encouragement. We would also like to thank Stephen Wright for his advice and feedback, and Stephen Becker and Michael Grant for useful discussions. E. C. would like to thank the ONR (grant N00014-09-1-0258) and the Broadcom Foundation for their support.

References