A Variational Perspective on Accelerated Methods in Optimization

Introduction

The phenomenon of acceleration plays an important role in theory and practice of convex optimization. Introduced by Nesterov in 1983 in the context of gradient descent, the acceleration idea has been extended to a wide range of other settings, including composite optimization , stochastic optimization , nonconvex optimization , and conic programming . There have been generalizations to non-Euclidean optimization and higher-order algorithms , and there have been numerous applications that further extend the reach of the idea . On the theoretical front, acceleration often improves the convergence rate of the underlying gradient-based procedure, and, under certain conditions, yields an optimal rate .

Despite this compelling evidence of the value of the idea of acceleration, it remains something of a conceptual mystery. Derivations of accelerated methods do not flow from a single underlying principle, but tend to rely on case-specific algebra . The basic Nesterov technique is often explained intuitively in terms of momentum, but this intuition does not easily carry over to non-Euclidean settings . In recent years, the number of explanations and interpretations of acceleration has increased , but these explanations have been focused on restrictive instances of acceleration, such as first-order algorithms, the Euclidean setting, or cases in which the objective function is strongly convex or quadratic. It is not yet clear what the natural scope of the acceleration concept is and indeed whether it is a single phenomenon.

In this paper we study acceleration from a continuous-time, variational point of view. We build on recent work by , who show that the continuous-time limit of Nesterov’s accelerated gradient descent is a second-order differential equation, and we take inspiration from continuous time analysis of mirror descent . In our approach, rather than starting from existing discrete-time accelerated gradient methods and deriving differential equations by taking limits, we take as our point of departure a variational formulation in which we define a functional on continuous-time curves that we refer to as a Bregman Lagrangian. Next, we calculate and discretize the Euler-Lagrange equation corresponding to the Bregman Lagrangian. It turns out that naive discretization (the Euler method) does not yield a stable discrete-time algorithm that retains the rate of the underlying differential equation; rather, a more elaborate discretization involving an auxiliary sequence is necessary. This auxiliary sequence is essentially that used by Nesterov in his constructions of accelerated mirror descent and accelerated cubic-regularized Newton’s method , and later generalized by Baes . Thus, from our perspective, Nesterov’s approach can be viewed as a methodology for the discretization of a certain class of differential equations. Given the complexities associated with the discretization of differential equations, it is perhaps not surprising that it has been difficult to perceive the generality and scope of the acceleration concept in a discrete-time framework.

The Bregman-Lagrangian framework permits a systematic understanding of the matching rates associated with higher-order gradient methods in discrete and continuous time. In the case of gradient descent, Su et al. show that the discrete and continuous-time dynamics have convergence rates of $O(1/(\epsilon k))$ and $O(1/t)$, respectively, and that these match using the identification $t=\epsilon k$; for accelerated gradient descent, the convergence rates are $O(1/(\epsilon k^{2}))$ and $O(1/t^{2})$ respectively, which match using the identification $t=\sqrt{\epsilon}k$ . This result has been extended to the non-Euclidean case by Krichene et al. . Higher-order gradient descent is a descent method which minimizes a regularized $(p-1)$-st order Taylor approximation of the objective function $f$, generalizing gradient descent ($p=2$) and Nesterov and Polyak’s cubic-regularized Newton’s method ($p=3$) . The $p$-th order gradient algorithm with a constant step size $\epsilon$ has convergence rate $O(1/(\epsilon k^{p-1}))$ when $\nabla^{p-1}f$ is $(1/\epsilon)$-Lipchitz and, in continuous time, as $\epsilon\to 0$, this algorithm corresponds to the $p$-th rescaled gradient flow, which is a first-order differential equation with a matching convergence rate $O(1/t^{p-1})$. Thus, the $p$-th order gradient algorithm can be seen as a discretization $t=\delta k$ of the rescaled gradient flow with time step $\delta=\epsilon^{1/(p-1)}$. Similarly, we show that the accelerated higher-order gradient algorithm achieves an improved convergence rate $O(1/(\epsilon k^{p}))$ under the same assumption (i.e., $\nabla^{p-1}f$ is $(1/\epsilon)$-Lipschitz). In continuous time, as $\epsilon\to 0$, this corresponds to the second-order Euler-Lagrange curve of the Bregman Lagrangian with a matching convergence rate $O(1/t^{p})$. Thus, the $p$-th order accelerated algorithm can be seen as a discretization $t=\delta k$ of the Euler-Lagrange equation of the Bregman Lagrangian with time step $\delta=\epsilon^{1/p}$.

In addition to its value in relating continuous-time and discrete-time acceleration, the study of the Bregman Lagrangian can provide further insights into the nature of acceleration. For instance, it is noteworthy that the Bregman Lagrangian is closed under time dilation. This means that if we take an Euler-Lagrange curve of a Bregman Lagrangian and reparameterize time so we travel the curve at a different speed, then the resulting curve is also the Euler-Lagrange curve of another Bregman Lagrangian, with appropriately modified parameters. Thus, the entire family of accelerated methods correspond to a single curve in spacetime and can be obtained by speeding up (or slowing down) any single curve. Another insight is obtained by noting that from the discrete-time point of view, an interpretation of acceleration starts with a base algorithm, which we can accelerate by coupling with a suitably weighted mirror descent step. From the continuous-time point of view, however, it is the weighted mirror descent step that is important since the base gradient algorithm operates on a smaller time scale. Thus, Nesterov’s accelerated gradient methods are but one possible implementation of second-order Bregman-Lagrangian curves as a discrete-time algorithm.

The remainder of the paper is organized as follows. In Section 2, we introduce the general family of Bregman Lagrangians and study its properties. In Section 3, we demonstrate how to discretize the Euler-Lagrange equations corresponding to the polynomial subfamily of Bregman Lagrangians to obtain discrete-time accelerated algorithms. In particular, we introduce the family of higher-order gradient methods which can be used to complete the discretization. In Section 4, we discuss additional properties of the Bregman Lagrangian, including gauge-invariance properties, connection to classical gradient flows, and the correspondence with a functional that we refer to as a Bregman Hamiltonian. Finally, we end in Section 5 with a brief discussion.

which is nonnegative since $h$ is convex. When $x$ and $y$ are nearby the Bregman divergence is an approximation to the Hessian metric,

The Euclidean setting is obtained when $h(x)=\frac{1}{2}\|x\|^{2}$, in which case the Bregman divergence and Hessian metric coincide since $\nabla^{2}h(x)$ is the identity matrix.

In continuous time, the Hessian metric is generally studied rather than the more general Bregman divergence; this is the case, for instance, in the case of natural gradient flow, which is the continuous-time limit of mirror descent . By way of contrast, we shall see that our continuous-time, Lagrangian framework crucially employs the Bregman divergence.

The Bregman Lagrangian

these conditions will be justified in the following section.

In this section we show that—under the ideal scaling assumption (2.2)—the Bregman Lagrangian (2.1) defines a variational problem the solutions to which minimize the objective function $f$ at an exponential rate.

Specifically, for the Bregman Lagrangian (2.1), the partial derivatives are

Thus, for general functions $\alpha_{t},\beta_{t},\gamma_{t}$, the Euler-Lagrange equation (2.3) for the Bregman Lagrangian (2.1) is a second-order differential equation given by

We now impose the ideal scaling condition (2.2b). In this case the last term in (2.5) vanishes, so the Euler-Lagrange equation simplifies to

In (2.6), we have assumed the Hessian matrix $\nabla^{2}h(X_{t}+e^{-\alpha_{t}}\dot{X}_{t})$ is invertible. But we can also write the equation (2.6) in the following way, which only requires that $\nabla h$ be differentiable,

To establish a convergence rate associated with solutions to the Euler-Lagrange equation—under the ideal scaling conditions—we take a Lyapunov function approach. Defining the following energy functional:

we immediately obtain a convergence rate, as shown in the following theorem.

If the ideal scaling (2.2) holds, then solutions to the Euler-Lagrange equation (2.7) satisfy

The time derivative of the energy functional is

If $X_{t}$ satisfies the Euler-Lagrange equation (2.7), then the time derivative simplifies to

