On the Curved Geometry of Accelerated Optimization

TL;DR

This paper introduces a differential geometric perspective on Nesterov's accelerated gradient method, modeling it as a proximal point method on a Riemannian manifold, and analyzes its convergence for quadratic problems.

Contribution

It offers a novel geometric interpretation of AGM using Riemannian manifolds and extends the analysis to continuous-time ODEs on curved spaces.

Findings

AGM can be viewed as a proximal point method on a Riemannian manifold.

The continuous-time limit of AGM corresponds to an ODE on the manifold.

Convergence rate analysis is provided for quadratic objectives.

Abstract

In this work we propose a differential geometric motivation for Nesterov's accelerated gradient method (AGM) for strongly-convex problems. By considering the optimization procedure as occurring on a Riemannian manifold with a natural structure, The AGM method can be seen as the proximal point method applied in this curved space. This viewpoint can also be extended to the continuous time case, where the accelerated gradient method arises from the natural block-implicit Euler discretization of an ODE on the manifold. We provide an analysis of the convergence rate of this ODE for quadratic objectives.

Figures

Tables

Table 1. Table 1: Equivalent forms of Nesterov’s method for μ 𝜇 \mu -strongly convex, L 𝐿 L -smooth f 𝑓 f . Proofs of the stated relations are available in the appendix.

Form Name	Algorithm	Relations
Nesterov (2013)
form I
Nesterov (2013)
form II		${\beta }_{\text{Nes}}=\frac{\sqrt{L}-\sqrt{\mu }}{\sqrt{L}+\sqrt{\mu }}$
Sutskever et al. (2013)
Modern Momentum111PyTorch & Tensorflow (for instance) implement this form. Evaluating the gradient and function at the current iterate $x^k$ instead of a shifted point makes it more consistent with gradient descent when wrapped in a generic optimization layer.
Auslender and Teboulle (2006)
Lan and Zhou (2017)

Equations

$$x^{k+1}=x^k-\gamma \nabla f\left(x^k\right)+\beta \left(x^k-x^{k-1}\right),$$

$$x^k=\arg \underset{x}{\min }\left\{f(x)+\frac{\eta }{2}{\Vert x-x^{k-1}\Vert }^2\right\}.$$

$$\frac{d^2{\gamma }^i}{d{t}^2}\doteq \frac{d^2{x}^i}{d{t}^2}+\sum _{j,k}{\Gamma }_{jk}^i(x)\frac{d{x}^j}{dt}\frac{d{x}^k}{dt}=0,$$

$$B_{\phi}(x,y)=\phi(x)-\phi(y)-\bigl{\langle}\nabla\phi(y),x-y\bigr{\rangle},$$

$$g_x(u,v)=u^{T}H(x)v.$$

$$y=\nabla \varphi (x).$$

$$u=H(x)v,$$

$$g_x(v_1,v_2)=v_1^{T}H(x)v_2=v_1^{T}H(x)\left(H(x{)}^{-1}{u}_2\right)=v_1^T{u}_2,$$

$$\nabla \varphi \left(\gamma (t)\right)=at+b,$$

$$\gamma (t)=\nabla {\varphi }^{\ast }\left(t\nabla \varphi (y)+\left(1-t\right)\nabla \varphi (x)\right).$$

$$y=\nabla {\varphi }^{\ast }\left(\nabla \varphi (x)+H(x^k)v\right).$$

$${\Gamma }_{ij}^{(D)k}(x)=2{\Gamma }_{ij}^{(R)k}={\left[H(x{)}^{-1}{\left(\nabla H(x)\right)}_i\right]}_{kj},$$

$$x^k=\arg \underset{x}{\min }\left\{f(x)+\rho B_{\varphi }(x,x^{k-1})\right\}.$$

$$x^k=\arg \underset{x}{\min }\left\{⟨x,\nabla f(x^{k-1})⟩+\rho B_{\varphi }(x,x^{k-1})\right\}.$$

$$-\frac{1}{\rho }\nabla f(x^k)=\nabla \varphi (x^k)-\nabla \varphi (x^{k-1}).$$

$${\textstyle x^{k}=\arg\min_{x}\bigl{\{}f(x)+\frac{\rho}{2}\left\|x-x^{k-1}\right\|^{2}\bigr{\}}.}$$

$$x^k=\arg \underset{x}{\min }\left\{f(x)+\frac{\eta }{2}{\Vert x-x^{k-1}\Vert }^2\right\}.$$

$$x^k=x^{k-1}-\frac{1}{\eta }\nabla f(x^k),$$

$$x^k+\frac{1}{\eta }\nabla f(x^k)=x^{k-1}.$$

$$\nabla f^{\ast }\left(g^k\right)+\frac{1}{\eta }{g}^k=x^{k-1},$$

$$g^k=\arg \underset{g}{\min }\left\{f^{\ast }(g)+\frac{1}{2\eta }{\Vert g-\eta x^{k-1}\Vert }^2\right\}.$$

