A Differential Equation for Modeling Nesterov's Accelerated Gradient Method: Theory and Insights

Introduction

In many fields of machine learning, minimizing a convex function is at the core of efficient model estimation. In the simplest and most standard form, we are interested in solving

First-order methods have regained popularity as data sets and problems are ever increasing in size and, consequently, there has been much research on the theory and practice of accelerated first-order schemes. Perhaps the earliest first-order method for minimizing a convex function $f$ is the gradient method, which dates back to Euler and Lagrange. Thirty years ago, however, in a seminal paper Nesterov proposed an accelerated gradient method (Nesterov 1983), which may take the following form: starting with $x_{0}$ and $y_{0}=x_{0}$, inductively define

For any fixed step size $s\leq 1/L$, where $L$ is the Lipschitz constant of $\nabla f$, this scheme exhibits the convergence rate

Above, $x^{\star}$ is any minimizer of $f$ and $f^{\star}=f(x^{\star})$. It is well-known that this rate is optimal among all methods having only information about the gradient of $f$ at consecutive iterates (Nesterov 2004). This is in contrast to vanilla gradient descent methods, which have the same computational complexity but can only achieve a rate of $O(1/k)$. This improvement relies on the introduction of the momentum term $x_{k}-x_{k-1}$ as well as the particularly tuned coefficient $(k-1)/(k+2)\approx 1-3/k$. Since the introduction of Nesterov’s scheme, there has been much work on the development of first-order accelerated methods, see Nesterov 2004; Nesterov 2005; Nesterov 2013 for theoretical developments, and Tseng 2008 for a unified analysis of these ideas. Notable applications can be found in sparse linear regression (Beck and Teboulle 2009; Qin and Goldfarb 2012), compressed sensing (Becker et al. 2011) and, deep and recurrent neural networks (Sutskever et al. 2013).

In a different direction, there is a long history relating ordinary differential equation (ODEs) to optimization, see Helmke and Moore 1996, Schropp and Singer 2000, and Fiori 2005 for example. The connection between ODEs and numerical optimization is often established via taking step sizes to be very small so that the trajectory or solution path converges to a curve modeled by an ODE. The conciseness and well-established theory of ODEs provide deeper insights into optimization, which has led to many interesting findings. Notable examples include linear regression via solving differential equations induced by linearized Bregman iteration algorithm (Osher et al. 2014), a continuous-time Nesterov-like algorithm in the context of control design (Dürr and Ebenbauer 2012; Dürr et al. 2012), and modeling design iterative optimization algorithms as nonlinear dynamical systems (Lessard et al. 2014).

In this work, we derive a second-order ODE which is the exact limit of Nesterov’s scheme by taking small step sizes in (1); to the best of our knowledge, this work is the first to use ODEs to model Nesterov’s scheme or its variants in this limit. One surprising fact in connection with this subject is that a first-order scheme is modeled by a second-order ODE. This ODE takes the following form:

Approximate equivalence between Nesterov’s scheme and the ODE is established later in various perspectives, rigorous and intuitive. In the main body of this paper, examples and case studies are provided to demonstrate that the homogeneous and conceptually simpler ODE can serve as a tool for understanding, analyzing and generalizing Nesterov’s scheme.

In the following, two insights of Nesterov’s scheme are highlighted, the first one on oscillations in the trajectories of this scheme, and the second on the peculiar constant 3 appearing in the ODE.

In general, Nesterov’s scheme is not monotone in the objective function value due to the introduction of the momentum term. Oscillations or overshoots along the trajectory of iterates approaching the minimizer are often observed when running Nesterov’s scheme. Figure 1 presents typical phenomena of this kind, where a two-dimensional convex function is minimized by Nesterov’s scheme. Viewing the ODE as a damping system, we obtain interpretations as follows.

Small $t$. In the beginning, the damping ratio $3/t$ is large. This leads the ODE to be an overdamped system, returning to the equilibrium without oscillating; Large $t$. As $t$ increases, the ODE with a small $3/t$ behaves like an underdamped system, oscillating with the amplitude gradually decreasing to zero.

As depicted in Figure 1(a), in the beginning the ODE curve moves smoothly towards the origin, the minimizer $x^{\star}$. The second interpretation “Large $t$’’ provides partial explanation for the oscillations observed in Nesterov’s scheme at later stage. Although our analysis extends farther, it is similar in spirit to that carried in O’Donoghue and Candès 2013. In particular, the zoomed Figure 1(b) presents some butterfly-like oscillations for both the scheme and ODE. There, we see that the trajectory constantly moves away from the origin and returns back later. Each overshoot in Figure 1(b) causes a bump in the function values, as shown in Figure 1(c). We observe also from Figure 1(c) that the periodicity captured by the bumps are very close to that of the ODE solution. In passing, it is worth mentioning that the solution to the ODE in this case can be expressed via Bessel functions, hence enabling quantitative characterizations of these overshoots and bumps, which are given in full detail in Section 3.

2 A Phase Transition

The constant 3, derived from $(k+2)-(k-1)$ in (3), is not haphazard. In fact, it is the smallest constant that guarantees $O(1/t^{2})$ convergence rate. Specifically, parameterized by a constant $r$, the generalized ODE

can be translated into a generalized Nesterov’s scheme that is the same as the original (1) except for $(k-1)/(k+2)$ being replaced by $(k-1)/(k+r-1)$. Surprisingly, for both generalized ODEs and schemes, the inverse quadratic convergence is guaranteed if and only if $r\geq 3$. This phase transition suggests there might be deep causes for acceleration among first-order methods. In particular, for $r\geq 3$, the worst case constant in this inverse quadratic convergence rate is minimized at $r=3$.

Figure 2 illustrates the growth of $t^{2}(f(X(t))-f^{\star})$ and $sk^{2}(f(x_{k})-f^{\star})$, respectively, for the generalized ODE and scheme with $r=1$, where the objective function is simply $f(x)=\frac{1}{2}x^{2}$. Inverse quadratic convergence fails to be observed in both Figures 2(a) and 2(b), where the scaled errors grow with $t$ or iterations, for both the generalized ODE and scheme.

3 Outline and Notation

The rest of the paper is organized as follows. In Section 2, the ODE is rigorously derived from Nesterov’s scheme, and a generalization to composite optimization, where $f$ may be non-smooth, is also obtained. Connections between the ODE and the scheme, in terms of trajectory behaviors and convergence rates, are summarized in Section 3. In Section 4, we discuss the effect of replacing the constant $3$ in (3) by an arbitrary constant on the convergence rate. A new restarting scheme is suggested in Section 5, with linear convergence rate established and empirically observed.

Derivation