For a given $\alpha_{t}$, which determines $\gamma_{t}$ by (2.2a), the optimal convergence rate is achieved by setting $\dot{\beta}_{t}=e^{\alpha_{t}}$, resulting in convergence rate $O(e^{-\beta_{t}})=O(\exp(-\int_{t_{0}}^{t}e^{\alpha_{s}}\,ds))$. In Section 3 we study a subfamily of Bregman Lagrangians that have a polynomial convergence rate, and we show how we can discretize the resulting Euler-Lagrange equations to obtain discrete-time methods that have a matching, accelerated convergence rate. In Section 4 we study another subfamily of Bregman Lagrangians that have an exponential convergence rate, and discuss its connection to a generalization of Nesterov’s restart scheme. In the Euclidean setting, our derivations simplify. We present these derivations in Appendix B.7, and comment on the insight that they provide into the question posed by Su et al. on the significance of the value $3$ in the damping coefficient for Nesterov’s accelerated gradient descent.

2 Time dilation

A notable property of the Bregman Lagrangian family is that it is closed under time dilation. This means if we take the Euler-Lagrange equation (2.5) of the Bregman Lagrangian (2.1) and reparameterize time to travel the curve at a different speed, the resulting curve is also the Euler-Lagrange equation of a Bregman Lagrangian with a suitably modified set of parameters.

That is, the new curve $Y$ is obtained by traversing the original curve $X$ at a new speed of time determined by $\tau$. If $\tau(t)>t$, then we say that $Y$ is the sped-up version of $X$, because the curve $Y$ at time $t$ has the same value as the original curve $X$ at the future time $\tau(t)$.

For clarity, we let $\mathcal{L}_{\alpha,\beta,\gamma}$ denote the Bregman Lagrangian (2.1) parameterized by $\alpha,\beta,\gamma$. Then we have the following result whose proof is provided in Appendix A.1.

In Section 3, we show that the Bregman Lagrangian generates the family of higher-order accelerated methods in discrete time. Thus, the time-dilation property means that the entire family of curves for accelerated methods in continuous time corresponds to a single curve in spacetime, which is traveled at different speeds. This suggests that the underlying solution curve has a more fundamental structure that is worth exploring further.

Polynomial convergence rates and accelerated methods

In this section, we study a subfamily of Bregman Lagrangians (2.1) with the following choice of parameters, indexed by a parameter $p>0$,

where $C>0$ is a constant. The parameters $\alpha,\beta,\gamma$ satisfy the ideal scaling condition (2.2) (with an equality on the first condition (2.2a)). The Euler-Lagrange equation (2.6) is given by

and, by Theorem 2.1, it has an $O(1/t^{p})$ rate of convergence. As direct result of the time-dilation property (Theorem 2.2), the entire family of curves (3.2) can be obtained by speeding up the curve in the case $p=2$ by the time-dilation function $\tau(t)=t^{p/2}$. In Appendix A.2 we discuss the issue of the existence and uniqueness of the solution to the differential equation (3.2).

The case $p=2$ of the equation (3.2) is the continuous-time limit of Nesterov’s accelerated mirror descent , and the case $p=3$ is the continuous-time limit of Nesterov’s accelerated cubic-regularized Newton’s method . The case $p=2$ has also been derived independently in a recent work of Krichene et al. ; in the Euclidean case, when the Hessian $\nabla^{2}h$ is the identity matrix, we recover the differential equation of Su et al. .

We now turn to the challenge of discretizing the differential equation in (3.2), with the goal of obtaining a discrete-time algorithm whose convergence rate matches that of the underlying differential equation. As we show in this section, a naive Euler method is not able to match the underlying rate. To match the rate a more sophisticated approach is needed, and it is at this juncture that Nesterov’s three-sequence idea makes its appearance.

We first write the second-order equation (3.2) as the following system of first-order equations:

Now we discretize $X_{t}$ and $Z_{t}$ into sequences $x_{k}$ and $z_{k}$ with time step $\delta>0$. That is, we make the identification $t=\delta k$ and set $x_{k}=X_{t}$, $x_{k+1}=X_{t+\delta}\approx X_{t}+\delta\dot{X}_{t}$ and $z_{k}=Z_{t}$, $z_{k+1}=Z_{t+\delta}\approx Z_{t}+\delta\dot{Z}_{t}$. Applying the forward-Euler method to (3.3a) gives the equation $z_{k}=x_{k}+\frac{\delta k}{p}\frac{1}{\delta}(x_{k+1}-x_{k})$, or equivalently,

Similarly, applying the backward-Euler method to equation (3.3b) gives $\frac{1}{\delta}(\nabla h(z_{k})-\nabla h(z_{k-1}))=-Cp(\delta k)^{p-1}\nabla f(x_{k})$, which we can write as the optimality condition of the following weighted mirror descent step:

with step size $\epsilon=\delta^{p}$. In principle, the two updates (3.4), (3.5) define an algorithm that implements the dynamics (3.3) in discrete time. However, we cannot establish a convergence rate for the algorithm (3.4), (3.5); indeed, empirically, we find that the algorithm is unstable. Even for the simple case in which $f$ is a quadratic function in two dimensions, the iterates of the algorithm initially approach and oscillate near the minimizer, but eventually the oscillation increases and the iterates shoot off to infinity.

2 A rate-matching discretization

We now discuss how to modify the naive discretization scheme (3.4), (3.5) into an algorithm whose rate matches that of the underlying differential equation. Our approach is inspired by Nesterov’s constructions of accelerated mirror descent and accelerated cubic-regularized Newton’s method , which maintain three sequences in the algorithms and use the estimate sequence technique to prove convergence. Indeed, from our point of view, Nesterov’s methodology can be viewed as a rate-matching discretization methodology.

Specifically, we consider the following scheme, in which we introduce a third sequence $y_{k}$ to replace $x_{k}$ in the updates,

where $k^{(p-1)}:=k(k+1)\cdots(k+p-2)$ is the rising factorial. A sufficient condition for the algorithm (3.6) to have an $O(1/(\epsilon k^{p}))$ convergence rate is that the new sequence $y_{k}$ satisfy the inequality

for some constant $M>0$. Note that in going from (3.4) to (3.6a) we have replaced the weight $\frac{p}{k}$ by $\frac{p}{k+p}$; this is only for convenience in the proof given below, and does not change the asymptotics since $\frac{p}{k}=\Theta(\frac{p}{k+p})$ as $k\to\infty$. Similarly, we replace $k^{p-1}$ in (3.5) by the rising factorial $k^{(p-1)}$ in (3.6b) to make the algebra easier, but we still have $k^{(p-1)}=\Theta(k^{p-1})$.

The following result also requires a uniform convexity assumption on the distance-generating function $h$. Recall that $h$ is $\sigma$-uniformly convex of order $p\geq 2$ if its Bregman divergence is lower bounded by the $p$-th power of the norm,

The case $p=2$ is the usual definition of strong convexity. An example of a uniformly convex function is the $p$-th power of the norm, $h(x)=\frac{1}{p}\|x-w\|^{p}$ for any $w\in\mathcal{X}$, which is $\sigma$-uniformly convex of order $p$ with $\sigma=2^{-p+2}$ [28, Lemma 4].

Assume $h$ is $1$-uniformly convex of order $p\geq 2$, and the sequence $y_{k}$ satisfies the inequality (3.7) for all $k\geq 0$. Then the algorithm (3.6) with the constant $C\leq M^{p-1}/p^{p}$ and initial condition $z_{0}=x_{0}\in\mathcal{X}$ has the convergence rate

The proof of Theorem 3.1 uses a generalization of Nesterov’s estimate sequence technique, and can be found in Appendix A.3. We note that with the scaling $\epsilon=\delta^{p}$ as in the previous section, the convergence rate $O(1/(\epsilon k^{p}))$ matches the $O(1/t^{p})$ rate in continuous time for the differential equation (3.2). We also note that the result in Theorem 3.1 does not require any assumptions on $f$ beyond the ability to construct a sequence $y_{k}$ satisfying (3.7). In the next section, we will see that we can satisfy (3.7) using the higher-order gradient method, which requires a higher-order smoothness assumption on $f$; the resulting algorithm is then the accelerated higher-order gradient method.

3 Higher-order gradient method

We study the higher-order gradient update, which minimizes a regularized higher-order Taylor approximation of the objective function $f$.

Recall that for an integer $p\geq 2$, the $(p-1)$-st order Taylor approximation of $f$ centered at $x\in\mathcal{X}$ is the $(p-1)$-st degree polynomial