$$\nabla f^{\ast }\left(g^k\right)+\frac{1}{\eta }{g}^k=\frac{1}{\eta }{g}^{k-1}+\left(x^{k-1}-\frac{1}{\eta }{g}^{k-1}\right),$$

$$g^k$$

$$x^k$$

$$g^k=\arg \underset{g}{\min }\left\{f^{\ast }(g)-⟨g,x^{k-1}-\frac{1}{\eta }{g}^{k-1}⟩+\tau B_{f^{\ast }}(g,g^{k-1})\right\}.$$

$$\nabla f^{\ast }(g^k)-\left[x^{k-1}-\frac{1}{\eta }{g}^{k-1}\right]+\tau \nabla f^{\ast }(g^k)-\tau \nabla f^{\ast }(g^{k-1})=0,$$

$$\text{therefore }{g}^k=\nabla f\left({\left(1+\tau \right)}^{-1}\left[x^{k-1}-\frac{1}{\eta }{g}^{k-1}+\tau \nabla f^{\ast }(g^{k-1})\right]\right).$$

$$y^k=\frac{x^{k-1}-\frac{1}{\eta }{g}^{k-1}+\tau y^{k-1}}{1+\tau }.$$

$$y^k$$

$$g^k$$

$$x^k$$

$$\eta =\sqrt{\mu L},\tau =\frac{L}{\eta },\text{α=1+ττ}.$$

Full text

1 Introduction

The core algorithms of convex optimization are gradient descent (GD) and the accelerated gradient method (AGM). These methods are rarely used directly, more often they occur as the building blocks for distributed, composite, or non-convex optimization. In order to build upon these components, it is helpful to understand not just how they work, but why. The gradient method is well understood in this sense. It is commonly viewed as following a direction of steepest descent or as minimizing a quadratic upper bound. These interpretations provide a motivation for the method as well as suggesting a potential convergence proof strategy.

The accelerated gradient method in contrast has an identity crisis. Its equational form is remarkably malleable, allowing for many different ways of writing the same updates. We list a number of these forms in Table 1. Nesterov’s original motivation for the AGM method used the concept of estimate sequences. Unfortunately, estimate sequences do not necessarily yield the simplest accelerated methods when generalized, such as for the composite case (Beck and Teboulle 2009; Nesterov 2007), and they have not been successfully applied in the important finite-sum (variance reduced) optimization setting.

Because of the complexity of estimate sequences, the AGM method is commonly motivated as a form of momentum. This view is flawed as a way of introducing the AGM method from first principles, as the most natural way of using momentum yields the heavy ball method instead:

[TABLE]

which arises from discretizing the physics of a particle in a potential well with additional friction. The heavy-ball method does not achieve an accelerated convergence rate on general convex problems, suggesting that momentum, per se, is not the reason for acceleration. Another contemporary view is the linear-coupling interpretation of Allen-Zhu and Orecchia (2017), which views the AGM method as an interpolation between gradient descent and mirror descent. We take a more geometric view in our interpretation.

In this work we motivate the AGM by introducing it as an application of the proximal-point method:

[TABLE]

The proximal point (PP) method is perhaps as foundational as the gradient descent method, although it sees even less direct use as each step requires solving a regularized subproblem, in contrast to the closed form steps for GD and AGM. The PP method gains remarkable convergence rate properties in exchange for the computational difficulty, including convergence for any positive step-size.

We construct the AGM by applying a dual form of the proximal point method in a curved space. Each step follows a geodesic on a manifold in a sense we make precise in Section 4. We use the term curved with respect to a coordinate system, rather than a coordinate free notion of curvature such as the Riemannian curvature. We first give a brief introduction to the concepts from differential geometry necessary to understand our motivation. The equational form that our argument yields is much closer to those that have been successfully applied in practice, particularly for the minimization of finite sums (Lan and Zhou, 2017; Zhang and Xiao, 2017).

2 Connections

An (affine) connection is a type of structure on a manifold that can be used to define and compute geodesics. Geodesics in this sense represent curves of zero acceleration. These geodesics are more general concepts than Riemannian geodesics induced by the Riemannian connection, not necessarily representing the shortest path in any metric. Indeed, we will define multiple connections on the same manifold that lead to completely different geodesics.

Given a $n$ dimensional coordinate system, a connection is defined by $n^{3}$ numbers at every point $x$ on the manifold, called the connection coefficients (or Christoffel symbols) $\varGamma_{ij}^{k}(x)$. A geodesic is a path $\gamma:[0,1]\rightarrow\mathcal{M}$ (in our case $\mathcal{M}=\mathbb{R}^{n}$) between two points $x$ and $y$ can then be computed as the unique solution $\gamma(t)=x(t)$ to the system of ordinary differential equations (Lee, 1997, Page 58, Eq 4.11):

[TABLE]

with boundary conditions $x(0)=x$ and $x(1)=y.$ Here $x^{i}$ denotes the $i\text{th}$ component of $x$ expressed in the same coordinate system as the connection.

3 Divergences induce Hessian manifold structure