First, we sketch an informal derivation of the ODE (3). Assume $f\in\mathcal{F}_{L}$ for $L>0$. Combining the two equations of (1) and applying a rescaling gives

Introduce the Ansatz $x_{k}\approx X(k\sqrt{s})$ for some smooth curve $X(t)$ defined for $t\geq 0$. Put $k=t/\sqrt{s}$. Then as the step size $s$ goes to zero, $X(t)\approx x_{t/\sqrt{s}}=x_{k}$ and $X(t+\sqrt{s})\approx x_{(t+\sqrt{s})/\sqrt{s}}=x_{k+1}$, and Taylor expansion gives

and $\sqrt{s}\nabla f(y_{k})=\sqrt{s}\nabla f(X(t))+o(\sqrt{s})$. Thus (4) can be written as

By comparing the coefficients of $\sqrt{s}$ in (5), we obtain

The first initial condition is $X(0)=x_{0}$. Taking $k=1$ in (4) yields

Hence, the second initial condition is simply $\dot{X}(0)=0$ (vanishing initial velocity).

One popular alternative momentum coefficient is $\theta_{k}(\theta_{k-1}^{-1}-1)$, where $\theta_{k}$ are iteratively defined as $\theta_{k+1}=\left(\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2}\right)/2$, starting from $\theta_{0}=1$ (Nesterov 1983; Beck and Teboulle 2009). Simple analysis reveals that $\theta_{k}(\theta_{k-1}^{-1}-1)$ asymptotically equals $1-3/k+O(1/k^{2})$, thus leading to the same ODE as (1).

The next theorem, in a rigorous way, guarantees the validity of the derivation of this ODE. The proofs of both theorems are deferred to the appendices.

For any $f\in\mathcal{F}_{\infty}$, as the step size $s\rightarrow 0$, Nesterov’s scheme (1) converges to the ODE (3) in the sense that for all fixed $T>0$,

Also note that minimizing $f/c^{2}$ is equivalent to minimizing $f$. Hence, the ODE is invariant under the time change. In fact, it is easy to see that time invariance holds if and only if the coefficient of $\dot{X}$ has the form $C/t$ for some constant $C$. Rotational Invariance. Nesterov’s scheme and other gradient-based schemes are invariant under rotations. As expected, the ODE is also invariant under orthogonal transformation. To see this, let $Y=QX$ for some orthogonal matrix $Q$. This leads to $\dot{Y}=Q\dot{X},\ddot{Y}=Q\ddot{X}$ and $\nabla_{Y}f=Q\nabla_{X}f$. Hence, denoting by $Q^{T}$ the transpose of $Q$, the ODE in the new coordinate system reads $Q^{T}\ddot{Y}+\frac{3}{t}Q^{T}\dot{Y}+Q^{T}\nabla_{Y}f=0$, which is of the same form as (3) once multiplying $Q$ on both sides. Initial Asymptotic. Assume sufficient smoothness of $X$ such that $\lim_{t\rightarrow 0}\ddot{X}(t)$ exists. The mean value theorem guarantees the existence of some $\xi\in(0,t)$ that satisfies $\dot{X}(t)/t=(\dot{X}(t)-\dot{X}(0))/t=\ddot{X}(\xi)$. Hence, from the ODE we deduce $\ddot{X}(t)+3\ddot{X}(\xi)+\nabla f(X(t))=0$. Taking the limit $t\rightarrow 0$ gives $\ddot{X}(0)=-\nabla f(x_{0})/4$. Hence, for small $t$ we have the asymptotic form:

This asymptotic expansion is consistent with the empirical observation that Nesterov’s scheme moves slowly in the beginning.

2 ODE for Composite Optimization

Connections and Interpretations

In this section, we explore the approximate equivalence between the ODE and Nesterov’s scheme, and provide evidence that the ODE can serve as an amenable tool for interpreting and analyzing Nesterov’s scheme. The first subsection exhibits inverse quadratic convergence rate for the ODE solution, the next two address the oscillation phenomenon discussed in Section 1.1, and the last subsection is devoted to comparing Nesterov’s scheme with gradient descent from a numerical perspective.

The original result from Nesterov 1983 states that, for any $f\in\mathcal{F}_{L}$, the sequence $\{x_{k}\}$ given by (1) with step size $s\leq 1/L$ satisfies

Our next result indicates that the trajectory of (3) closely resembles the sequence $\{x_{k}\}$ in terms of the convergence rate to a minimizer $x^{\star}$. Compared with the discrete case, this proof is shorter and simpler.

For any $f\in\mathcal{F}_{\infty}$, let $X(t)$ be the unique global solution to (3) with initial conditions $X(0)=x_{0},\dot{X}(0)=0$. Then, for any $t>0$,

Substituting $3\dot{X}/2+t\ddot{X}/2$ with $-t\nabla f(X)/2$, the above equation gives

where the inequality follows from the convexity of $f$. Hence by monotonicity of $\mathcal{E}$ and non-negativity of $2\|X+t\dot{X}/2-x^{\star}\|^{2}$, the gap satisfies

Making use of the approximation $t\approx k\sqrt{s}$, we observe that the convergence rate in (6) is essentially a discrete version of that in (7), providing yet another piece of evidence for the approximate equivalence between the ODE and the scheme.

2 Quadratic ff and Bessel Functions

For quadratic $f$, the ODE (3) admits a solution in closed form. This closed form solution turns out to be very useful in understanding the issues raised in the introduction.

with $X_{i}(0)=x_{0,i},\dot{X}_{i}(0)=0$. Introduce $Y_{i}(u)=uX_{i}(u/\sqrt{\lambda_{i}})$, which satisfies

This is Bessel’s differential equation of order one. Since $Y_{i}$ vanishes at $u=0$, we see that $Y_{i}$ is a constant multiple of $J_{1}$, the Bessel function of the first kind of order one. Up to a constant multiplier, $J_{1}$ is the unique solution to the Bessel’s differential equation $u^{2}\ddot{J}_{1}+u\dot{J}_{1}+(u^{2}-1)J_{1}=0$ that is finite at the origin. In the analytic expansion of $J_{1}$, $m!!$ denotes the double factorial defined as $m!!=m\times(m-2)\times\cdots\times 2$ for even $m$, or $m!!=m\times(m-2)\times\cdots\times 1$ for odd $m$. It has an analytic expansion:

when $u\rightarrow 0$. Requiring $X_{i}(0)=x_{0,i}$, hence, we obtain

For large $t$, the Bessel function has the following asymptotic form (Watson 1995, see e.g.):

This asymptotic expansion yields (note that $f^{\star}=0$)

On the other hand, (9) and (10) give a lower bound:

In view of (10), Nesterov’s scheme might possibly exhibit $O(1/k^{3})$ convergence rate for strongly convex functions. This convergence rate is consistent with the second inequality in Theorem 6. In Section 4.3, we prove the $O(1/t^{3})$ rate for a generalized version of (3). However, (11) rules out the possibility of a higher order convergence rate.

Recall that the function considered in Figure 1 is $f(x)=0.02x_{1}^{2}+0.005x_{2}^{2}$, starting from $x_{0}=(1,~1)$. As the step size $s$ becomes smaller, the trajectory of Nesterov’s scheme converges to the solid curve represented via the Bessel function. While approaching the minimizer $x^{\star}$, each trajectory displays the oscillation pattern, as well-captured by the zoomed Figure 1(b). This prevents Nesterov’s scheme from achieving better convergence rate. The representation (8) offers excellent explanation as follows. Denote by $T_{1},T_{2}$, respectively, the approximate periodicities of the first component $|X_{1}|$ in absolute value and the second $|X_{2}|$. By (9), we get $T_{1}=\pi/\sqrt{\lambda_{1}}=5\pi$ and $T_{2}=\pi/\sqrt{\lambda_{2}}=10\pi$. Hence, as the amplitude gradually decreases to zero, the function $f=2x_{0,1}^{2}J_{1}(\sqrt{\lambda_{1}}t)^{2}/t^{2}+2x_{0,2}^{2}J_{1}(\sqrt{\lambda_{2}}t)^{2}/t^{2}$ has a major cycle of $10\pi$, the least common multiple of $T_{1}$ and $T_{2}$. A careful look at Figure 1(c) reveals that within each major bump, roughly, there are $10\pi/T_{1}=2$ minor peaks.

3 Fluctuations of Strongly Convex ff

The analysis carried out in the previous subsection only applies to convex quadratic functions. In this subsection, we extend the discussion to one-dimensional strongly convex functions. The Sturm-Picone theory (Hinton 2005, see e.g.) is extensively used all along the analysis.

which, apparently, admits a solution $Y(t)=X(t)$. To apply the Sturm-Picone comparison theorem, consider

To obtain a similar result in the opposite direction, consider

Denote by $0<t_{1}<t_{2}<\cdots$ all the roots of $X(t)-x^{\star}$. Then these roots satisfy, for all $i\geq 1$,

4 Nesterov’s Scheme Compared with Gradient Descent

The ansatz $t\approx k\sqrt{s}$ in relating the ODE and Nesterov’s scheme is formally confirmed in Theorem 2. Consequently, for any constant $t_{c}>0$, this implies that $x_{k}$ does not change much for a range of step sizes $s$ if $k\approx t_{c}/\sqrt{s}$. To empirically support this claim, we present an example in Figure 3(a), where the scheme minimizes $f(x)=\|y-Ax\|^{2}/2+\|x\|_{1}$ with $y=(4,~2,~0)$ and $A(:,1)=(0,~2,~4),~A(:,2)=(1,~1,~1)$ starting from $x_{0}=(2,~0)$ (here $A(:,j)$ is the $j$th column of $A$). From this figure, we are delight to observe that $x_{k}$ with the same $t_{c}$ are very close to each other.

This interesting square-root scaling has the potential to shed light on the superiority of Nesterov’s scheme over gradient descent. Roughly speaking, each iteration in Nesterov’s scheme amounts to traveling $\sqrt{s}$ in time along the integral curve of (3), whereas it is known that the simple gradient descent $x_{k+1}=x_{k}-s\nabla f(x_{k})$ moves $s$ along the integral curve of $\dot{X}+\nabla f(X)=0$. We expect that for small $s$ Nesterov’s scheme moves more in each iteration since $\sqrt{s}$ is much larger than $s$. Figure 3(b) illustrates and supports this claim, where the function minimized is $f=|x_{1}|^{3}+5|x_{2}|^{3}+0.001(x_{1}+x_{2})^{2}$ with step size $s=0.05$ (The coordinates are appropriately rotated to allow $x_{0}$ and $x^{\star}$ lie on the same horizontal line). The circles are the iterates for $k=1,10,20,30,45,60,90,120,150,190,250,300$. For Nesterov’s scheme, the seventh circle has already passed $t=15$, while for gradient descent the last point has merely arrived at $t=15$.

A second look at Figure 3(b) suggests that Nesterov’s scheme allows a large deviation from its limit curve, as compared with gradient descent. This raises the question of the stable step size allowed for numerically solving the ODE (3) in the presence of accumulated errors. The finite difference approximation by the forward Euler method is

Assuming $f$ is sufficiently smooth, we have $\nabla f(x+\delta x)\approx\nabla f(x)+\nabla^{2}f(x)\delta x$ for small perturbations $\delta x$, where $\nabla^{2}f(x)$ is the Hessian of $f$ evaluated at $x$. Identifying $k=t/\Delta t$, the characteristic equation of this finite difference scheme is approximately

The numerical stability of (14) with respect to accumulated errors is equivalent to this: all the roots of (16) lie in the unit circle (Leader 2004, see e.g.). When $\nabla^{2}f\preceq LI_{n}$ (i.e. $LI_{n}-\nabla^{2}f$ is positive semidefinite), if $\Delta t/t$ small and $\Delta t<2/\sqrt{L}$, we see that all the roots of (16) lie in the unit circle. On the other hand, if $\Delta t>2/\sqrt{L}$, (16) can possibly have a root $\lambda$ outside the unit circle, causing numerical instability. Under our identification $s=\Delta t^{2}$, a step size of $s=1/L$ in Nesterov’s scheme (1) is approximately equivalent to a step size of $\Delta t=1/\sqrt{L}$ in the forward Euler method, which is stable for numerically integrating (14).

As a comparison, note that the finite difference scheme of the ODE $\dot{X}(t)+\nabla f(X(t))=0$, which models gradient descent with updates $x_{k+1}=x_{k}-s\nabla f(x_{k})$, has the characteristic equation $\det(\lambda-(1-\Delta t\nabla^{2}f))=0$. Thus, to guarantee $-I_{n}\preceq 1-\Delta t\nabla^{2}f\preceq I_{n}$ in worst case analysis, one can only choose $\Delta t\leq 2/L$ for a fixed step size, which is much smaller than the step size $2/\sqrt{L}$ for (14) when $\nabla f$ is very variable, i.e., $L$ is large.

The Magic Constant 3