We say that $f$ is $L$-smooth of order $p-1$ if $f$ is $p$-times continuously differentiable and $\nabla^{p-1}f$ is $L$-Lipschitz, which means for all $x,y\in\mathcal{X}$,

For a constant $N>0$ and step size $\epsilon>0$, we define the update operator $G_{p,\epsilon,N}\colon\mathcal{X}\to\mathcal{X}$ by

When $f$ is smooth of order $p-1$, the operator $G_{p,\epsilon,N}$ has the following property, which generalizes [28, Lemma 6]. We provide the proof in Appendix A.4.

Let $x\in\mathcal{X}$, $y=G_{p,\epsilon,N}(x)$, and $N>1$. If $f$ is $L=\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, then

The inequality (3.12) means that we can use the update operator $G_{p,\epsilon,N}$ to produce a sequence $y_{k}$ satisfying the requirement (3.7) under a higher-order order smoothness condition on $f$. We state the resulting algorithm in the next section.

In this section, we study the following higher-order gradient algorithm defined by the update operator $G_{p,\epsilon,N}$:

The case $p=2$ is the usual gradient descent algorithm, and the case $p=3$ is Nesterov and Polyak’s cubic-regularized Newton’s method .

If $f$ is smooth of order $p-1$, then the algorithm (3.14) is a descent method. Furthermore, we can prove the following rate of convergence, which generalizes the results for gradient descent and the cubic-regularized Newton’s method. We provide the proof in Appendix A.5.

If $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, then the algorithm (3.14) with constant $N>0$ and initial condition $x_{0}\in\mathcal{X}$ has the convergence rate

where $R=\sup_{x\colon f(x)\leq f(x_{0})}\|x-x^{*}\|$ is the radius of the level set of $f$ from the initial point $x_{0}$.

We can take the continuous-time limit of the higher-order gradient algorithm as the step size $\epsilon\to 0$. The resulting curve is a first-order differential equation that is a rescaled version of gradient flow. We show that it minimizes $f$ with a matching convergence rate. In the following, we take $N=1$ in (3.14) for simplicity (the general $N$ simply scales the vector field by a constant). We provide the proof of Theorem 3.4 in Appendix A.6.

The continuous-time limit of the algorithm (3.14) is the rescaled gradient flow

where we define the right-hand side to be the zero if $\nabla f(X_{t})=0$. Furthermore, the rescaled gradient flow has convergence rate

where $R=\sup_{x\colon f(x)\leq f(X_{0})}\|x-x^{*}\|$ is the radius of the level set of $f$ from the initial point $X_{0}$.

Equivalently, we can interpret the higher-order gradient algorithm (3.14) as a discretization of the rescaled gradient flow (3.16) with time step $\delta=\epsilon^{\frac{1}{p-1}}$, so $t=\delta k=\epsilon^{\frac{1}{p-1}}k$. With this identification, the convergence rates in discrete time, $O(1/(\epsilon k^{p-1}))$, and in continuous time, $O(1/t^{p-1})$, match. The convergence rate for the continuous-time dynamics does not require any assumption beyond the convexity and differentiability of $f$ (as in the case of the Lagrangian flow (2.6)), whereas the convergence rate for the discrete-time algorithm requires the higher-order smoothness assumption on $f$. We note that the limiting case $p\to\infty$ of (3.16) is the normalized gradient flow, which has been shown to converge to the minimizer of $f$ in finite time . We also note that unlike the Lagrangian flow, the family of rescaled gradient flows is not closed under time dilation.

4 Accelerated higher-order gradient method

By the result of Lemma 3.2, we see that we can use the higher-order gradient update $G_{p,\epsilon,N}$ to produce a sequence $y_{k}$ satisfying the inequality (3.7), to complete the algorithm (3.14) that implements the polynomial family of the Bregman-Lagrangian flow (3.2). Explicitly, the resulting algorithm is as follows,

By Theorem 3.1 and Lemma 3.2, we have the following guarantee for this algorithm.

Assume $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, and $h$ is $1$-uniformly convex of order $p$. Then the algorithm (3.18) with constants $N>1$ and $C\leq(N^{2}-1)^{\frac{p-2}{2}}/((2N)^{p-1}p^{p})$ and initial conditions $z_{0}=x_{0}\in\mathcal{X}$ has an $O(1/(\epsilon k^{p}))$ convergence rate.

The resulting algorithm (3.18) and its convergence rate recovers the results of Baes , who studied a generalization of Nesterov’s estimate sequence technique to higher-order algorithms. We note that the convergence rate $O(1/(\epsilon k^{p}))$ of algorithm (3.18) is better than the $O(1/(\epsilon k^{p-1}))$ rate of the higher-order gradient algorithm (3.14), under the same assumption of the $(p-1)$-st order smoothness of $f$. This gives the interpretation of the algorithm (3.18) as “accelerating” the higher-order gradient method. Indeed, in this view the “base algorithm” that we start with is the higher-order gradient algorithm in the $y$-sequence (3.18b), and the acceleration is obtained by coupling it with a suitably weighted mirror descent step in (3.18a) and (3.18c).

However, from the continuous-time point of view, where our starting point is the polynomial Lagrangian flow (3.2), we see that the algorithm (3.18) is only one possible implementation of the flow as a discrete-time algorithm. As we saw in Section 3.2, it is only the $x$- and $z$-sequences (3.18a) and (3.18c) that play a role in the correspondence between the continuous-time dynamics and its discrete-time implementation, and the requirement (3.7) in the $y$-update is only needed to complete the convergence proof. Indeed, the higher-order gradient update (3.18b) does not change the continuous-time limit, since from (3.13) in Lemma 3.2 we have that $\|x_{k}-y_{k}\|=\Theta(\epsilon^{\frac{1}{p-1}})$, which is smaller than the $\delta=\epsilon^{\frac{1}{p}}$ time step in the discretization of (3.2). Therefore, the $x$ and $y$ sequences in (3.18) coincide in continuous time as $\epsilon\to 0$. Thus, from this point of view, Nesterov’s accelerated methods (for the cases $p=2$ and $p=3$) are one of possibly many discretizations of the polynomial Lagrangian flow (3.2). For instance, in the case $p=2$, Krichene et al. [17, Section 4.1] show that we can use a general regularizer in the gradient step (3.18b) under some additional smoothness assumptions. If there are other implementations, it would be interesting to see if the higher-gradient methods have some distinguishing property, such as computational efficiency.

Further explorations of the Bregman Lagrangian

In addition to providing a unifying framework for the generation of accelerated gradient-based algorithms, the Bregman Lagrangian has mathematical structure that can be investigated directly. In this section we briefly discuss some of the additional perspective that can be obtained from the Bregman Lagrangian. See Appendices B.1–B.6 for technical details of the results discussed here.

It is important to note the presence of the Bregman divergence in the Bregman Lagrangian (2.1). In the non-Euclidean setting, intuition might suggest using the Hessian metric $\nabla^{2}h$ to measure a “kinetic energy,” and thereby obtain a Hessian Lagrangian. This approach turns out to be unsatisfying, however, because the resulting differential equation does not yield a convergence rate and the Euler-Lagrange equation involves the third-order derivative $\nabla^{3}h$, posing serious difficulties for discretization. As we have seen, the Bregman Lagrangian, on the other hand, readily provides a rate of convergence via a Lyapunov function; moreover, the resulting discrete-time algorithm in (3.18) involves only the gradient $\nabla h$ via the weighted mirror descent update.

In the Euclidean case, it is known classically that we can view gradient flow as the strong-friction limit of a damped Lagrangian flow [34, p. 646]. We show that the same interpretation holds for natural gradient flow and rescaled gradient flow. In particular, we show in Appendix B.3 that we can recover natural gradient flow as the strong-friction limit of a Bregman Lagrangian flow with an appropriate choice of parameters. Similarly, we can recover the rescaled gradient flow (3.16) as the strong-friction limit of a Lagrangian flow that uses the $p$-th power of the norm as the kinetic energy. Therefore, the general family of second-order Lagrangian flows is more general, and includes first-order gradient flows in its closure. From this point of view, a particle with gradient-flow dynamics is operating in the regime of high friction. The particle simply rolls downhill and stops at the equilibrium point as soon as the force $-\nabla f$ vanishes; there is no oscillation since it is damped by the infinitely strong friction. Thus, the effect of moving from a first-order gradient flow to a second-order Lagrangian flow is to reduce the friction from infinity to a finite amount; this permits oscillation , but also allows faster convergence.