Let $\phi$ be a smooth strongly convex function defined on $\mathbb{R}^{n}$. The Bregman divergence generated by $\phi$:

[TABLE]

and its derivatives can be used to define a remarkable amount of structure on the domain of $\phi$. In particular, we can define a Riemannian manifold, together with two dually flat connections with biorthogonal coordinate systems (Amari and Nagaoka, 2000; Shima, 2007). This structure is also known as a Hessian manifold. Topologically it is $\mathcal{M}=\mathbb{R}^{n}$ with the following additional geometric structures.

Riemannian structure

Riemannian manifolds have the additional structure of a metric tensor (a generalized dot-product), defined on their tangent spaces. We denote the vector space of tangent vectors at a point $x$ as $T_{x}\mathcal{M}$. If we express the tangent vectors with respect to the Euclidean basis, the metric at a point $x$ is a quadratic form with the Hessian matrix $H(x)=\nabla_{x}^{2}B(x,y)=\nabla^{2}\phi(x)$ of $\phi$ at $x$:

[TABLE]

Biorthogonal coordinate systems

Central to the notion of a manifold is the invariance to the choice of coordinate system. We can express a point on the manifold as well as a point in the tangent space using any coordinate system that is most convenient. Of course, when we wish to perform calculations on the manifold we must be careful to express all quantities in that coordinate system. Euclidean coordinates $e_{i}$ are the most natural on our Hessian manifold, however there is another coordinate system which is naturally dual to $e_{i},$ and ties the manifold structure directly to duality theory in optimization.

Recall that for a convex function $\phi$ we may define the convex conjugate $\phi^{*}(y)=\max_{x}\left\{\left\langle x,y\right\rangle-\phi(x)\right\}.$ The dual coordinate system we define simply identifies each point $x$, when expressed in Euclidean (“primal”) coordinates, with the vector of “dual” coordinates:

[TABLE]

Our assumptions of smoothness and strong convexity imply this is a one-to-one mapping, with inverse given by $x=\nabla\phi^{*}(y).$ The remarkable fact that the gradient of the conjugate is the inverse of the gradient is a key building block of the theory in this paper.

The notion of biorthogonality refers to natural tangent space coordinates of these two systems. A tangent vector $v$ at a point $x$ can be converted to a vector $u$ of dual (tangent space) coordinates using matrix multiplication with the Hessian (Shima, 2007):

[TABLE]

Given the definition of the metric above, it is easy to see that if we have two vectors $v_{1}$ and $v_{2}$, we may express $v_{2}$ in dual coordinates $u_{2}$ so that the metric tensor takes the simple form:

[TABLE]

which is the biorthogonal relation between the two tangent space coordinate systems.

Dual Flat Connections

There is an obvious connection $\varGamma^{(E)}$ we can apply to the Hessian manifold, the Euclidean connection that trivially identifies straight lines in $\mathbb{R}^{n}$ as geodesics. Normally when we perform gradient descent in $\mathbb{R}^{n}$ we are implicitly following a geodesic of this connection. The connection coefficients $\varGamma_{ij}^{(E)k}$ are all zero when this connection is expressed in Euclidean coordinates. A connection that has $\varGamma_{ij}^{k}=0$ with respect to some coordinate system is a flat connection.

The Hessian manifold admits another flat connection, which we will call the dual connection, as it corresponds to straight lines in the dual coordinate system established above. In particular each dual geodesic can be expressed in primal coordinates $\gamma(t)$ as a solution to the equation:

[TABLE]

for vectors $a$, $b$ representing the initial velocity and point respectively (both represented in dual coordinates) that depend on the boundary conditions. This is quite easy to solve using the relation $\nabla\phi^{-1}=\nabla\phi^{*}$ discussed above. For instance, a geodesic $\gamma:[0,1]\rightarrow\mathcal{M}$ between two arbitrary points $x$ and $y$ under the dual connection could be computed explicitly in Euclidean coordinates as:

[TABLE]

If we instead know the initial velocity we can find the endpoint with:

[TABLE]

The flatness of the dual connection $\varGamma^{(D)}$ is crucial to its computability in practice. If we instead try to compute the geodesic in Euclidean coordinates using the geodesic ODE, we have to work with the connection coefficients which involve third derivatives of $\phi$ (taking the form of double those of the Riemannian connection $\varGamma^{(R)}$):

[TABLE]

The Riemannian connection’s geodesics are similarly difficult to compute directly from the ODE (they also can’t generally be expressed in a simpler form).

4 Bregman proximal operators follow geodesics

Bregman divergences arise in optimization primarily through their use in proximal steps. A Bregman proximal operation balances finding a minimizer of a given function $f$ with maintaining proximity to a given point $y$, measured using a Bregman divergence instead of a distance metric:

[TABLE]

A core application of this would be the mirror descent step (Nemirovski and Yudin, 1983; Beck and Teboulle, 2003), where the operation is applied to a linearized version of $f$ instead of $f$ directly:

[TABLE]