Recall that the constant 3 appearing in the coefficient of $\dot{X}$ in (3) originates from $(k+2)-(k-1)=3$. This number leads to the momentum coefficient in (1) taking the form $(k-1)/(k+2)=1-3/k+O(1/k^{2})$. In this section, we demonstrate that 3 can be replaced by any larger number, while maintaining the $O(1/k^{2})$ convergence rate. To begin with, let us consider the following ODE parameterized by a constant $r$:

with initial conditions $X(0)=x_{0},\dot{X}(0)=0$. The proof of Theorem 1, which seamlessly applies here, guarantees the existence and uniqueness of the solution $X$ to this ODE.

Interpreting the damping ratio $r/t$ as a measure of friction In physics and engineering, damping may be modeled as a force proportional to velocity but opposite in direction, i.e. resisting motion; for instance, this force may be used as an approximation to the friction caused by drag. In our model, this force would be proportional to $-\frac{r}{t}\dot{X}$ where $\dot{X}$ is velocity and $\frac{r}{t}$ is the damping coefficient. in the damping system, our results say that more friction does not end the $O(1/t^{2})$ and $O(1/k^{2})$ convergence rate. On the other hand, in the lower friction setting, where $r$ is smaller than 3, we can no longer expect inverse quadratic convergence rate, unless some additional structures of $f$ are imposed. We believe that this striking phase transition at 3 deserves more attention as an interesting research challenge.

Here, we study the convergence rate of (17) with $r>3$ and $f\in\mathcal{F}_{\infty}$. Compared with (3), this new ODE as a damping suffers from higher friction. Following the strategy adopted in the proof of Theorem 3, we consider a new energy functional defined as

By studying the derivative of this functional, we get the following result.

Noting $r\dot{X}+t\ddot{X}=-t\nabla f(X)$, we get $\dot{\mathcal{E}}$ equal to

where the inequality follows from the convexity of $f$. Since $f(X)\geq f^{\star}$, the last display implies that $\mathcal{E}$ is non-increasing. Hence

yielding the first inequality of this theorem. To complete the proof, from (18) it follows that

as desired for establishing the second inequality. ∎

The first inequality is the same as (7) for the ODE (3), except for a larger constant $(r-1)^{2}/2$. The second inequality measures the error $f(X(t))-f^{\star}$ in an average sense, and cannot be deduced from the first inequality.

Parametrizing by a constant $r$, we propose the generalized Nesterov’s scheme,

starting from $y_{0}=x_{0}$. The discrete analog of Theorem 5 is below.

The sequence $\{x_{k}\}$ given by (19) with $0<s\leq 1/L$ satisfies

The first inequality suggests that the generalized Nesterov’s schemes still achieve $O(1/k^{2})$ convergence rate. However, if the error bound satisfies $f(x_{k^{\prime}})-f^{\star}\geq c/k^{\prime 2}$ for some arbitrarily small $c>0$ and a dense subsequence $\{k^{\prime}\}$, i.e., $|\{k^{\prime}\}\cap\{1,\ldots,m\}|\geq\alpha m$ for all $m\geq 1$ and some $\alpha>0$, then the second inequality of the theorem would be violated. To see this, note that if it were the case, we would have $(k^{\prime}+r-1)(f(x_{k^{\prime}})-f^{\star})\gtrsim\frac{1}{k^{\prime}}$; the sum of the harmonic series $\frac{1}{k^{\prime}}$ over a dense subset of $\{1,2,\ldots\}$ is infinite. Hence, the second inequality is not trivial because it implies the error bound is, in some sense, $O(1/k^{2})$ suboptimal.

Now we turn to the proof of this theorem. It is worth pointing out that, though based on the same idea, the proof below is much more complicated than that of Theorem 5.

where $z_{k}=(k+r-1)y_{k}/(r-1)-kx_{k}/(r-1)$. If we have

then it would immediately yield the desired results by summing (20) over $k$. That is, by recursively applying (20), we see

Noting that the left-hand side of (22) is lower bounded by $2s(k+r-2)^{2}(f(x_{k})-f^{\star})/(r-1)$, we thus obtain the first inequality of the theorem. Since $\mathcal{E}(k)\geq 0$, the second inequality is verified via taking the limit $k\rightarrow\infty$ in (22) and replacing $(r-3)(i+r-1)+1$ by $(r-3)(i+r-1)$.

We now establish (20). For $s\leq 1/L$, we have the basic inequality,

for any $x$ and $y$. Note that $y_{k-1}-sG_{s}(y_{k-1})$ actually coincides with $x_{k}$. Summing of $(k-1)/(k+r-2)\times\eqref{eq:prox_ineq}$ with $x=x_{k-1},y=y_{k-1}$ and $(r-1)/(k+r-2)\times\eqref{eq:prox_ineq}$ with $x=x^{\star},y=y_{k-1}$ gives

where we use $z_{k-1}-s(k+r-2)G_{s}(y_{k-1})/(r-1)=z_{k}$. Rearranging the above inequality and multiplying by $2s(k+r-2)^{2}/(r-1)$ gives the desired (20).

In closing, we would like to point out this new scheme is equivalent to setting $\theta_{k}=(r-1)/(k+r-1)$ and letting $\theta_{k}(\theta_{k-1}^{-1}-1)$ replace the momentum coefficient $(k-1)/(k+r-1)$. Then, the equal sign $``="$ in the update $\theta_{k+1}=(\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2})/2$ has to be replaced by an inequality sign $``\geq"$. In examining the proof of Theorem 1(b) in Tseng 2010, we can get an alternative proof of Theorem 6.

2 Low Friction

Now we turn to the case $r<3$. Then, unfortunately, the energy functional approach for proving Theorem 5 is no longer valid, since the left-hand side of (18) is positive in general. In fact, there are counterexamples that fail the desired $O(1/t^{2})$ or $O(1/k^{2})$ convergence rate. We present such examples in continuous time. Equally, these examples would also violate the $O(1/k^{2})$ convergence rate in the discrete schemes, and we forego the details.

Let $f(x)=\frac{1}{2}\|x\|^{2}$ and $X$ be the solution to (17). Then, $Y=t^{\frac{r-1}{2}}X$ satisfies

With the initial condition $Y(t)\approx t^{\frac{r-1}{2}}x_{0}$ for small $t$, the solution to the above Bessel equation in a vector form of order $(r-1)/2$ is $Y(t)=2^{\frac{r-1}{2}}\Gamma((r+1)/2)J_{(r-1)/2}(t)x_{0}$. Thus,

For large $t$, the Bessel function $J_{(r-1)/2}(t)=\sqrt{2/(\pi t)}\big(\cos(t-(r-1)\pi/4-\pi/4)+O(1/t)\big)$. Hence,