One way to understand a Lagrangian is to study its Hamiltonian, which is the Legendre conjugate (dual function) of the Lagrangian. Typically, when the Lagrangian takes the form of the difference between kinetic and potential energy, the Hamiltonian is the sum of the kinetic and potential energy. The Hamiltonian is often easier to study than the Lagrangian, since its second-order Euler-Lagrangian equation is transformed into a pair of first-order equations. In our case, the Hamiltonian corresponding to the Bregman Lagrangian (2.1) is the following Bregman Hamiltonian,

which indeed has the form of the sum of the kinetic and potential energy. Here the kinetic energy is measured using the Bregman divergence of $h^{*}$, which is the convex dual function of $h$. See Appendix B.4 for further discussion.

The Euler-Lagrange equation of a Lagrangian is gauge-invariant, which means it does not change when we add a total time derivative to the Lagrangian. For the Bregman Lagrangian with the ideal scaling condition (2.2b), this property implies that we can replace the Bregman divergence $D_{h}(X+e^{-\alpha_{t}}V,X)$ in (2.1) by its first term $h(X+e^{-\alpha_{t}}V)$. This might suggest a different interpretation of the role of $h$ in the Lagrangian.

The natural motion of the Bregman Lagrangian (i.e., the motion when there is no force, $-\nabla f\equiv 0$) is given by $X_{t}=ae^{-\gamma_{t}}+b$, for some constants $a,b\in\mathcal{X}$. Notice that even though the Bregman Lagrangian still involves the distance-generating function $h$, its natural motion is actually independent of $h$. Thus, the effect of $h$ is felt only via its interaction with $f$—this can also be seen in (2.6) where $h$ and $f$ only appear together in the final term. Furthermore, assuming $e^{\gamma_{t}}\to\infty$, the natural motion always converges to a limit point, which a priori can be anything. However, as we see from Theorem 2.1, as soon as we introduce a convex potential function $f$, all motions converge to the minimizer $x^{*}$ of $f$.

In addition to the polynomial family in Section 3, we can also study the subfamily of Bregman Lagrangians that have exponential convergence rates $O(e^{-ct})$, $c>0$. As we discuss in Appendix B.1, in this case the link to discrete-time algorithms is not as clear. Using the same discretization technique as in Section 3 suggests that to get a matching convergence rate, constant progress is needed at each iteration.

From the discrete-time perspective, we show that the higher-order gradient algorithm (3.14) achieves an exponential convergence rate when the objective function $f$ is uniformly convex. Furthermore, we show that a restart scheme applied to the accelerated method (3.18) achieves a better dependence on the condition number; this generalizes Nesterov’s restart scheme for the case $p=3$ [28, Section 5].

It is an open question to understand if there is a better connection between the discrete-time restart algorithms and the continuous-time exponential Lagrangian flows. In particular, it is of interest to consider whether a restart scheme is necessary to achieve exponential convergence in discrete time; we know it is not needed for the special case $p=2$, since a variant of Nesterov’s accelerated gradient descent that incorporates the condition number also achieves the optimal convergence rate.

Discussion

In this paper, we have presented a variational framework for understanding accelerated methods from a continuous-time perspective. We presented the general family of Bregman Lagrangian, which generates a family of second-order Lagrangian dynamics that minimize the objective function at an accelerated rate compared to gradient flows. These dynamics are related to each other by the operation of speeding up time, because the Bregman Lagrangian family is closed under time dilation. In the polynomial case, we showed how to discretize the second-order Lagrangian dynamics to obtain an accelerated algorithm with a matching convergence rate. The resulting algorithm accelerates a base algorithm by coupling it with a weighted mirror descent step. An example of a base algorithm is a higher-order gradient method, which in continuous time corresponds to a first-order rescaled gradient flow with a matching convergence rate. Our continuous-time perspective makes clear that it is the mirror descent coupling that is more important for the acceleration phenomenon rather than the base algorithm. Indeed, the higher-order gradient algorithm operates on a smaller timescale than the enveloping mirror descent coupling step, so it makes no contribution in the continuous-time limit, and in principle we can use other base algorithms.

Our work raises many questions for further research. First, the case $p=2$ is worthy of further investigation. In particular, the assumptions needed to show convergence of the discrete-time algorithm ($\nabla^{p-1}f$ is Lipschitz) are different than those required to show existence and uniqueness of solutions of the continuous-time dynamics ($\nabla f$ is Lipschitz). In the case $p=2$ however, these assumptions match. This suggests a strong link between the discrete- and continuous-time dynamics that might help us understand why several results seem to be unique to the special case $p=2$. Second, in discrete time, Nesterov’s accelerated methods have been extended to various settings, for example to the stochastic setting. An immediate question is whether we can extend our Lagrangian framework to these settings. Third, we would like to understand better the transition from continuous-time dynamics to discrete-time algorithms, and whether we can establish general assumptions that preserve desirable properties (e.g., convergence rate). In Section 3 we saw that the polynomial convergence rate requires a higher-order smoothness assumption in discrete time, and in Section 4 we discussed whether the exponential case requires a uniform convexity assumption. Finally, our work to date focuses on the convergence rates of the function values rather than the iterates. Recently there has been some work extending to study the convergence of the iterates and some perturbative aspects ; it would be interesting to extend these results to the general Bregman Lagrangian.

At an abstract level, the general family of Bregman Lagrangian has a rich mathematical structure that deserves further study; we discussed some of these properties in Section 4. We hope that doing so will give us new insights into the nature of the optimization problem in continuous time, and help us design better dynamics with matching discrete-time algorithms. For example, we can study how to use some of the appealing properties of the Hamiltonian formalism (e.g., volume preservation in phase space) to help us discretize the dynamics. We also wish to understand where the Bregman Lagrangian itself comes from, why it works so well, and whether there are other Lagrangian families with similarly favorable properties.

References

Appendix A Proofs of results

The velocity and acceleration of the reparameterized curve $Y_{t}=X_{\tau(t)}$ are given by

By assumption, the original curve $X_{t}$ satisfies the Euler-Lagrange equation (2.5) for the Bregman Lagrangian $\mathcal{L}_{\alpha,\beta,\gamma}$. At time $\tau(t)$, this equation reads

We now use the relations (A.1). After multiplying by $\dot{\tau}(t)^{2}$ and collecting terms, we get

Furthermore, suppose $\alpha,\beta,\gamma$ satisfy the ideal scaling (2.2). Then

A.2 Existence and uniqueness of solution to the polynomial family

In this section we discuss the existence and uniqueness of solution to the differential equation (3.2) arising from the polynomial family of Bregman Lagrangian. We begin by writing the second-order equation (3.2) as the pair of first-order equations (3.3). We also write $W_{t}=\nabla h(Z_{t})$, so we can write (3.3) as

where $\mathcal{X}^{*}$ is the dual space of $\mathcal{X}$, i.e., the space of all linear functionals over $\mathcal{X}$. Under the assumption that $h$ be essentially smooth, the supremum in (A.3) is achieved by $z=\nabla h^{*}(w)$, and we have the relation that $\nabla h$ and $\nabla h^{*}$ are inverses of each other, i.e., $z=\nabla h^{*}(w)\Leftrightarrow w=\nabla h(z)$. Thus, with the definition $W_{t}=\nabla h(Z_{t})$, we can write $Z_{t}=\nabla h^{*}(W_{t})$, which gives us (A.2).

Now assume $\nabla f$ and $\nabla h^{*}$ are Lipschitz continuous functions. Then over any bounded time intervals $[t_{0},t_{1}]$ with $0<t_{0}<t_{1}$, the right-hand side of (A.2) is a Lipschitz continuous vector field. Thus, by the Cauchy-Lipschitz theorem, for any given initial conditions $(X_{t_{0}},W_{t_{0}})=(x_{0},w_{0})$ at time $t=t_{0}$, the system of differential equations (A.2) has a unique solution over the time interval $[t_{0},t_{1}]$. Furthermore, the solution does not blow up in any finite time, since from Theorem 2.1 we know that the energy functional $\mathcal{E}_{t}$ (2.8) is non-increasing, so in particular, the Bregman divergence $D_{h}(x^{*},X_{t}+\frac{t}{p}\dot{X}_{t})$ is bounded above by a constant. Since $t_{1}$ is arbitrary, this shows that (A.2) has a unique maximal solution, i.e., $t_{1}$ can be extended to $t_{1}\to+\infty$.