Bregman proximal operations can be interpreted as geodesic steps with respect to the dual connection. The key idea is that given an input point $x^{k-1}$, they output a point $x$ such that the velocity of the connecting geodesic is equal to $-\nabla\frac{1}{\rho}f(x)$ at $x$. This velocity is measured in the flat coordinate system of the connection, the dual coordinates. To see why, consider a geodesic $\gamma(t)=(1-t)\nabla\phi(x^{k-1})+t\nabla\phi(x^{k}).$ Here $x^{k-1}$ and $x^{k}$ are in primal coordinates and $\gamma(t)$ is in dual coordinates. The velocity is $\frac{d}{dt}\gamma(t)=\nabla\phi(x^{k})-\nabla\phi(x^{k-1}).$ Contrast to the optimality condition of the Bregman prox (Equation 3):

[TABLE]

For instance, when using the Euclidean penalty the step is:

[TABLE]

The final velocity is just $x^{k}-x^{k-1}$, and so $x^{k}-x^{k-1}=-\frac{1}{\rho}\nabla f(x^{k})$, which is the solution of the proximal operation.

5 Primal-Dual form of the proximal point method

The proximal point method is the building block from which we will construct the accelerated gradient method. Consider the basic form of the proximal point method applied to a strongly convex function $f$. At each step, the iterate $x^{k}$ is constructed from $x^{k-1}$ by solving the proximal operation subproblem given an inverse step size parameter $\eta$:

[TABLE]

This step can be considered an implicit form of the gradient step, where the gradient is evaluated at the end-point of the step instead of the beginning:

[TABLE]

which is just the optimality condition of the subproblem in Equation 4, found by taking the derivative $\nabla f(x)+\eta x-\eta x^{k-1}$ to be zero. A remarkable property of the proximal operation becomes apparent when we rearrange this formula, namely that the solution to the operation is not a single point but a primal-dual pair, whose weighted sum is equal to the input point:

[TABLE]

If we define $g^{k}=\nabla f(x^{k})$, the primal-dual pair obeys a duality relation: $g^{k}=\nabla f(x^{k})$ and $x^{k}=\nabla f^{*}(g^{k})$, which allows us to interchange primal and dual quantities freely. Indeed we may write the condition in a dual form as:

[TABLE]

which is the optimality condition for the proximal operation:

[TABLE]

Our goal in this section is to express the proximal point method in terms of a dual step, and while this equation involves the dual function $f^{*}$, it is not a step in the sense that $g^{k}$ is formed by a proximal operation from $g^{k-1}.$

We can manipulate this formula further to get an update of the form we want, by simply adding and subtracting $g^{k-1}$ from 5:

[TABLE]

Which gives the updates:

[TABLE]

We call this the primal-dual form of the proximal point method.

6 Acceleration as a change of geometry

The proximal point method is rarely used in practice due to the difficulty of computing the solution to the proximal subproblem. It is natural then to consider modifications of the subproblem to make it more tractable. The subproblem becomes particularly simple if we replace the proximal operation with a Bregman proximal operation with respect to $f^{*}$,

[TABLE]

We have additionally changed the penalty parameter to a new constant $\tau$, which is necessary as the change to the Bregman divergence changes the scaling of distances. We discuss this further below.

Recall from Section 4 that Bregman proximal operations follow geodesics. The key idea is that we are now following a geodesic in the dual connection of $\phi=f^{*}$, using the notation of Section 1, which is a straight-line in the primal coordinates of $f$ due to the flatness of the connection (Section 1). Due to the flatness property, a simple closed-form solution can be derived by equating the derivative to 0:

[TABLE]

This formula gives $g^{k}$ in terms of the derivative of known quantities, as $\nabla f^{*}(g^{k-1})$ is known from the previous step as the point at which we evaluated the derivative at. We will denote this argument to the derivative operation $y$, so that $g^{k}=\nabla f(y^{k})$. It no longer holds that $g^{k}=\nabla f(x^{k})$ after the change of divergence. Using this relation, $y$ can be computed each step via the update:

[TABLE]

In order to match the accelerated gradient method exactly we need some additional flexibility in the step size used in the $y^{k}$ update. To this end we introduce an additional constant $\alpha$ in front of $g^{k-1}$, which is 1 for the proximal point variant. The full method is as follows:

Bregman form of the accelerated gradient method

$\displaystyle y^{k}$ $\displaystyle=\frac{x^{k-1}-\frac{\alpha}{\eta}g^{k-1}+\tau y^{k-1}}{1+\tau},$

$\displaystyle g^{k}$ $\displaystyle=\nabla f(y^{k}),$

$\displaystyle x^{k}$ $\displaystyle=x^{k-1}-\frac{1}{\eta}g^{k}.$

(6)

This is very close to the equational form of Nesterov’s method explored by Lan and Zhou (2017), with the change that they assume an explicit regularizer is used, whereas we assume strong convexity of $f$. Indeed we have chosen our notation so that the constants match. This form is algebraically equivalent to other known forms of the accelerated gradient method for appropriate choice of constants. Table 1 shows the direct relation between the many known ways of writing the accelerated gradient method in the strongly-convex case (Proofs of these relations are in the Appendix). When $f$ is $\mu$-strongly convex and $L$-smooth, existing theory implies an accelerated geometric convergence rate of at least $1-\sqrt{\frac{\mu}{L}}$ for the parameter settings (Nesterov, 2013):