where the exponent $r$ is tight. This rules out the possibility of inverse quadratic convergence of the generalized ODE and scheme for all $f\in\mathcal{F}_{L}$ if $r<2$. An example with $r=1$ is plotted in Figure 2.

Next, we consider the case $2\leq r<3$ and let $f(x)=|x|$ (this also applies to multivariate $f=\|x\|$). This function does not have a Lipschitz continuous gradient. However, a similar pattern as in Figure 2 can be also observed if we smooth $|x|$ at an arbitrarily small vicinity of 0. Starting from $x_{0}>0$, we get $X(t)=x_{0}-\frac{t^{2}}{2(1+r)}$ for $t\leq\sqrt{2(1+r)x_{0}}$. Requiring continuity of $X$ and $\dot{X}$ at the change point 0, we get

for $\sqrt{2(1+r)x_{0}}<t\leq\sqrt{2c^{\star}(1+r)x_{0}}$, where $c^{\star}$ is the positive root other than 1 of $(r-1)c+4c^{-\frac{r-1}{2}}=r+3$. Repeating this process solves for $X$. Note that $t^{1-r}$ is in the null space of $\ddot{X}+r\dot{X}/t$ and satisfies $t^{2}\times t^{1-r}\rightarrow\infty$ as $t\rightarrow\infty$. For illustration, Figure 4 plots $t^{2}(f(X(t))-f^{\star})$ and $sk^{2}(f(x_{k})-f^{\star})$ with $r=2,2.5$, and $r=4$ for comparison For Figures 4(d), 4(e) and 4(f), if running generalized Nesterov’s schemes with too many iterations (e.g. $10^{5}$), the deviations from the ODE will grow. Taking a sufficiently small $s$ can solve this issue.. It is clearly that inverse quadratic convergence does not hold for $r=2,2.5$, that is, (2) does not hold for $r<3$. Interestingly, in Figures 4(a) and 4(d), the scaled errors at peaks grow linearly, whereas for $r=2.5$, the growth rate, though positive as well, seems sublinear.

However, if $f$ possesses some additional property, inverse quadratic convergence is still guaranteed, as stated below. In that theorem, $f$ is assumed to be a continuously differentiable convex function.

Suppose $1<r<3$ and let $X$ be a solution to the ODE (17). If $(f-f^{\star})^{\frac{r-1}{2}}$ is also convex, then

Since $(f-f^{\star})^{\frac{r-1}{2}}$ is convex, we obtain

which can be simplified to $\frac{2}{r-1}(f(X)-f^{\star})\leq\langle X-x^{\star},\nabla f(X)\rangle$. This inequality combined with (18) leads to the monotonically decreasing of $\mathcal{E}(t)$ defined for Theorem 5. This completes the proof by noting $f(X)-f^{\star}\leq(r-1)\mathcal{E}(t)/(2t^{2})\leq(r-1)\mathcal{E}(0)/(2t^{2})=(r-1)^{2}\|x_{0}-x^{\star}\|^{2}/(2t^{2})$. ∎

3 Strongly Convex ff

Strong convexity is a desirable property for optimization. Making use of this property carefully suggests a generalized Nesterov’s scheme that achieves optimal linear convergence (Nesterov 2004). In that case, even vanilla gradient descent has a linear convergence rate. Unfortunately, the example given in the previous subsection simply rules out such possibility for (1) and its generalizations (19). However, from a different perspective, this example suggests that $O(t^{-r})$ convergence rate can be expected for (17). In the next theorem, we prove a slightly weaker statement of this kind, that is, a provable $O(t^{-\frac{2r}{3}})$ convergence rate is established for strongly convex functions. Bridging this gap may require new tools and more careful analysis.

When clear from the context, $\mathcal{E}(t;\alpha)$ is simply denoted as $\mathcal{E}(t)$. For $r>3$, taking $\alpha=2r/3$ in the theorem stated below gives $f(X(t))-f^{\star}\lesssim\|x_{0}-x^{\star}\|^{2}/t^{\frac{2r}{3}}$.

for any $t>0$. Above, the constant $C$ only depends on $\alpha$ and $r$.

Note that $\dot{\mathcal{E}}(t;\alpha)$ equals

By the strong convexity of $f$, the second term of the right-hand side of (24) is bounded below as

Substituting the last display into (24) with the awareness of $r\geq 3\alpha/2$ yields

Hence, if $t\geq t_{\alpha}:=\sqrt{(\alpha-2)(2r-\alpha)/(2\mu)}$, we obtain

Integrating the last inequality on the interval $(t_{\alpha},t)$ gives

Making use of (26), we apply induction on $\alpha$ to finish the proof. First, consider $2<\alpha\leq 4$. Applying Theorem 5, from (26) we get that $\mathcal{E}(t)$ is upper bounded by

Then, we bound $\mathcal{E}(t_{\alpha})$ as follows.

where in the second inequality we use the decreasing property of the energy functional defined for Theorem 5. Combining (27) and (28), we have

For $t\geq t_{\alpha}$, it suffices to apply $f(X(t))-f^{\star}\leq\mathcal{E}(t)/t^{3}$ to the last display. For $t<t_{\alpha}$, by Theorem 5, $f(X(t))-f^{\star}$ is upper bounded by

which verify the induction for $t\geq t_{\alpha}$. As for $t<t_{\alpha}$, the validity of the induction follows from Theorem 5, similarly to (29). Thus, combining the base and induction steps, the proof is completed. ∎

It should be pointed out that the constant $C$ in the statement of Theorem 8 grows with the parameter $r$. Hence, simply increasing $r$ does not guarantee to give a better error bound. While it is desirable to expect a discrete analogy of Theorem 8, i.e., $O(1/k^{\alpha})$ convergence rate for (19), a complete proof can be notoriously complicated. That said, we mimic the proof of Theorem 8 for $\alpha=3$ and succeed in obtaining a $O(1/k^{3})$ convergence rate for the generalized Nesterov’s schemes, as summarized in the theorem below.

Suppose $f$ is written as $f=g+h$, where $g\in\mathcal{S}_{\mu,L}$ and $h$ is convex with possible extended value $\infty$. Then, the generalized Nesterov’s scheme (19) with $r\geq 9/2$ and $s=1/L$ satisfies

This theorem states that the discrete scheme (19) enjoys the error bound $O(1/k^{3})$ without any knowledge of the condition number $L/\mu$. In particular, this bound is much better than that given in Theorem 6 if $k\gg\sqrt{L/\mu}$. The strategy of the proof is fully inspired by that of Theorem 8, though it is much more complicated and thus deferred to the Appendix. The relevant energy functional $\mathcal{E}(k)$ for this Theorem 9 is equal to

4 Numerical Examples