In the above argument we have started at time $t_{0}>0$, because the vector field in (A.2) has a singularity at $t=0$. For $p=2$, Su et al. and Krichene et al. treat the case when we start at $t=0$ with initial condition $(X_{0},W_{0})=(x_{0},\nabla h(x_{0}))$, so that $\dot{X}_{0}=0$. In that case, they show that the system (A.2) still has a unique solution for all time $[0,\infty)$, by replacing the $p/t$ coefficient by the approximation $p/\max\{t,\delta\}$ for $\delta>0$ and letting $\delta\to 0$. We can adapt this technique to the more general case (A.2); alternatively, we can appeal to the time dilation property and state that since the general system (A.2) is the result of speeding up the $p=2$ case by time dilation function $\tau(t)=t^{p/2}$, once we know a unique solution exists for $p=2$, we can also conclude that it exists for all $p>0$.

A.3 Proof of Theorem 3.1

We define the following function, which is a generalization of Nesterov’s estimate function from ,

The estimate function $\psi_{k}$ arises as the objective function that the sequence $z_{k}$ is optimizing in (3.6b). Indeed, the optimality condition for the $z_{k}$ update (3.6b) is

and since $x_{0}=z_{0}$, we can write this equation as $\nabla\psi_{k}(z_{k})=0$. Since $\psi_{k}$ is a convex function, this means $z_{k}$ is the minimizer of $\psi_{k}$. Thus, we can equivalently write the update for $z_{k}$ as

For proving the convergence rate for the algorithm (3.6), we have the following property.

We proceed via induction on $k\geq 0$. The base case $k=0$ is true since both sides equal zero. Now assume (A.6) holds for some $k\geq 0$; we will show it also holds for $k+1$.

Since $h$ is $1$-uniformly convex of order $p$, the rescaled Bregman divergence $\frac{1}{\epsilon}D_{h}(x,x_{0})$ is $(\frac{1}{\epsilon})$-uniformly convex. Thus, the estimate function $\psi_{k}$ (A.4) is also $(\frac{1}{\epsilon})$-uniformly convex of order $p$. Since $z_{k}$ is the minimizer of $\psi_{k}$, $\nabla\psi_{k}(z_{k})=0$, so for all $x\in\mathcal{X}$ we have

Applying the inductive hypothesis (A.6) and using the convexity of $f$ gives us

We now add $Cp(k+1)^{(p-1)}[f(y_{k+1})+\langle\nabla f(y_{k+1}),x-y_{k+1}\rangle]$ to both sides of the equation to obtain

where $\tau_{k}=\frac{p(k+1)^{(p-1)}}{(k+1)^{(p)}}=\frac{p}{k+p}$, and where we have also used the definition of $x_{k+1}$ as a convex combination of $y_{k}$ and $z_{k}$ with weight $\tau_{k}$ (3.6a).

Note that the first term in (A.7) gives our desired inequality (A.6) for $k+1$. So to finish the proof, we have to prove the remaining terms in (A.7) are nonnegative. We do so by applying two inequalities. We first apply the inequality (3.7) to the term $\langle\nabla f(y_{k+1}),x_{k+1}-y_{k+1}\rangle$, so from (A.7) we have

Next, we apply the Fenchel-Young inequality [28, Lemma 2]

with the choices $u=\epsilon^{-\frac{1}{p}}(x-z_{k})$ and $s=\epsilon^{\frac{1}{p}}Cp(k+1)^{(p-1)}\nabla f(y_{k+1})$. Then from (A.8), we obtain

Notice that $\{(k+1)^{(p-1)}\}^{\frac{p}{p-1}}\leq(k+1)^{(p)}$. Then from the assumption $C\leq M^{p-1}/p^{p}$, we see that the second term inside the parentheses is nonnegative. Hence we conclude the desired inequality $\psi_{k+1}(x)\geq C(k+1)^{(p)}f(y_{k+1})$. Since $x\in\mathcal{X}$ is arbitrary, it also holds for the minimizer $x=z_{k+1}$ of $\psi_{k+1}$, finishing the induction. ∎

With Lemma A.1 in hand, we can complete the proof of Theorem 3.1.

Since $f$ is convex, we can bound the estimate sequence $\psi_{k}$ by

This holds for all $x\in\mathcal{X}$, and in particular for the minimizer $x^{*}$ of $f$. Combining the bound with the result of Lemma A.1, and recalling that $z_{k}$ is the minimizer of $\psi_{k}$, we get

Rearranging and dividing by $Ck^{(p)}$ gives us the desired convergence rate (3.9). ∎

A.4 Proof of Lemma 3.2

We follow the approach of [28, Lemma 6]. Since $y$ solves the optimization problem (3.11), it satisfies the optimality condition

Furthermore, since $\nabla^{p-1}f$ is $\frac{(p-1)!}{\epsilon}$-Lipschitz, we have the following error bound on the $(p-2)$-nd order Taylor expansion of $\nabla f$,

Substituting (A.10) to (A.11) and writing $r=\|y-x\|$, we obtain

Squaring both sides, expanding, and rearranging the terms, we get the inequality

Note that if $p=2$, then the first term in (A.13) already implies the desired bound (3.12). Now assume $p\geq 3$. The right-hand side of (A.13) is of the form $A/r^{p-2}+Br^{p}$, which is a convex function of $r>0$ and minimized by $r^{*}=\left\{\frac{(p-2)}{p}\frac{A}{B}\right\}^{\frac{1}{2p-2}}$, yielding a minimum value of

Substituting the values $A=\frac{\epsilon}{2N}\|\nabla f(y)\|_{*}^{2}$ and $B=\frac{1}{2N\epsilon}(N^{2}-1)$ from (A.13), we obtain

To obtain the first inequality of (3.13), we use Cauchy-Schwarz inequality on (3.12),

and cancel out $\|\nabla f(y)\|_{*}$ from both sides. For the second inequality of (3.13), we use triangle inequality on the left hand side of (A.12),

Rearranging the terms and taking the $(p-1)$-st root of both sides gives us the result (3.13). ∎

A.5 Proof of Theorem 3.3

This proof follows the approach in the proof of [28, Theorem 1]. We first prove the following lemma. Here $\delta_{k}=f(x_{k})-f(x^{*})\geq 0$ denotes the residual value at iteration $k$.

Under the setting of Theorem 3.3, we have

Since $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, by the Taylor remainder theorem we have the bound

Then from the definition of $x_{k+1}$ (3.14), we have

Plugging in $x=x_{k}$ on the right-hand side of (A.15) shows that $f(x_{k+1})\leq f(x_{k})$; that is, the algorithm (3.14) is a descent method. In particular, for all $k\geq 0$ we have $\|x_{k}-x^{*}\|\leq R$, where $R=\sup_{x\colon f(x)\leq f(x_{0})}\|x-x^{*}\|$ is the radius of the level set as defined in Theorem 3.3. Moreover, plugging in $x=x^{*}$ on the right-hand side of (A.15) gives us

Now for any $\lambda\in$, consider the midpoint

By Jensen’s inequality, $f(x_{\lambda})\leq\lambda f(x^{*})+(1-\lambda)f(x_{k})$. We also have $\|x_{\lambda}-x_{k}\|=\lambda\|x_{k}-x^{*}\|\leq\lambda R$. Plugging in the point $x_{\lambda}$ to the right-hand side of (A.15) gives

With the notation $\delta_{k}=f(x_{k})-f(x^{*})$, we can write the last inequality as

The right-hand side is a convex function of $\lambda$, which is minimized at $\lambda^{*}=\left\{\frac{\epsilon}{N+1}\frac{\delta_{k}}{R^{p}}\right\}^{\frac{1}{p-1}}$. Note that $\lambda^{*}\in$ by (A.16). Plugging in $\lambda^{*}$ to (A.17) yields the desired bound (A.14). ∎

With Lemma A.2, we can complete the proof of Theorem 3.3.

Define the energy functional $e_{k}=\delta_{k}^{-\frac{1}{p-1}}$. We can write