[TABLE]

In contrast, the primal-dual form of the proximal point method achieves at least that convergence rate for parameters:

[TABLE]

The difference in $\tau$ arises from the difference in the scaling of the Bregman penalty compared to the Euclidean penalty. The Bregman generator $f^{*}$ is strongly convex with constant $1/L$ whereas the Euclidean generator $\frac{1}{2}\left\|\cdot\right\|^{2}$ is strongly convex with constant 1, so the change in scale requires rescaling by $L$.

6.1 Interpretations

After the change in geometry, the $g$ update no longer gives a dual point that is directly the gradient of the primal iterate. However, notice that the term we are attempting to minimize in the $g$ step:

[TABLE]

has a fixed point of $\nabla f^{*}\left(g^{k}\right)=x^{k-1}-\frac{\alpha}{\eta}g^{k},$ which is precisely an $\alpha$-weight version of the proximal point’s key property from Equation 5. Essentially we have relaxed the proximal-point method. Instead of this relation holding precisely at every step, we are instead constantly taking steps which pull $g$ closer to satisfying it.

6.2 Inertial form

The primal-dual view of the proximal point method can also be written in terms of the quantity $z^{k-1}=x^{k-1}-\frac{\alpha}{\eta}g^{k-1}$ instead of $x^{k-1}$. This form is useful for the construction of ODEs that model the discrete dynamics. Under this change of variables the updates are:

[TABLE]

6.3 Relation to the heavy ball method

Consider Equation 6 with $\alpha=0$, which removes the over-extrapolation before the proximal operation. If we define $\beta=\frac{\tau}{1+\tau}$ we may write the method as:

[TABLE]

We can eliminate $x^{k}$ from the $y^{k}$ update above by plugging in the $x^{k}$ step equation, then using the $y^{k}$ update from the previous step in the form $\left(1-\beta\right)x^{k-1}=y^{k-1}-\beta y^{k-2}:$

[TABLE]

This has the exact form of the heavy ball method with step size $\left(1-\beta\right)/\eta$ and momentum $\beta$. We can also derive the heavy ball method by starting from the saddle-point expression for $f$:

[TABLE]

The alternating-block gradient descent/ascent method on the objective $\left\langle x,g\right\rangle-f^{*}(g)$ with step-size $\gamma$ is simply:

[TABLE]

If we instead perform a Bregman proximal update in the dual geometry for the $g$ part, we arrive at the same equations as we had for the primal-dual proximal point method but with $\alpha=0$, yielding the heavy ball method. In order to get the accelerated gradient method instead of the heavy ball method, the extra inertia that arises from starting from the proximal point method instead of the saddle point formulation appears to be crucial.

7 Dual geometry in continuous time

The inertial form (Equation 7) of the proximal point method can be formulated as an ODE in a very natural way, by mapping $z^{k}-z^{k-1}\rightarrow\dot{z}$ and $g^{k}-g^{k-1}\rightarrow\dot{g}$, and taking $x$ and $g$ to be at time $t$. This is the inverse of the Euler class of discretizations applied separately to the two terms, which is the most natural way to discretize an ODE. The resulting proximal point ODE is:

[TABLE]

We have suppressed the dependence on $t$ of each quantity for notational simplicity. We can treat $g$ more formally as a point $g\in\mathcal{M}$ on a Hessian manifold $\mathcal{M}$. Then the solution for the $g$ variable of the ODE is a curve $\gamma(t):I\rightarrow\mathcal{TM}$ from an interval $I$ to the tangent bundle on the manifold so the velocity $\dot{\gamma}(t)\in T_{g}\mathcal{M}$ obeys the ODE: $\dot{\gamma}(t)=f_{g}(z,g,t).$ The right hand side of the ODE is a point in the tangent space of the manifold at $\gamma(t)$, expressed in Euclidean coordinates.

We can now apply the same partial change of geometry that we used in the discrete case. We will consider the quantity $f_{g}(z,g,t)$ to be a tangent vector in dual tangent space coordinates For the $\phi=f^{*}$ Hessian manifold, instead of its primal tangent space coordinates (which would leave the ODE unchanged). The variable $g$ remains in primal coordinates with respect to $\phi$, so we must add to the ODE a change of coordinates for the tangent vector, yielding:

[TABLE]

where we have used the inverse of Equation 1, with $\phi=f^{*}$. We can rewrite this as:

[TABLE]

giving the AGM ODE system:

[TABLE]

It is now easily seen that the implicit Euler update for the $g$ variable with $z$ fixed now corresponds to the solution of the Bregman proximal operation considered in the discrete case. So this ODE is a natural continuous time analogue to the accelerated gradient method.

Convergence in continuous time