We study six synthetic examples to compare (19) with the step sizes are fixed to be $1/L$, as illustrated in Figure 5. The error rates exhibits similar patterns for all $r$, namely, decreasing while suffering from local bumps. A smaller $r$ introduces less friction, thus allowing $x_{k}$ moves towards $x^{\star}$ faster in the beginning. However, when sufficiently close to $x^{\star}$, more friction is preferred in order to reduce overshoot. This point of view explains what we observe in these examples. That is, across these six examples, (19) with a smaller $r$ performs slightly better in the beginning, but a larger $r$ has advantage when $k$ is large. It is an interesting question how to choose a good $r$ for different problems in practice.

Lasso with fat design. Minimize $f(x)=\frac{1}{2}\|Ax-b\|^{2}+\lambda\|x\|_{1}$, in which $A$ a $100\times 500$ random matrix with i.i.d. standard Gaussian $\mathcal{N}(0,1)$ entries, $b$ generated independently has i.i.d. $\mathcal{N}(0,25)$ entries, and the penalty $\lambda=4$. The plot is Figure 5(a).

Lasso with square design. Minimize $f(x)=\frac{1}{2}\|Ax-b\|^{2}+\lambda\|x\|_{1}$, where $A$ a $500\times 500$ random matrix with i.i.d. standard Gaussian entries, $b$ generated independently has i.i.d. $\mathcal{N}(0,9)$ entries, and the penalty $\lambda=4$. The plot is Figure 5(b).

Nonnegative least squares (NLS) with fat design. Minimize $f(x)=\|Ax-b\|^{2}$ subject to $x\succeq 0$, with the same design $A$ and $b$ as in Figure 5(a). The plot is Figure 5(c).

Nonnegative least squares with sparse design. Minimize $f(x)=\|Ax-b\|^{2}$ subject to $x\succeq 0$, in which $A$ is a $1000\times 10000$ sparse matrix with nonzero probability $10\%$ for each entry and $b$ is given as $b=Ax^{0}+\mathcal{N}(0,I_{1000})$. The nonzero entries of $A$ are independently Gaussian distributed before column normalization, and $x^{0}$ has 100 nonzero entries that are all equal to 4. The plot is Figure 5(d).

Restarting

The example discussed in Section 4.2 demonstrates that Nesterov’s scheme and its generalizations (19) are not capable of fully exploiting strong convexity. That is, this example suggests evidence that $O(1/\mathtt{poly}(k))$ is the best rate achievable under strong convexity. In contrast, the vanilla gradient method achieves linear convergence $O((1-\mu/L)^{k})$. This drawback results from too much momentum introduced when the objective function is strongly convex. The derivative of a strongly convex function is generally more reliable than that of non-strongly convex functions. In the language of ODEs, at later stage a too small $3/t$ in (3) leads to a lack of friction, resulting in unnecessary overshoot along the trajectory. Incorporating the optimal momentum coefficient $\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}}$ (This is less than $(k-1)/(k+2)$ when $k$ is large), Nesterov’s scheme has convergence rate of $O((1-\sqrt{\mu/L})^{k})$ (Nesterov 2004), which, however, requires knowledge of the condition number $\mu/L$. While it is relatively easy to bound the Lipschitz constant $L$ by the use of backtracking, estimating the strong convexity parameter $\mu$, if not impossible, is very challenging.

Among many approaches to gain acceleration via adaptively estimating $\mu/L$ (Nesterov 2013, see), O’Donoghue and Candès 2013 proposes a procedure termed as gradient restarting for Nesterov’s scheme in which (1) is restarted with $x_{0}=y_{0}:=x_{k}$ whenever $f(x_{k+1})>f(x_{k})$. In the language of ODEs, this restarting essentially keeps $\langle\nabla f,\dot{X}\rangle$ negative, and resets $3/t$ each time to prevent this coefficient from steadily decreasing along the trajectory. Although it has been empirically observed that this method significantly boosts convergence, there is no general theory characterizing the convergence rate.

In this section, we propose a new restarting scheme we call the speed restarting scheme. The underlying motivation is to maintain a relatively high velocity $\dot{X}$ along the trajectory, similar in spirit to the gradient restarting. Specifically, our main result, Theorem 10, ensures linear convergence of the continuous version of the speed restarting. More generally, our contribution here is merely to provide a framework for analyzing restarting schemes rather than competing with other schemes; it is beyond the scope of this paper to get optimal constants in these results. Throughout this section, we assume $f\in\mathcal{S}_{\mu,L}$ for some $0<\mu\leq L$. Recall that function $f\in\mathcal{S}_{\mu,L}$ if $f\in\mathcal{F}_{L}$ and $f(x)-\mu\|x\|^{2}/2$ is convex.

We first define the speed restarting time. For the ODE (3), we call

the speed restarting time. In words, $T$ is the first time the velocity $\|\dot{X}\|$ decreases. Back to the discrete scheme, it is the first time when we observe $\|x_{k+1}-x_{k}\|<\|x_{k}-x_{k-1}\|$. This definition itself does not directly imply that $0<T<\infty$, which is proven later in Lemmas 13 and 25. Indeed, $f(X(t))$ is a decreasing function before time $T$; for $t\leq T$,

There exist positive constants $c_{1}$ and $c_{2}$, which only depend on the condition number $L/\mu$, such that for any $f\in\mathcal{S}_{\mu,L}$, we have

2 Proof of Linear Convergence

First, we collect some useful estimates. Denote by $M(t)$ the supremum of $\|\dot{X}(u)\|/{u}$ over $u\in(0,t]$ and let

It is guaranteed that $M$ defined above is finite, for example, see the proof of Lemma 18. The definition of $M$ gives a bound on the gradient of $f$,

Hence, it is easy to see that $I$ can also be bounded via $M$,

To fully facilitate these estimates, we need the following lemma that gives an upper bound of $M$, whose proof is deferred to the appendix.

Next we give a lemma which claims that the objective function decays by a constant through each speed restarting.

There is a universal constant $C>0$ such that

where $C_{1}>0$ is an absolute constant and the second inequality follows from Lemma 25 in the appendix. Consequently,

where $C=16C_{1}/25$ and in the last inequality we use the $\mu$-strong convexity of $f$. Thus we have

To complete the proof, note that $f(X(T))\leq f(X(4/(5\sqrt{L})))$ by Lemma 25.

With each restarting reducing the error $f-f^{\star}$ by a constant a factor, we still need the following lemma to ensure sufficiently many restartings.

Plugging $t=4/(5\sqrt{L})$ into (32) gives $\|\dot{X}(4/(5\sqrt{L}))\|\geq\frac{C_{1}}{\sqrt{L}}\|\nabla f(x_{0})\|$ for some universal constant $C_{1}>0$. Hence, from the last display we get