Since $\delta_{k+1}\leq\delta_{k}$, we can upper bound the summation in the denominator of (A.18) by $(p-1)\delta_{k}^{\frac{p-2}{p-1}}$. We use Lemma A.2 to lower bound $\delta_{k}-\delta_{k-1}$, obtaining

Summing (A.19) and telescoping the terms, we get

which gives us the desired conclusion (3.15). ∎

A.6 Proof of Theorem 3.4

We write the higher-order gradient algorithm (3.14) (with $N=1$) as

Our goal is to express the sequence $x_{k}$ as a discretization $x_{k}=X_{t}$, $x_{k+1}=X_{t+\delta}\approx X_{t}+\delta\dot{X}_{t}$ of some continuous-time curve $X_{t}$ with time step $\delta>0$, which will be a function of $\epsilon$. To that end, we write $u=\delta v$, so (A.20) becomes

Eliminating the constant term $f(x_{k})$ from the right-hand side, which does not change the minimizer, and canceling a factor of $\delta$, we get

We see that the first term in the objective function does not depend on $\delta$. As $\epsilon\to 0$, for the equation to have a meaningful limit, we have to set $\delta^{p-1}=\epsilon$, so the last term in the objective function becomes a constant. On the other hand, the middle terms all have dependence on $\delta=\epsilon^{\frac{1}{p-1}}$, so as $\epsilon\to 0$, those terms vanish. Thus, the limit as $\epsilon\to 0$ is

Equivalently, $\dot{X}_{t}$ satisfies the optimality condition

This gives us the relation $\|\nabla f(X_{t})\|_{*}=\|\dot{X}_{t}\|^{p-1}$, so we can also write (A.22) as

which is the rescaled gradient flow as claimed in (3.16).

We note that the rescaled gradient flow (3.16) is a descent method, since

Now to establish the convergence rate of the rescaled gradient flow (3.16), we consider the energy functional

which is the same energy functional as in the discrete-time convergence proof in Appendix A.5. The energy functional $\mathcal{E}_{t}$ has time derivative

If $X_{t}$ satisfies the rescaled gradient flow equation (3.16), then $\dot{\mathcal{E}}_{t}$ simplifies to

By the convexity of $f$ and the Cauchy-Schwarz inequality, we have

Since the rescaled gradient flow is a descent method, we have $\|X_{t}-x^{*}\|\leq R$. Therefore, from (A.24) we get the bound

This means that $\mathcal{E}_{t}$ increases at least linearly, so

which gives us the desired result (3.17). ∎

From the proof above, we see that rescaled gradient flow (3.16) is a generalization of the usual gradient flow (the case $p=2$) which is obtained by replacing the squared norm by the $p$-th power of the norm in the variational formulation (A.21). It turns out that when the objective function is the $p$-th power of the norm, $f(x)=\frac{1}{p}\|x\|^{p}$, the rescaled gradient flow (3.16) reduces to an explicit equation. Specifically, in this case we have $\nabla f(x)=\|x\|^{p-2}x$, so $\|\nabla f(x)\|_{*}=\|x\|^{p-1}$. Therefore, the rescaled gradient flow equation (3.16) becomes

which is now independent of $p$, and has an explicit solution $X_{t}=e^{-t}X_{0}$.

In the proof above, we can also use the following alternative energy functional,

where (A.26a) follows from the convexity of $f$, and in (A.26b) we have substituted the rescaled gradient flow dynamic (3.16). We now apply the Fenchel-Young inequality (A.9) with $s=t^{p-1}\nabla f(X_{t})$ and $u=-(p-1)(X_{t}-x^{*})$, to obtain

which is exactly the same bound as claimed in (3.17).

Appendix B Further properties

Similar to the polynomial case in Section 3, in this section we study the subfamily of Bregman Lagrangian (2.1) with the following choice of parameters, parameterized by $c>0$,

The parameters (B.1) satisfy the ideal scaling condition (2.2), with an equality on the first condition (2.2a). The Euler-Lagrange equation (2.6) in this case is given by

and by Theorem 2.1, it has an $O(e^{-ct})$ rate of convergence. Thus, whereas the polynomial Lagrangian flow (3.2) has a polynomial rate of convergence, the exponential Lagrangian flow (B.2) has an exponential rate of convergence. Furthermore, from the time-dilation property in Theorem 2.2, we see that we can obtain the exponential curve (B.2) by speeding up the polynomial curve (3.2) using a time-dilation function $\tau(t)=e^{ct/p}$.

However, unlike the polynomial Lagrangian flow (3.2), the process of discretizing the exponential Lagrangian flow (B.2) is not as straightforward. Following the same approach as the polynomial family, we write the second-order equation (B.2) as the following pair of first-order equations:

Now we discretize $X_{t}$ and $Z_{t}$ into sequences $x_{k}$ and $z_{k}$ with time step $\delta>0$, so that $t=\delta k$ as before. In doing so, we can write (B.3) as the following discrete-time equations similar to (3.4) and (3.5):

Note that the weight in (B.4a) is independent of time, but depends on $\delta$, and (B.4b) suggests the step size $\epsilon=\delta$ in the algorithm. If our analogy between continuous and discrete-time convergence holds, then given the $O(e^{-ct})$ convergence rate in continuous time, we expect a matching $O(\frac{1}{\epsilon}e^{-ck})$ convergence rate in discrete time. However, it is not clear how to obtain that rate via (B.4). If we try to adapt the proof of Theorem 3.1, we find that in order to conclude a convergence rate $O(\delta e^{-c\delta k})$, we need to introduce a sequence $y_{k}$ satisfying the following analog of inequality (3.7) (with the ideal choice $p=\infty$):

Notice that the rates are consistent if we set $\epsilon=\delta=1$. However, the condition (B.5) means we need to make a constant improvement in each iteration from $x_{k}$ to $y_{k}$, although we are also free on how we choose to construct $y_{k}$ and impose any assumptions on $f$.

In the remainder of this section, we approach this problem from a discrete-time perspective, and study the performance of the higher-order gradient algorithm (3.14) and its accelerated variant (3.18) when $f$ is uniformly convex.

In this section we show that the higher-order gradient algorithm (3.14) has an exponential convergence rate when the objective function $f$ is uniformly convex of order $p\geq 2$; this generalizes the results in [28, Section 5] for the case $p=3$, and the classical result of gradient descent for the case $p=2$ .

Specifically, we have the following result. Recall the definition of smoothness in (3.10), and the definition of uniform convexity in (3.8).

Suppose $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, and $\sigma$-uniformly convex of order $p$. Then the $p$-th order gradient algorithm (3.14) with $N>1$ has convergence rate

where $L=(N^{2}-1)^{\frac{p-2}{2p-2}}/(2N)$, and $\kappa=\epsilon\sigma$ is the inverse condition number (which we assume is small).

By inequality (3.12) from Lemma 3.2, we know that since $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$,

where $L=(N^{2}-1)^{\frac{p-2}{2p-2}}/(2N)$. Since $f$ is convex, we have $f(x_{k})-f(x_{k+1})\geq\langle\nabla f(x_{k+1}),x_{k}-x_{k+1}\rangle$. Furthermore, since $f$ is $\sigma$-uniformly convex of order $p$, from [28, Lemma 3] we also have

Combining these inequalities and recalling the definition $\kappa=\epsilon\sigma$ gives us

Note that by the smoothness of $f$, as in (A.15), we can write $f(x_{1})\leq\min_{x}\{f(x)+\frac{N+1}{\epsilon p}\|x-x_{0}\|^{p}\}\leq f(x^{*})+\frac{N+1}{\epsilon p}\|x_{0}-x^{*}\|^{p}$. Furthermore, since we assume the inverse condition number $\kappa=\epsilon\sigma$ is small, we can write $1+L\kappa^{\frac{1}{p-1}}\approx\exp(L\kappa^{\frac{1}{p-1}})$. Therefore, (B.8) yields the desired convergence rate (B.6). ∎

Notice that the result of Theorem B.1 matches the desired convergence rate $O(\frac{1}{\epsilon}e^{-ck})$ discussed in Appendix B.1, with $c=L\kappa^{\frac{1}{p-1}}$.

As a side remark, we note that the rescaled gradient flow also has an exponential convergence rate when the objective function $f$ is uniformly convex. However, notice that the following continuous-time convergence rate only depends on the uniform convexity constant of $f$, whereas the discrete-time convergence rate above also depends on the Lipschitz constant for the higher-order smoothness of $f$.