The natural analogy to convergence in continuous time is known as the decay rate of the ODE. A sufficient condition for an ODE with parameters $u=[z;g]$ to decay with constant $\rho$ is:

[TABLE]

where $u^{*}$ is a fixed point. We can relate this to the discrete case by noting that $\exp(-t\rho)=\lim_{k\rightarrow\infty}(1-\frac{t}{k}\rho)^{k}$, so given our discrete-time convergence rate is proportional to $(1-\sqrt{\mu/L})^{k},$ we would expect values of $\rho$ proportional to $\sqrt{\mu/L}$ if the ODE behaves similarly to the discrete process. We have been able to establish this result for both the proximal and AGM ODEs for quadratic objectives (proof in the Appendix in the supplementary material).

Theorem 1**.**

The proximal and AGM ODEs decay with at least the following rates for $\mu$-strongly convex and $L$-smooth quadratic objective functions when using the same hyper-parameters as in the discrete case:

[TABLE]

Figure 2 contrasts the convergence of the discrete and continuous variants. The two methods have quite distinct paths whose shape is shared by their ODE counterparts.

8 Related Work

The application of Bregman divergence to the analysis of continuous time views of the accelerated gradient method has recently been explored by Wibisono et al. (2016) and Wilson et al. (2018). Their approaches do not use the Bregman divergence of $f^{*}$, a key factor of our approach. The Bregman divergence of a function $\phi$ occurs explicitly as a term in a Hamiltonian, in contrast to our view of $\phi$ as curving space. The accelerated gradient method has been shown to be modeled by a momentum of the form ODE $\ddot{X}+c(t)\dot{X}+\nabla f(x)=0$ by Su et al. (2014). Natural discretizations of their ODE result in the heavy-ball method instead of the accelerated gradient method, unlike our form which can produce both based on the choice of $\alpha$. The fine-grained properties of momentum ODEs have also been studied in the quadratic case by Scieur et al. (2017).

A primal-dual form of the regularized accelerated gradient method appears in Lan and Zhou (2017). Our form can be seen as a special case of their form when the regularizer is zero. Our work extends theirs, providing an understanding of the role that geometry plays in unifying acceleration and implicit steps.

The Riemannian connection induced by a function has been heavily explored in the optimization literature as part of the natural gradient method (Amari, 1998), although other connections on this manifold are less explored. The dual-flat connections have primarily seen use in the information-geometry setting for optimization over distributions (Amari and Nagaoka, 2000).

The accelerated gradient method is not the only way to achieve accelerated rates among first order methods. Other techniques include the Geometric descent method of Bubeck et al. (2015), where a bounding ball is updated at each step that encloses two other balls, a very different approach. The method described by Nemirovski and Yudin (1983) is also notable as being closer to the conjugate gradient method than other accelerated approaches, but at the expense of requiring a 2D search instead of a 1D line search at each step.

9 Conclusion

We believe the tools of differential geometry may provide a new and insightful avenue for the analysis of accelerated optimization methods. The analysis we provide in this work is a first step in this direction. The advantage of the differential geometric approach is that it provides high level tools that make the derivation of acceleration easier to state. This derivation, from the proximal point method to the accelerated gradient method, is in our opinion not nearly as mysterious as the other known approaches to understanding acceleration.

Appendix

Appendix A Reformulations of the accelerated gradient method

Form II

This simplification is described in Nesterov [2013] which we reproduce for completeness. Recall that Form I is given by the updates:

[TABLE]

Nesterov specifies the requirement that $\gamma=(1-\alpha)\gamma+\alpha\mu$. When $\alpha=\sqrt{\mu/L}$ then $\gamma$ must then satisfy:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We may rewrite the $y^{k}$ definition as:

[TABLE]

[TABLE]

Plugging this into the $v$ step:

[TABLE]

For $y^{k}$, we start with using $\gamma=\mu$, then apply the $v$ simplification;

[TABLE]

Note that by multiplying by $\sqrt{L}/\mu$ we get:

[TABLE]

Sutskever’s form

Recall Sutskever’s form:

[TABLE]

and Nesterov’s form:

[TABLE]

We will show that using the substitutions:

[TABLE]

[TABLE]

applied to Sutskever’s form gives Nesterov’s form. We start with the momentum term:

[TABLE]

[TABLE]

[TABLE]

Defining: $y_{\text{Nes}}^{k}:=x_{\text{Nes}}^{k+1}+\frac{1}{L}\nabla f\left(y_{\text{Nes}}^{k}\right)$ and applying on the right gives Nesterov’s $y$ update:

[TABLE]

Modern form

We want:

[TABLE]

Starting from Sutskever’s form,

[TABLE]

Define $x_{\text{Mod}}^{k}=x_{\text{Sut}}^{k}+\beta p_{\text{Sut}}^{k}$. Note that this is equal to $y_{\text{Nes}}^{k}$ by definition. So we have:

[TABLE]

Plugging that into the Sutskever step

[TABLE]

[TABLE]

Then define $p_{\text{mod}}^{k}=-Lp_{\text{Sut}}^{k},$ so the update becomes:

[TABLE]

The momentum update changes from:

[TABLE]

to:

[TABLE]

[TABLE]

Auslender & Teboulle form

$y_{\text{AT}}^{k}=(1-\theta_{\text{AT}})\hat{x}_{\text{AT}}^{k}+\theta_{\text{AT}}z_{\text{AT}}^{k},$

$z_{\text{AT}}^{k+1}=z_{\text{AT}}^{k}-\frac{\gamma_{\text{AT}}}{\theta_{\text{AT}}}\nabla f(y_{\text{AT}}^{k}),$

$\hat{x}_{\text{AT}}^{k+1}=(1-\theta_{\text{AT}})\hat{x}_{\text{AT}}^{k}+\theta_{\text{AT}}z_{\text{AT}}^{k+1}.$

We first eliminate $y_{\text{AT}}^{k}$ from the $\hat{x}_{\text{AT}}^{k+1}$ update:

[TABLE]

This notational similarity with Nesterov’s $x_{\text{Nes}}^{k+1}=y_{\text{Nes}}^{k}-\frac{1}{L}\nabla f(y_{\text{Nes}}^{k})$ suggests matching

[TABLE]

as well as $\gamma_{\text{AT}}=\frac{1}{L}$ and

[TABLE]

Note that using using the step for $\hat{x}$ we can rearrange to get:

[TABLE]

Now to determine $\theta$ we simplify using this substitution:

[TABLE]

So my matching constants against Nesterov’s method we have $\beta_{\text{Nes}}=1-\theta_{\text{AT}}.$

Lan form

[TABLE]

Now we can eliminate $\tilde{x}$ from the updates to give

[TABLE]

Also consider the $y^{k}$ update for the AT method:

$y_{\text{AT}}^{k}=(1-\theta_{\text{AT}})\hat{x}_{\text{AT}}^{k}+\theta_{\text{AT}}z_{\text{AT}}^{k},$

We can write this to not involve $\hat{x}$, giving:

[TABLE]

This form suggests matching the iterates with:

[TABLE]

Under this matching, the constants need to satisfy the following relations:

[TABLE]

[TABLE]

[TABLE]

The settings:

[TABLE]

result in these three constraints being satisfied:

[TABLE]

[TABLE]

[TABLE]

Matching the $x_{\text{Lan}}^{k}$ step $x_{\text{Lan}}^{k}=x_{\text{Lan}}^{k-1}-\frac{1}{\eta_{\text{Lan}}}g_{\text{Lan}}^{k}.$ requires $\frac{1}{\eta_{\text{Lan}}}=\frac{\gamma_{\text{AT}}}{\theta_{\text{AT}}}$ also.

Appendix B Continuous time theory

Lemma 2**.**

Consider a linear ODE of the form:

[TABLE]

Where $A$ is real and diagonalizable but not necessarily normal or symmetric, and $x^{*}$ is the fixed point. Then

[TABLE]

where $\lambda_{j}$ are eigenvalues of $A$ with real part $Re[\lambda_{j}]$.

Proof.

Because the inhomogeneous term $Au^{*}$ is constant w.r.t time, without loss of generality we can reduce our problem to a homogenous ODE by shifting the origin:

[TABLE]

The solution of a linear ODE of this form is given in closed form using the matrix exponential:

[TABLE]

we can use this formula to bound the norm of $u(t)$:

[TABLE]

where $\left\|\cdot\right\|$ is the spectral norm. If $A$ is diagonalizable with $A=U\text{diag}(\lambda_{1},\dots\lambda,d)U^{-1}$ then the matrix exponential can be expressed as:

[TABLE]

The spectral norm is given by the largest absolute value of the eigenvalues, so we must consider the interaction of the real and complex parts. For an eigenvalue $\lambda=a+bi$ of $A$, the norm takes the simple form:

[TABLE]

So the spectral norm is given by the maximum over the eigenvalues $\lambda_{j}$ of $\exp\left(tRe[\lambda_{j}]\right)$. ∎

Theorem 3**.**

Consider the following linear ODE:

[TABLE]

[TABLE]

where $H:n\times n$ is a real, positive definite and symmetric matrix with minimum eigenvalue $\mu$ and maximum eigenvalue $L$. This corresponds to the proximal ODE for a quadratic function $f(x)=\frac{1}{2}\left(x-x^{*}\right)^{T}H\left(x-x^{*}\right)$, with $u=[x;g]$ and $u^{*}=[x^{*};0]$. Then the decay rate of the ODE towards the origin $u^{*}=0$ can be bounded as follows for $\eta=\sqrt{\mu L}$:

[TABLE]

[TABLE]

Proof.

We will take the approach of bounding the real parts of the eigenvalues of $A$, so that we can directly apply Lemma 2.

Let $U\Lambda U^{T}=H$ be the eigen-decomposition of $H$. Note that the operation of conjugation by $U$ leaves the identity matrix unchanged ($UIU^{T}=I$). Each block of A is just a weighted combination of the identity matrix and $H^{-1}$, so this implies that conjugation of a block by it’s self gives a diagonal matrix. We can use this idea to define a similarity transform that converts $A$ into a matrix where each of the four blocks are diagonal. In particular we have:

[TABLE]

where we have written A in terms of its 4 constituent $n\times n$ blocks. Each block is a weighted sum of the identity matrix and the diagonal matrix of inverse eigenvalues of H, for instance:

[TABLE]

Next we construct a permutation matrix $\Pi:2n\times 2n$ with the goal of converting $D$ into a block diagonal matrix with $2\times 2$ blocks along the diagonal, where each block has the structure of $A$ as if it was applied to a 1D optimization problem, with $H$ being replaced by one of the $d$ eigenvalues of $H$. This is achieved with the permutation matrix $\Pi$ that is zero except for:

[TABLE]

For instance, in the $n=2$ case the matrix is:

[TABLE]

This matrix has the effect of interleaving the primal dual pairs (per coordinate) instead of having all the primal coordinates together followed by all the dual coordinates. So when we conjugate using $\Pi$ we get:

[TABLE]

where each $T$ is a $2\times 2$ matrix of the described form:

[TABLE]

The eigenvalues of a block diagonal matrix are just the eigenvalues of the blocks concatenated, and since there is a similarity transform between $A$ and this block diagonal matrix, we have effectively reduced our problem to considering 1D quadratics, with curvature between $\mu$ and $L$, for fixed $\eta$.

Recall that for a matrix \left[\begin{array}[]{cc}a&b\\ c&d\end{array}\right] the eigenvalues are given by the two roots of a quadratic, namely:

[TABLE]

We use the notation $\nu$ to avoid confusion between the eigenvalues of the $T$ blocks and those of H . For a block $T_{i}$, this expression is

[TABLE]

Suppose that the discriminate (the quantity under the square root) is negative, then

[TABLE]

this is obviously at least as small as $-\rho=-\frac{\sqrt{\mu}}{\sqrt{\mu}+\sqrt{L}},$ since the largest value of $\rho$ possible is when $\mu=L$, in which case $\rho=\frac{1}{2}$. So consider instead the case where the discriminate is positive. We need only consider the $v_{+}$ root as it is strictly larger. Then we will use the concavity of the square root function to bound $\nu_{+}$,

[TABLE]

For $h=\sqrt{\cdot},$ with $y=\left(-1-\eta\lambda_{i}^{-1}\right)^{2}$ and $x=\left(-1-\eta\lambda_{i}^{-1}\right)^{2}-4$. We get:

[TABLE]

Therefore:

[TABLE]

Thus we have shown that the real parts of all eigenvalues of $A$ are less than $-\rho$. ∎

Theorem 4**.**

Consider the following linear ODE:

[TABLE]

[TABLE]

where $H:d\times d$ is a real, positive definite and symmetric matrix with minimum eigenvalue $\mu$ and maximum eigenvalue $L$. This corresponds to the AGM ODE for a quadratic objective$f(x)=\frac{1}{2}\left(x-x^{*}\right)^{T}H\left(x-x^{*}\right)$, with $u=[x;g]$ and $u^{*}=[x^{*};0]$. The decay rate of this ODE towards the origin $u^{*}=0$ can be bounded as follows for $\eta=\sqrt{\mu L}$, $\tau=L/\eta$, and $\alpha\in[0,1]$:

[TABLE]

[TABLE]

Proof.

We can apply the same proof technique as for Theorem 3, we omit the details. The $2\times 2$ block diagonal matrices are:

[TABLE]

For eigenvalues $\lambda_{i}$ of $H$. The two eigenvalues of $T_{i}$ are:

[TABLE]

The discriminate here is always non-positive. To see why, we plug in the constants $\eta$ and $\tau$:

[TABLE]

So the real part is given by the quantity outside the square root, namely:

[TABLE]

Using Lemma 2 gives the result. ∎

Appendix C The standard heavy ball ODE

The standard heavy-ball ODE for a quadratic $f(x)=\frac{1}{2}x^{T}Hx$ is

[TABLE]

Which can be written in first-order form in terms of a momentum parameter $p$ as:

[TABLE]

The constants that result in optimal convergence rate for the discretize heavy ball method:

[TABLE]

which can also be written as:

[TABLE]

are:

[TABLE]

Theorem 5**.**

Consider the following ODE:

[TABLE]

[TABLE]

This is the ODE in Equation 8 written in matrix form for a combined iterate $u$. For the parameters given in Equation 9, this ODE has decay rate at least:

[TABLE]

[TABLE]

Proof.

We can reduce the problem to considering the eigenvalues of $2\times 2$ matrices as we did for the proximal and AGM ODEs. We have matrices of the form:

[TABLE]

whose eigenvalues are given by the general formula

[TABLE]

Simplifying:

[TABLE]

Note that:

[TABLE]

The choice of step size ensures that the discriminant is always non-positive:

[TABLE]

Therefore the decay rate is bounded by the real part of the eigenvalues, which is $-\frac{1}{2}(1-\beta).$ ∎