Now, we are ready to prove Theorem 10 by applying Lemmas 12 and 13.

In closing, we remark that we believe that estimate in Lemma 12 is tight, while not for Lemma 13. Thus we conjecture that for a large class of $f\in\mathcal{S}_{\mu,L}$, if not all, $T=O(\sqrt{L}/\mu)$. If this is true, the exponent constant $c_{2}$ in Theorem 10 can be significantly improved.

3 Numerical Examples

Below we present a discrete analog to the restarted scheme. There, $k_{\min}$ is introduced to avoid having consecutive restarts that are too close. To compare the performance of the restarted scheme with the original (1), we conduct four simulation studies, including both smooth and non-smooth objective functions. Note that the computational costs of the restarted and non-restarted schemes are the same.

Quadratic. $f(x)=\frac{1}{2}x^{T}Ax+b^{T}x$ is a strongly convex function, in which $A$ is a $500\times 500$ random positive definite matrix and $b$ a random vector. The eigenvalues of $A$ are between $0.001$ and $1$. The vector $b$ is generated as i.i.d. Gaussian random variables with mean 0 and variance 25.

where $n=50,m=200,\rho=20$. The matrix $A=(a_{ij})$ is a random matrix with i.i.d. standard Gaussian entries, and $b=(b_{i})$ has i.i.d. Gaussian entries with mean $0$ and variance $2$. This function is not strongly convex.

Lasso. $f(x)=\frac{1}{2}\|Ax-b\|^{2}+\lambda\|x\|_{1}$, where $A$ is a $1000\times 500$ random matrix and $b$ is given as $b=Ax^{0}+z$ for 20-sparse $x^{0}$ and Gaussian noise $z$. We set $\lambda=1.5\sqrt{2\log p}$.

In these examples, $k_{\min}$ is set to be 10 and the step sizes are fixed to be $1/L$. If the objective is in composite form, the Lipschitz bound applies to the smooth part. Figure 6 presents the performance of the speed restarting scheme, the gradient restarting scheme, the original Nesterov’s scheme and the proximal gradient method. The objective functions include strongly convex, non-strongly convex and non-smooth functions, violating the assumptions in Theorem 10. Among all the examples, it is interesting to note that both speed restarting scheme empirically exhibit linear convergence by significantly reducing bumps in the objective values. This leaves us an open problem of whether there exists provable linear convergence rate for the gradient restarting scheme as in Theorem 10. It is also worth pointing out that compared with gradient restarting, the speed restarting scheme empirically exhibits more stable linear convergence rate.

Discussion

This paper introduces a second-order ODE and accompanying tools for characterizing Nesterov’s accelerated gradient method. This ODE is applied to study variants of Nesterov’s scheme and is capable of interpreting some empirically observed phenomena, such as oscillations along the trajectories. Our approach suggests (1) a large family of generalized Nesterov’s schemes that are all guaranteed to converge at the rate $O(1/k^{2})$, and (2) a restarting scheme provably achieving a linear convergence rate whenever $f$ is strongly convex.

In this paper, we often utilize ideas from continuous-time ODEs, and then apply these ideas to discrete schemes. The translation, however, involves parameter tuning and tedious calculations. This is the reason why a general theory mapping properties of ODEs into corresponding properties for discrete updates would be a welcome advance. Indeed, this would allow researchers to only study the simpler and more user-friendly ODEs.

As evidenced by many examples, the viewpoint of regarding the ODE as a surrogate for Nesterov’s scheme would allow a new perspective for studying accelerated methods in optimization. The discrete scheme and the ODE are closely connected by the exact mapping between the coefficients of momentum (e.g. $(k-1)/(k+2)$) and velocity (e.g. $3/t$). The derivations of generalized Nesterov’s schemes and the speed restarting scheme are both motivated by trying a different velocity coefficient, in which the surprising phase transition at 3 is observed. Clearly, such alternatives are endless, and we expect this will lead to findings of many discrete accelerated schemes. In a different direction, a better understanding of the trajectory of the ODEs, such as curvature, has the potential to be helpful in deriving appropriate stopping criteria for termination, and choosing step size by backtracking.

Appendix A. Proof of Theorem 1

The proof is divided into two parts, namely, existence and uniqueness.

Below, some preparatory lemmas are given before turning to the proof of this lemma. To begin with, for any $\delta>0$ consider the smoothed ODE

with $X(0)=x_{0},\dot{X}(0)=0$. Denoting by $Z=\dot{X}$, then (34) is equivalent to

with $X(0)=x_{0},Z(0)=0$. As functions of $(X,Z)$, both $Z$ and $-3Z/\max(\delta,t)-\nabla f(X))$ are $\max(1,L)+3/\delta$-Lipschitz continuous. Hence by standard ODE theory, (34) has a unique global solution in $C^{2}[0,\infty)$, denoted by $X_{\delta}$. Note that $\ddot{X}_{\delta}$ is also well defined at $t=0$. Next, introduce $M_{\delta}(t)$ to be the supremum of $\|\dot{X}_{\delta}(u)\|/u$ over $u\in(0,t]$. It is easy to see that $M_{\delta}(t)$ is finite because $\|\dot{X}_{\delta}(u)\|/u=(\|\dot{X}_{\delta}(u)-\dot{X}_{\delta}(0)\|)/u=\|\ddot{X}_{\delta}(0)\|+o(1)$ for small $u$. We give an upper bound for $M_{\delta}(t)$ in the following lemma.

The proof of Lemma 15 relies on a simple lemma.

For any $u>0$, the following inequality holds

For $0<t\leq\delta$, the smoothed ODE takes the form

Taking the supremum of $\|\dot{X}_{\delta}(t)\|/t$ over $0<t\leq\delta$ and rearranging the inequality give the desired result. ∎

Next, we give an upper bound for $M_{\delta}(t)$ when $t>\delta$.

For $\delta<\sqrt{6/L}$ and $\delta<t<\sqrt{12/L}$, we have

For $t>\delta$, the smoothed ODE takes the form

Hence, by integration, $t^{3}\dot{X}_{\delta}(t)$ is equal to

where the last expression is an increasing function of $t$. So for any $\delta<t^{\prime}<t$, it follows that

which also holds for $t^{\prime}\leq\delta$. Taking the supremum over $t^{\prime}\in(0,t)$ gives

The desired result follows from rearranging the inequality. ∎

By Lemmas 15 and 17, for any $t\in[0,\sqrt{6/L}],\delta\in(0,\sqrt{3/L})$ the gradient is uniformly bounded as