If $f$ is $\sigma$-uniformly convex of order $p$, then the rescaled gradient flow (3.16) has convergence rate

As we saw in (B.7), the uniform convexity of $f$ implies the inequality

Using this inequality and plugging in the rescaled gradient flow equation (3.16), we have

Dividing both sides by $f(X_{t})-f(x^{*})$ and integrating, we get the desired convergence rate (B.9). ∎

B.1.2 Exponential convergence rate of accelerated method with restart scheme

We now show that a variant of the accelerated gradient method (3.18) with a restart scheme also attains an exponential convergence rate, with a better dependence on the condition number $\kappa$ than the higher-order gradient method as in Appendix B.1.1.

Specifically, we consider the following variant of the accelerated gradient method (3.18),

In (B.10), for simplicity we have explicitly set the constant $N$ in (3.18b) to be $N=2$, and set $C$ in (3.18c) to be $C=1/(4p)^{p}$, which satisfies the condition $C\leq(N^{2}-1)^{\frac{p-2}{2}}/((2N)^{p-1}p^{p})$. Furthermore, for the $z$-update (B.10c) we have used the equivalent version (A.5) where we unroll the recursion, and we have also replaced the Bregman divergence in the $z$-update (3.18c) by the rescaled $p$-th power $d_{p}(z)=\frac{2^{p-2}}{p}\|z-x_{0}\|^{p}$, which is $1$-uniformly convex of order $p$. The proof of Theorem 3.1 still holds in this case, so we have the guarantee

Then we define the following restart scheme, which proceeds by running the accelerated method (B.10) for some number of iterations at each step,

Suppose $f$ is $\frac{(p-1)!}{\epsilon}$-smooth of order $p-1$, and $\sigma$-uniformly convex of order $p$. Let $\hat{x}_{k}$ be the output of running the restart scheme (B.12) for $k/m$ times with $m=8p/\kappa^{\frac{1}{p}}$, where $\kappa=\epsilon\sigma$ is the inverse condition number, and let $\hat{y}_{k}=G_{p,\epsilon,2}(\hat{x}_{k})$ be the output of running one step of the gradient update (3.11) with input $\hat{x}_{k}$. Then we have the convergence rate

Since $f$ is $\sigma$-uniformly convex of order $p$, and by the bound (B.11), we have

where the last inequality follows from our choice of $m$. Thus, an execution of (B.12) with $m$ iterations of the accelerated method reduces the distance to optimum by a factor of at least $1/e$. Iterating (B.14), we obtain $\|\hat{x}_{k}-x^{*}\|^{p}\leq e^{-k/m}\|\hat{x}_{0}-x^{*}\|^{p}$. To convert this into a bound on the function value, we use the smoothness of $f$. As noted in (A.15), since $\hat{y}_{k}$ is the output of one step of the gradient update (3.11) with input $\hat{x}_{k}$, we have $f(\hat{y}_{k})-f(x^{*})\leq\frac{3}{\epsilon p}\|\hat{x}_{k}-x^{*}\|^{p}$. This gives the desired bound (B.13). ∎

The result of Theorem B.3 matches the desired convergence rate $O(\frac{1}{\epsilon}e^{-ck})$ as discussed in Appendix B.1 with $c=\frac{1}{8p}\kappa^{\frac{1}{p}}$. Note that this convergence rate has a better dependence on the inverse condition number $\kappa=\epsilon\sigma$ than the higher-order gradient algorithm as in Theorem B.1, because $\kappa^{\frac{1}{p}}>\kappa^{\frac{1}{p-1}}$ for small $\kappa$. This generalizes the conclusion of [28, Section 5] for the case $p=3$. However, as noted previously, the link to continuous time is not as clear as that of the polynomial family.

B.2 Hessian vs. Bregman Lagrangian

In a Hessian manifold, the metric is generated by the Hessian $\nabla^{2}h$ of the distance-generating function $h$. So for example, the gradient flow equation in the Euclidean case, $\dot{X}_{t}=-\nabla f(X_{t})$, which can be written as

in general becomes the natural gradient flow $\dot{X}_{t}=-[\nabla^{2}h(X_{t})]^{-1}\nabla f(X_{t})$, or equivalently,

which is obtained by replacing the Euclidean squared norm $\|v\|^{2}=\langle v,v\rangle$ by the Hessian metric

At the Lagrangian level, recall that a starting point of our work is the differential equation $\ddot{X}_{t}+\frac{3}{t}\dot{X}_{t}+\nabla f(X_{t})=0$ for accelerated gradient descent , which we observe is the Euler-Lagrange equation for the damped Lagrangian

How should we generalize this Lagrangian to the non-Euclidean case? From our discussion on natural gradient flow, a natural guess is to replace the Euclidean metric in (B.15) by the Hessian metric. Thus, we are led to consider the following family of Hessian Lagrangians:

involves the third-order derivative $\nabla^{3}h$ (which comes from being the derivative of the metric tensor $\nabla^{2}h$). This makes the analysis difficult, preventing us from obtaining a convergence rate for (B.17). Furthermore, the presence of $\nabla^{3}h$ in the equation makes it difficult to implement as an efficient discrete-time algorithm.

On the other hand, our work shows that the “correct” way to generalize (B.15) to the non-Euclidean case is to use the Bregman divergence, rather than Hessian metric. This results in the general Bregman Lagrangian family (2.1), which requires an additional parameter $\alpha_{t}$ controlling the amount of interaction between the position $X$ and velocity $V$. When the parameters are coupled in an ideal scaling, the Bregman Lagrangian produces dynamics that converge at a provable rate. This is achieved via the design of a corresponding Lyapunov function (the energy functional $\mathcal{E}_{t}$ (2.8)), whose form is intimately tied to the use of the Bregman divergence in the Lagrangian. Furthermore, for the polynomial family, we can discretize the resulting dynamics as a discrete-time algorithm (3.18) that does not require the Hessian $\nabla^{2}h$, but only the gradient $\nabla h$.

It is interesting to consider whether the Hessian Lagrangian (B.16) has useful properties, and how it relates to the Bregman Lagrangian. For a small displacement $\varepsilon>0$ we know that Bregman divergence approximates the Hessian metric, i.e., $D(x+\varepsilon v,x)\approx\frac{\varepsilon^{2}}{2}\|v\|_{\nabla^{2}h(x)}^{2}$. Setting $\varepsilon=e^{-\alpha_{t}}$, this suggests that the Bregman Lagrangian (2.1) is approximating the Hessian Lagrangian $\mathcal{L}_{\text{Hess}}(X,V,t)=e^{\gamma_{t}-\alpha_{t}}\left(\frac{1}{2}\|V\|_{\nabla^{2}h(X)}^{2}-e^{2\alpha_{t}+\beta_{t}}f(X)\right)$. However, this argument assumes $\epsilon$ is small, whereas in our particular case of interest (the polynomial subfamily in Section 3) the value of $\epsilon=e^{-\alpha_{t}}=\frac{t}{p}$ is growing over time.

B.3 Gradient vs. Lagrangian flows

In the Euclidean case, we can think of gradient flow as describing the behavior of a damped Lagrangian system “in an asymptotic regime in which dissipative effects play such an important role, that the effects of forcing and dissipation compensate each other” [34, p. 646]. That is, the gradient flow equation $\dot{X}_{t}=-\nabla f(X_{t})$ can be seen as the strong-friction limit $\lambda\to\infty$ of the equation $\ddot{X}_{t}+\lambda\dot{X}_{t}+\lambda\nabla f(X_{t})=0$.

This is perhaps more apparent if we define $m=1/\lambda$ to be the “mass” of the fictitious particle, so the equation of motion becomes

which is the Euler-Lagrange equation of the damped Lagrangian

where the damping factor $e^{t/m}$ also scales with $m$. In the massless limit $m\to 0$, we indeed recover gradient flow from (B.18). In the following, we show that this result also holds more generally, both for natural gradient flow (as the massless limit of a Bregman Lagrangian flow) and for the rescaled gradient flow (as the massless limit of a Lagrangian flow which uses the $p$-th power of the norm).