Thus it immediately implies that $\mathcal{F}$ is equicontinuous. To establish the uniform boundedness, note that

We are now ready for the proof of Lemma 14.

as $i\rightarrow\infty$. However, by definition we have

The first condition $\breve{X}(0)=x_{0}$ is a direct consequence of $X_{\delta_{m_{i}}}(0)=x_{0}$. To check the second, pick a small $t>0$ and note that

where $\xi_{i}\in(0,t)$ is given by the mean value theorem. The desired result follows from taking $t\rightarrow 0$. ∎

Next, we aim to prove the uniqueness of the solution to (3).

For any $f\in\mathcal{F}_{\infty}$, the ODE (3) has at most one local solution in a neighborhood of $t=0$.

By Lipschitz continuity of the gradient, one has

Given all of the aforementioned lemmas, the proof of Theorem 1 is simply combining 14 and 19.

Appendix B. Proof of Theorem 2

Identifying $\sqrt{s}=\Delta t$, the comparison between (4) and (15) reveals that Nesterov’s scheme is a discrete scheme for numerically integrating the ODE (3). However, its singularity of the damping coefficient at $t=0$ leads to the nonexistence of off-the-shelf ODE theory for proving Theorem 2. To address this difficulty, we use the smoothed ODE (34) to approximate the original one; then bound the difference between Nesterov’s scheme and the forward Euler scheme of (34), which may take the following form:

with $X_{0}^{\delta}=x_{0}$ and $Z_{0}^{\delta}=0$.

With step size $\Delta t=\sqrt{s}$, for any $T>0$ we have

Let $z_{k}=(x_{k+1}-x_{k})/\sqrt{s}$. Then Nesterov’s scheme is equivalent to

Denote by $a_{k}=\|X_{k}^{\delta}-x_{k}\|,\quad b_{k}=\|Z_{k}^{\delta}-z_{k}\|$, whose initial values are $a_{0}=0$ and $b_{0}=\|\nabla f(x_{0})\|\sqrt{s}$. The idea of this proof is to bound $a_{k}$ via simultaneously estimating $a_{k}$ and $b_{k}$. By comparing (37) and (38), we get the iterative relationship for $a_{k}$: $a_{k+1}\leq a_{k}+\sqrt{s}b_{k}$. Denoting by $S_{k}=b_{0}+b_{1}+\cdots+b_{k}$, this yields

Similarly, for sufficiently small $s$ we get

To upper bound $\|z_{k}\|$, denoting by $C_{1}$ the supremum of $\sqrt{2L(f(y_{k})-f^{\star})}$ over all $k$ and $s$, we have

which gives $\|z_{k}\|\leq C_{1}(k+1)\sqrt{s}$. Hence,

By induction on $k$, for $k\leq\delta/\sqrt{s}$ it holds that

Letting $k^{\star}=\lfloor\delta/\sqrt{s}\rfloor$, we get

Next, we bound $b_{k}$ for $k>k^{\star}=\lfloor\delta/\sqrt{s}\rfloor$. To this end, we consider the worst case of (40), that is,

for $k>k^{\star}$ and $S_{k^{\star}}=S_{k^{\star}+1}=C_{3}\delta^{2}/\sqrt{s}+o_{s}(1/\sqrt{s})$ for some sufficiently large $C_{3}$. In this case, $C_{2}s/\delta<sS_{k-1}$ for sufficiently small $s$. Hence, the last display gives

Letting $k^{\diamond}=\lfloor T/\sqrt{s}\rfloor$, we further get

for $k^{\star}<k\leq k^{\diamond}$. Last, combining (41) and the last display, we get the desired result.

where $X_{\delta}(\cdot)$ is the solution to the smoothed ODE (34). The proof of Lemma 14 implies that, we can choose a sequence $\delta_{m}\rightarrow 0$ such that

The second term $\|X^{\delta_{m}}_{k}-X_{\delta_{m}}(k\sqrt{s})\|$ will uniformly vanish as $s\rightarrow 0$ and so does the first term $\|x_{k}-X^{\delta_{m}}_{k}\|$ if first $s\rightarrow 0$ and then $\delta_{m}\rightarrow 0$. This completes the proof. ∎

Appendix C. ODE for Composite Optimization

In analogy to (3) for smooth $f$ in Section 2, we develop an ODE for composite optimization,

where $g\in\mathcal{F}_{L}$ and $h$ is a general convex function possibly taking on the value $+\infty$. Provided it is easy to evaluate the proximal of $h$, Beck and Teboulle 2009 propose a proximal gradient version of Nesterov’s scheme for solving (42). It is to repeat the following recursion for $k\geq 1$,

where the proximal subgradient $G_{s}$ has been defined in Section 4.1. If the constant step size $s\leq 1/L$, it is guaranteed that (Beck and Teboulle 2009)

which in fact is a special case of Theorem 6.

Compared to the smooth case, it is not as clear to define the driving force as $\nabla f$ in (3). At first, it might be a good try to define

if it exists. However, as implied in the proof of Theorem 24 stated below, this definition fails to capture the directional aspect of the subgradient. To this end, we define the subgradients through the following lemma.

Note that the directional derivative is semilinear in $p$ because

and $G(x,0)$ to be any element in $\partial f(x)$. Now we give the main theorem. However, note that we do not guarantee the existence of solution to (43).

Given a convex function $f(x)$ with directional subgradient $G(x,p;f)$, assume that the second order ODE

admits a solution $X(t)$ on $[0,\alpha)$ for some $\alpha>0$. Then for any $0<t<\alpha$, we have

It suffices to establish that $\mathcal{E}$, first defined in the proof of Theorem 3, is monotonically decreasing. The difficulty comes from that $\mathcal{E}$ may not be differentiable in this setting. Instead, we study $(\mathcal{E}(t+\Delta t)-\mathcal{E}(t))/\Delta t$ for small $\Delta t>0$. In $\mathcal{E}$, the second term $2\|X+t\dot{X}/2-x^{\star}\|^{2}$ is differentiable, with derivative $4\langle X+\frac{t}{2}\dot{X}-x^{\star},\frac{3}{2}\dot{X}+\frac{t}{2}\ddot{X}\rangle$. Hence,

Since $f$ is locally Lipschitz, $o(\Delta t)$ term does not affect the function in the limit,

Combining all of (44), (46) and (47), we obtain

where the last inequality follows from the convexity of $f$. Thus,

which along with the continuity of $\mathcal{E}$, concludes that $\mathcal{E}(t)$ is a non-increasing function of $t$.

We give a simple example as follows. Consider the Lasso problem

To encourage sparsity, for any index $i$ with $x_{i}=0,p_{i}=0$, we let