However, notice that in all these cases, the momentum variable $P=\frac{\partial\mathcal{L}}{\partial V}$ becomes infinite as $m\to 0$. For instance, $P=me^{t/m}V$ for (B.19), and $me^{t/m}\to\infty$. This means as $m\to 0$, the particle also becomes more massive and has more inertia. Thus, gradient flow is the limiting case where the infinitely massive particle simply rolls downhill and stops at the minimum $x^{*}$ as soon as the force $-\nabla f$ vanishes, without oscillation (which is damped by the infinitely strong friction). In this view, moving from a first-order gradient algorithm to a second-order Lagrangian (accelerated) algorithm does not amount to preventing oscillation; rather, it is the opposite, by unwinding the curve to finite momentum where it can travel faster, albeit with some oscillation.

which is the Bregman Lagrangian (2.1) with parameters $\alpha_{t}=-\log m$, $\beta_{t}=\log m$, and $\gamma_{t}=t/m$ (which satisfy the ideal scaling (2.2)). Note that (B.20) recovers (B.19) in the Euclidean case. The Euler-Lagrange equation (2.6) for the Lagrangian (B.20) is given by

Multiplying the equation by $m$ and letting $m\to 0$, we recover

which is the natural gradient flow equation. In this case the momentum variable is $P=\frac{\partial\mathcal{L}}{\partial V}=e^{t/m}(\nabla h(X+mV)-\nabla h(X))\approx me^{t/m}\nabla^{2}h(X)V$, so we still have $P\to\infty$ as $m\to 0$.

where we use the $p$-th power of the norm to measure the kinetic energy. Note that (B.21) recovers (B.19) in the case $p=2$. The Euler-Lagrange equation is

which is equivalent to the rescaled gradient flow (3.16). In this case the momentum variable is $P=\frac{\partial\mathcal{L}}{\partial V}=me^{t/m}\|V\|^{p-2}V$, which still goes to infinity as $m\to 0$.

B.4 Bregman Hamiltonian

In this section we define and compute the Bregman Hamiltonian corresponding to the Bregman Lagrangian. In general, given a Lagrangian $\mathcal{L}(X,V,t)$, its Hamiltonian is defined by

where $P=\frac{\partial\mathcal{L}}{\partial V}$ is the momentum variable conjugate to position.

For the Bregman Lagrangian (2.1), the momentum variable is given by

We can invert this equation to solve for the velocity $V$,

where $h^{*}$ is the conjugate function to $h$ (recall the definition in (A.3)), and we have used the property that $\nabla h^{*}=[\nabla h]^{-1}$. So for the first term in the definition (B.23) we have

Next, we write the Bregman Lagrangian $\mathcal{L}(X,V,t)$ in terms of $(X,P,t)$. We can directly substitute (B.25) to the definition (2.1) and calculate the result. Alternatively, we can use the property that the Bregman divergences of $h$ and $h^{*}$ satisfy $D_{h}(y,x)=D_{h^{*}}(\nabla h(x),\nabla h(y))$. Therefore, we can write the Bregman Lagrangian (2.1) as

where in the second step we have used the relation $\nabla h(X+e^{-\alpha_{t}}V)=\nabla h(X)+e^{-\gamma_{t}}P$ from (B.24), and in the last step we have expanded the Bregman divergence.

Substituting these calculations into (B.23) and simplifying, we get the Hamiltonian

Since $X=\nabla h^{*}(\nabla h(X))$, we can also write this result in terms of the Bregman divergence of $h^{*}$,

We call the Hamiltonian (B.26) the Bregman Hamiltonian. Notice that whereas the Bregman Lagrangian takes the form of the difference between the kinetic and potential energy, the Bregman Hamiltonian takes the form of the sum of the kinetic and potential energy. (However, note that the kinetic energy is slightly different: it is $D_{h^{*}}(\nabla h(X)+e^{-\gamma_{t}}P,\,\nabla h(X))=D_{h}(X,X+e^{-\alpha_{t}}V)$ in the Hamiltonian (B.26), while it is $D_{h}(X+e^{-\alpha_{t}}V,X)$ in the Lagrangian (2.1).)

The second-order Euler-Lagrange equation of a Lagrangian can be equivalently written as a pair of first-order equations

For the Bregman Hamiltonian (B.26), the equations of motion are given by

Notice that the first equation (B.28a) recovers the definition of momentum (B.24). Furthermore, when $\dot{\gamma}_{t}=e^{\alpha_{t}}$, by substituting (B.28a) to (B.28b) we can write (B.28) as

Since $\nabla h(X_{t})+e^{-\gamma_{t}}P_{t}=\nabla h(X_{t}+e^{-\alpha_{t}}\dot{X}_{t})$ by (B.28a), this indeed recovers the Euler-Lagrange equation (2.7).

A Lyapunov function for the Hamiltonian equations of motion (B.28) is the following, which is simply the energy functional (2.8) written in terms of $(X_{t},P_{t},t)$,

The Hamiltonian formulation of the dynamics has appealing properties that seem worthy of further exploration. For example, Hamiltonian flow preserves volume in phase space (Liouville’s theorem); this property has been used in the context of sampling to develop the technique of Hamiltonian Markov chain Monte-Carlo, and may also be useful to help us design better algorithms for optimization. Furthermore, the Hamilton-Jacobi-Bellman equation (which is a reformulation of the Hamiltonian dynamics) is a central object of study in the field of optimal control theory, and it would be interesting to study the Bregman Hamiltonian framework from that perspective.

B.5 Gauge invariance

The Euler-Lagrange equation of a Lagrangian is gauge-invariant, which means it does not change when we transform the Lagrangian by adding a total time derivative,

for any smooth function $G$. We can show this by directly checking that the Euler-Lagrange equation of $\mathcal{L}^{\prime}$ is the same as that of $\mathcal{L}$. Alternatively, this follows from the formulation of the principle of least action, where we fix two points $(x_{0},t_{0})$ and $(x_{1},t_{1})$, and ask for a curve $X$ joining the two endpoints ($X_{t_{0}}=x_{0}$ and $X_{t_{1}}=x_{1}$) that minimizes the action $J(X)=\int_{t_{0}}^{t_{1}}\mathcal{L}(X_{t},\dot{X}_{t},t)dt$. Thus, when the Lagrangian transforms as (B.29), the action only changes to $J^{\prime}(X)=J(X)+\int_{t_{0}}^{t_{1}}\frac{d}{dt}G(X_{t},t)dt=J(X)+G(x_{1},t_{1})-G(x_{0},t_{0})$. Since $(x_{0},t_{0})$ and $(x_{1},t_{1})$ are fixed, this means the new action only differs from the old action by a constant; this implies that the optimal least action curve—namely, the Euler-Lagrange equation—does not change.

In our case, under the ideal scaling condition $\dot{\gamma}_{t}=e^{\alpha_{t}}$ (2.2b), this property implies that the Bregman Lagrangian (2.1) is equivalent to the following Lagrangian

where we have replaced the Bregman divergence $D_{h}(X+e^{-\alpha_{t}}V,X)$ by its first term $h(X+e^{-\alpha_{t}}V)$. Indeed, we can check that the difference between the Bregman Lagrangian (2.1) and the reduced form (B.30) is a total time derivative,

where the last step follows from the ideal scaling $e^{\alpha_{t}}=\dot{\gamma}_{t}$.

The reduced Lagrangian (B.30) is slightly simpler than the Bregman Lagrangian (2.1), and in a sense it makes the roles of $h$ and $f$ more symmetric. It also suggests that the role of $h$ is not so much as measuring the distance via the Hessian metric or Bregman divergence, but rather, as evaluating the extrapolated future point $X_{t}+e^{-\alpha_{t}}\dot{X}_{t}$.

B.6 Natural motion

A natural motion is the motion of a particle when it experiences no force. In the physical world, the natural motion of a particle is a straight-line motion with constant velocity. But for the Bregman Lagrangian, which describes a dissipative system, the natural motion always converges.

Specifically, the Bregman Lagrangian (2.1) in the case of zero (or constant) potential function $f\equiv 0$ is $\mathcal{L}(X,V,t)=e^{\alpha_{t}+\gamma_{t}}D_{h}(X+e^{-\alpha_{t}}V,X)$. Assuming the ideal scaling $\dot{\gamma}_{t}=e^{\alpha_{t}}$ (2.2b), its Euler-Lagrange equation is given by (2.7), which in this case is