# Approximate backwards differentiation of gradient flow

Yushen Huang Taehoon Lee Yifan Sun

(Stony Brook University)

###### Abstract

The gradient flow (GF) is an ODE for which its explicit Euler’s discretization
is the gradient descent method. In this work, we investigate a family of
methods derived from _approximate implicit discretizations_ of (GF), drawing
the connection between larger stability regions and less sensitive
hyperparameter tuning. We focus on the implicit $\tau$-step backwards
differentiation formulas (BDFs), approximated in an inner loop with a few
iterations of vanilla gradient descent, and give their convergence rate when
the objective function is convex, strongly convex, or nonconvex. Numerical
experiments show the wide range of effects of these different methods on
extremely poorly conditioned problems, especially those brought about in
training deep neural networks.

##  1 INTRODUCTION

Recently, there has been an increased interest in the dynamical systems
interpretation of optimization methods (Su et al.,, 2014; Shi et al.,, 2019;
Zhang et al.,, 2018; Romero and Benosman,, 2020); e.g., when an iterative
method is viewed as a discretization of an ODE. In particular, it has been
shown that having a method better mimic the ODE can have favorable
optimization properties, e.g. acceleration under extreme landscape conditions
(Zhang et al.,, 2018) and better quality step directions (Chen et al.,, 2021).
The most famous example of such a dynamical system is gradient flow, defined
by the ODE

$\dot{x}(t)=-\nabla f(x(t)),\quad x(0)={\mathbf{x}}^{0}$

whose explicit Euler discretization is the gradient descent method

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

In fact, it can be shown that very poorly conditioned strongly convex
functions $f$ can be better minimized by explicit multistep discretization
schemes (e.g. those determined by Runge-Kutta steps) via a widening stability
region for step size choices. This simple idea opens a floodgate of new
possibilities, in which the techniques and theories developed by the dynamical
systems community can be used to improve existing optimization methods.

In this paper, we explore approximate implicit discretizations of (GF), and
their effect on badly conditioned problems. We specifically focus on
_backwards differentiation formula (BDF) methods_ , which are used in modeling
and control applications. In particular, we investigate the impact of higher
order implicit methods on convergence stability, and their effect on step size
sensitivity; this is a significant advantage in large-scale machine learning,
where sensitive hyperparameter tuning leads not only to increased
computational cost, but also hints at model overfitting and instability.

###  1.1 Related Work

#### Dynamical systems.

The gradient system is a special kind of dynamical system, a topic which is
well-studied in science and applied mathematics (Strogatz,, 2018; Hirsch et
al.,, 2012; Arnold,, 1992). In general, dynamical systems do not have closed-
form solutions, and in practice we rely a large variety of numerical
discretization methods (Butcher,, 2016), each with different stability and
convergence properties.

#### Gradient flow and higher order discretizations.

In this paper we focus on the dynamical system (GF), of which its Euler
discretization is gradient descent (GD). For example, Su et al., (2014)
analyzed Nesterov accelerated gradient descent method as the discretization of
a second order ordinary differential equation, and (Shi et al.,, 2019)
analyzed the method formed by using a higher order symplectic discretization
of a related ODE (Shi et al.,, 2021), to provide acceleration. Others have
attempted the same using Runge-Kutta explicit discretizations (Zhang et al.,,
2018). A special example of discretization-inspired optimization improvements
is the extragradient method Korpelevich, (1976), which is widely used in min-
max optimization problems (Du and Pardalos,, 1995) and variational inequality
problems (Facchinei and Pang,, 2003).

#### Implicit methods and the proximal point method.

Implicit methods are used in many applications, such as the discontinuous
Garlekin method (Su et al.,, 2020; Rasetarinera and Hussaini,, 2001; Qiu et
al.,, 2005; Trahan and Dawson,, 2012; Wang et al.,, 2015), as well as in
several stochastic differential equation methods (Kloeden and Platen,, 1992;
Tian and Burrage,, 2001; Buckwar and Winkler, 2007b, ; Buckwar and Winkler,
2007a, ; Li et al.,, 2019). Unlike explicit higher-order discretizations,
implicit methods usually involve solving a system of equations at each step,
which may be much more expensive than their explicit counterparts. This paper
will focus on a simple class of implicit methods known as the _backwards
differentiation methods ( BDF)_, of varying order. The first order BDF method,
applied to the gradient flow, yields the familiar _proximal point method_
(PPM) (Moreau,, 1965). PPM is widely used in optimization, and has close
relationship to the augmented Lagarangian methods (Iusem,, 1999) in the dual
update. Other variants of the PPM include (Nesterov,, 2020), and stochastic
PPM (Asi and Duchi,, 2019).

#### Reduction methods.

Related to this work is the subject of _reduction methods_ (Allen-Zhu and
Hazan,, 2016), which take a function which is not strongly convex, apply a
small term that forces strong convexity, and solves a sequence of such
problems for an improved overall convergence rate. The added proximity term in
proximal methods can be interpreted in this way.

###  1.2 Our contribution

The contribution of our papers are the following:

  * •

We investigate a family of multistep implicit optimization algorithms, and
give convergence rates when $f$ is strongly convex, convex, and nonconvex.

  * •

We give further refined convergence rates for _approximate versions_ of these
implicit methods, and in particular refine limits on allowable step sizes, to
illustrate the stability advantage of such methods.

  * •

We show numerically that approximate implicit discretization methods offer
stability and lower hyperparameter tuning in training deep neural networks, of
which are notoriously ill-conditioned problems.

##  2 METHODS

For a first-order ODE of form

$\dot{x}(t)=G(x,t),$

the _backwards differentiation formula_ of order 1 computes discretizations as

${\mathbf{x}}^{k+1}={\mathbf{x}}^{k}+\Delta G({\mathbf{x}}^{k+1},\Delta k)$

where $\Delta>0$ is a discretization step length ($t=\Delta k$). Note that the
right-hand side requires knowledge of ${\mathbf{x}}^{k+1}$ which is
unavailable at iteration $k$, thus requiring an iterative approach to realize
this relation.

###  2.1 Approximate proximal point method (PPM)

Now consider the unconstrained problem

$\underset{{\mathbf{x}}\in\mathbb{R}^{n}}{\mathbf{minimize}}\quad
f({\mathbf{x}})$ (1)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is an everywhere differentiable
function. Applying (BDF-1) to (GF) gives the _proximal point method_ Moreau,
(1965) (PPM):

${\mathbf{x}}^{k+1}={\mathbf{x}}^{k}-\alpha\nabla f({\mathbf{x}}^{k+1})$

which is the optimality condition of the following proximal operation:

$\displaystyle{\mathbf{x}}^{k+1}$
$\displaystyle=\underset{{\mathbf{x}}}{\mathbf{argmin}}\;f({\mathbf{x}})+\frac{1}{2\alpha}\|{\mathbf{x}}-{\mathbf{x}}^{k}\|^{2}$
(ProxOp) $\displaystyle=:\mathbf{prox}_{\alpha f}({\mathbf{x}}^{k}).$

(In all cases, we use $\|\cdot\|=\|\cdot\|_{2}$ the 2-norm.) This operator is
more commonly known as the _proximal operator_ (Moreau,, 1965). Note that in
general, solving (ProxOp) completely may require as much computational
complexity as minimizing $f$ directly, so in general, we only ever consider
approximate versions of this method. For example, applying $m$ steps of
gradient descent on (ProxOp) using a step size $\beta$ gives an approximate
solution ${\tilde{\mathbf{x}}}^{k+1}\approx{\mathbf{x}}^{k+1}$.

We now consider the convergence behavior of (PPM) under different types of
$f$.

######  Assumption 2.1 (Convexity).

A function $f$ is $\mu$-convex if for all ${\mathbf{x}}$ and ${\mathbf{y}}$,

$f({\mathbf{y}})\geq f({\mathbf{x}})+\langle\nabla
f({\mathbf{x}}),{\mathbf{y}}-{\mathbf{x}}\rangle+\frac{\mu}{2}\|{\mathbf{y}}-{\mathbf{x}}\|^{2},$

Here, it may be that $\mu<0$ (weakly convex), $\mu=0$ (convex), or $\mu>0$
(strongly convex).

######  Theorem 2.1 (PPM, Moreau, (1965)).

Assume $f({\mathbf{x}})$ is $\mu$-convex and continuously differentiable
everywhere. The (PPM) converges at the following rates:

  * •

When $\mu>0$ ($f$ is strongly convex)

$\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|_{2}\leq\left(\frac{1}{1+\alpha\mu}\right)^{k}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|_{2}.$

  * •

When $f$ is $\mu$-weakly convex ($\mu\leq 0$), then

$\min_{0\leq s\leq k}\|\nabla
f({\mathbf{x}}^{k+1})\|^{2}\leq\frac{(f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*}))}{\alpha
k}$ (1)

provided

$0\leq\alpha<\frac{2}{|\mu|}.$

  * •

When $\mu=0$, by convexity, (1) holds, but there is no upper-bound restriction
on step size $\alpha>0$.

###### Proof.

The strongly convex case is a well-known result, and the other two are special
cases of Theorem 2.3. ∎

At first glance, it looks like the big-O rate for an implicit method is not
fundamentally improved over an explicit method like (GD). In fact, the main
benefit of implicit methods for (GF) is the _insensitivity to parameter
choices_ , rather than speedup. Note that when $f$ is convex, then
$L$-smoothness is not actually influencing the performance of (PPM). In the
case that $f$ is nonconvex, then step size choice is based on the _level of
nonconvexity_ , which may be less than the $L$-smoothness parameter. However,
as previously alluded to before, this is not a realistic convergence analysis,
since (ProxOp) is not usually cheaply solvable, and usually is harder if
$\alpha$ is chosen too large. Thus, we must limit realistic convergence rates
to tractable approximate versions of (PPM).

######  Assumption 2.2 (Smoothness).

A function $f$ is $L$-smooth if for all ${\mathbf{x}}$ and ${\mathbf{y}}$,

$\|\nabla f({\mathbf{y}})-\nabla f({\mathbf{x}})\|\leq
L\|{\mathbf{y}}-{\mathbf{x}}\|.$

######  Assumption 2.3 ($\gamma$-contractive).

We say a method for solving (ProxOp) is _$\gamma$ -contractive_ if, given
${\mathbf{x}}$, it produces $\tilde{\mathbf{x}}$ where

$\|\mathbf{prox}_{\alpha
f}({\mathbf{x}})-{\tilde{\mathbf{x}}}\|\leq\gamma\|\mathbf{prox}_{\alpha
f}({\mathbf{x}})-{\mathbf{x}}\|.$

As an example, if $f$ is $\mu$-strongly convex and $L$-smooth, then $m$
iterations of (GD) satisfies ($\gamma$C) with
$\gamma=(1-\frac{\alpha\mu+1}{\alpha L+1})^{m}$. Using an accelerated method
(Nesterov,, 1983) improves $\gamma=(1-\frac{\alpha\mu+1}{\alpha L+1})^{m/2}$.

######  Theorem 2.2 (Approx PPM).

Assume $f$ is $L$-smooth and $\mu$-convex. Suppose that a method for solving
(ProxOp) is $\gamma$-contractive. Then the following convergence results hold.

  * •

When $\mu>0$ ($f$ is strongly convex), then taking

$0\leq\alpha\leq\frac{2\gamma}{\mu(1-\gamma)}$

gives

$\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|\leq\left(\gamma+\frac{1+\gamma}{1+\alpha\mu}\right)^{k}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|_{2}.$

  * •

For general $\mu\leq 0$, the same convergence rate (1) holds, whenever

$\alpha<\min\left({\frac{2-16\gamma^{2}}{|\mu|},\frac{1}{|\mu|}}\right).$

  * •

For convex functions $f$, $\mu=0$ and there is no upper bound on the step size
$\alpha$.

###### Proof.

These are special instances of Theorem 2.4. ∎

From here we can see a clear tradeoff between $\gamma$, the requirement of
contraction (smaller is better) and $\alpha$, the step size for faster
convergence (bigger is better). The flexibility offered by this tradeoff is an
advantage of these implicit methods; in a sense, in cases where it may be
tolerable to allow a computation to run for longer at the benefit of less
chance of catastrophic divergence, this new family of methods may prove
advantageous.

###  2.2 Multistep explicit methods

Before moving on to multistep implicit methods, we review briefly the
framework and benefits of multistep _explicit_ methods, to frame our
expectations. A commonly used family of explicit higher order discretization
methods are the $q$-order Runge-Kutta methods (Butcher,, 2016), in which at
each iteration, we update

$\displaystyle\xi_{i}$ $\displaystyle=$
$\displaystyle\displaystyle\dot{x}\Big{(}(k+\omega_{i})\Delta,\;{\mathbf{x}}^{k}+\sum_{j=1}^{q}A_{ij}\xi_{j}\Big{)},$
$\displaystyle{\mathbf{x}}^{k+1}$ $\displaystyle=$
$\displaystyle{\mathbf{x}}^{k}+\sum_{i=1}^{q}\beta_{i}\xi_{i}.$

for some choice of constants $A_{ij}$, $\xi_{j}$, and $\beta_{i}$ for
$i,j=1,...,q$. Note that here, each update of ${\mathbf{x}}^{k}$ requires an
average over 4 intermediate “samples” of $\dot{x}(t)$, and are averaged
precisely to decrease truncation error (distance between $x(t)$ and
${\mathbf{x}}^{k}$) and widen the stability region (larger choice of $\Delta$
to maintain convergence). Explicit multistep discretizations of (GF) and its
variants have been previously studied Shi et al., (2019), with demonstrated
improved convergence behavior.

###  2.3 Multistep implicit methods

We now consider multistep _implicit_ discretizations of (GF), which can be
expressed as

$\displaystyle{{\mathbf{x}}}^{k+\tau+1}$
$\displaystyle=\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\bar{\alpha}\bar{\xi}\nabla
f({{\mathbf{x}}}^{k+\tau+1})$
$\displaystyle=\mathbf{prox}_{\bar{\alpha}\bar{\xi}f}\left(\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}\right).$
($\tau$-BDM)

for some choice of integer $\tau$ (order of method) and constants $\xi_{i}$
where $\sum_{i}\xi_{i}=1$. As reference, specific choices for the famous
backwards differentiation formulas are given in Table 1. For ease of notation,
we will use $\alpha=\bar{\alpha}\bar{\xi}$ from this point on. Note that
unlike Runge-Kutta explicit methods, the BDMs merge past iterates in a
“running fashion”, and there is no constant $\tau$-factor overhead in
computation. And, when $\tau=1$, the above updates correspond to the classical
(PPM).

| $\bar{\xi}$ | $\xi_{1}$ | $\xi_{2}$ | $\xi_{3}$ | $\xi_{4}$ | $\delta$
---|---|---|---|---|---|---
BDF1 | 1 | 1 |  |  |  | 0
BDF2 | 2/3 | -1/3 | 4/3 |  |  | 1/9
BDF3 | 6/11 | 2/11 | -9/11 | 18/11 |  | 192/121
BDF4 | 12/25 | -3/25 | 16/25 | -36/25 | 48/25 | 6603/625
Table 1: Summary of constants for BDF methods. A longer table is in Appendix A

For higher order BDMs, we also define a method-specific constant

$\delta=(\tau-1)\sum_{j=1}^{\tau-1}\left(\sum_{i=1}^{j}(\tau-i)\xi_{i}^{2}\right)$

whose values we also include in Table 1.

######  Theorem 2.3 (BDM).

Assume $f$ is $\mu$-convex and differentiable everywhere. The ($\tau$-BDM)
converges at the following rates:

  * •

When $\mu>0$ ($f$ is strongly convex),

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq
C_{0}\left(\frac{1}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}$

where
$C_{0}=\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|.$

  * •

For general $\mu\leq 0$, for BDMs where $\delta<\frac{1}{\tau-1}$,

$\min_{0\leq s\leq k}\|\nabla
f({\mathbf{x}}^{k+\tau+1})\|^{2}\leq\frac{C_{1}}{\alpha
k}+\frac{C_{2}\delta}{\alpha^{2}k}$ (1)

provided

$0\leq\alpha<\frac{2-2(\tau-1)\delta}{|\mu|}.$

Here,

$C_{1}=f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*}),\;C_{2}=2\sum_{s=0}^{\tau-1}\|{\mathbf{x}}^{s+1}-{\mathbf{x}}^{s}\|^{2}.$

  * •

When $f$ is convex ($\mu=0$), no upper bound restriction on $\alpha$ is
needed.

The proof is in Appendix B.

Now we give the convergence rates of their approximate forms.

######  Theorem 2.4 (Approx BDM).

Consider $f$ everywhere differentiable and $\mu$-convex. Suppose that a method
solving (ProxOp) is $\gamma$-contractive. Then the ($\tau$-BDM) converges at
the following rates.

  * •

When $\mu>0$ ($f$ is strongly convex) then for

$\alpha<\frac{2\gamma}{\mu(1-\gamma)},$

then

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq
C_{0}\left(\gamma+\frac{1+\gamma}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}$

where
$C_{0}=\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|$.

  * •

For general $\mu\leq 0$, when $\delta<\frac{1}{2\tau-2}$, the same convergence
rate (1) holds, whenever

$\alpha<\min\left({\frac{2-2(\tau-1)\delta-16\gamma^{2}}{|\mu|},\frac{1}{|\mu|}}\right).$

  * •

When $\mu=0$, there is no upper bound on $\alpha$.

The proof is in Appendix B.

######  Remark 2.4.

While there may be no step size restriction on $\alpha$, there is still a
dependency on the method used to solve (ProxOp). For example, if (GD) is used,
there will still be reliance between $\gamma$, $\beta$, and $m$. (See Fig. 1.)

##  3 QUADRATIC ANALYSIS

Note that in all the above cases, the convergence result was significantly
better (linear) when $\mu>0$, e.g. the problem is strongly convex. If we
further simplify analysis to when $f(x)=\tfrac{1}{2}x^{T}Qx$ a simple
quadratic function, where $Q$ is symmetric positive semidefinite with
eigenvalues $\mu\leq\lambda_{i}\leq L$, then we may offer much more precise
rates, and give a refined analysis of the step size sensitivity. Specifically,
each method has linear convergence, so we define a radius of convergence
$0\leq\rho\leq 1$ where

$\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|_{2}\leq\rho^{k}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|_{2},$
(2)

where $\rho$ will depend on $\alpha$.

#### Gradient descent.

As a review, the iterates of (GD) can be summarized as

${\mathbf{x}}^{k}=(I-\alpha Q)^{k}{\mathbf{x}}^{0}$

and therefore

$\|{\mathbf{x}}^{k+1}-{\mathbf{x}}^{*}\|_{2}\leq\|I-\alpha
Q\|_{2}\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|_{2}$

and the bound in (2) holds for $\rho_{{\textsc{GD}}}:=\|I-\alpha Q\|_{2}$ the
spectral norm of matrix $I-\alpha Q$. The requirement for stability is that
$\rho_{\textsc{GD}}\leq 1$. Since

$\rho_{\textsc{GD}}\geq\max\\{1-\alpha\mu,\alpha L-1\\},$

stability is satisfied when $\alpha\leq 2/L$. The optimal $\alpha$ is obtained
as the solution to

$\min_{\alpha}\,\max_{i}|1-\alpha\lambda_{i}|=\min_{\alpha}\max\\{\alpha
L-1,1-\alpha\mu\\}$

which is achieved when $\alpha L-1=1-\alpha\mu$, e.g.
$\alpha=\frac{2}{L+\mu}$. In this regime,

$\rho_{\textsc{GD}}^{*}=\frac{L-\mu}{L+\mu}.$

This standard analysis can be found in many textbooks, e.g. (Polyak,, 1987).

#### Explicit Runge-Kutta.

Let us now compare this with an _explicit_ multistep method; in particular,
the Runge Kutta 4 method uses the following update

${\mathbf{x}}^{k+1}={\mathbf{x}}^{k}-\frac{\alpha}{6}(\xi_{1}+2\xi_{2}+2\xi_{3}+\xi_{4})$

where

$\displaystyle\xi_{1}=-\nabla f({\mathbf{x}}^{k}),\qquad\xi_{2}=-\nabla
f({\mathbf{x}}^{k}+\frac{\alpha\xi_{1}}{2}),$ $\displaystyle\xi_{3}=-\nabla
f({\mathbf{x}}^{k}+\frac{\alpha\xi_{2}}{2}),\qquad\xi_{4}=-\nabla
f({\mathbf{x}}^{k}+\alpha\xi_{3}).$

The resulting iteration scheme

${\mathbf{x}}^{k+1}=(I-\alpha
Q+\frac{\alpha^{2}}{2}Q^{2}-\frac{\alpha^{3}}{6}Q^{3}+\frac{\alpha^{4}}{24}Q^{4}){\mathbf{x}}^{k}.$

Using similar analysis, we see that

$\rho_{\textsc{RK}}=\|I-\alpha
Q+\frac{\alpha^{2}}{2}Q^{2}-\frac{\alpha^{3}}{6}Q^{3}+\frac{\alpha^{4}}{24}Q^{4}\|$

and we can numerically ascertain that $\rho_{\textsc{RK}}\leq 1$ whenever

$\alpha\lesssim\frac{2.785}{L}.$

Note that this is a slightly larger range than for GD, which explains the
larger range of step sizes allowable while still avoiding divergence.

The exact value of the _best_ $\alpha$ and resulting $\rho_{\textsc{RK}}$
occurs at the point where

$1-\alpha\mu+\frac{\alpha^{2}}{2}\mu^{2}-\frac{\alpha^{3}}{6}\mu^{3}+\frac{\alpha^{4}}{24}\mu^{4}\\\
=1-\alpha
L+\frac{\alpha^{2}}{2}L^{2}-\frac{\alpha^{3}}{6}L^{3}+\frac{\alpha^{4}}{24}L^{4}.$

This value is numerically computed for three different values of $(\mu,L)$ in
Table 3.

#### Implicit proximal point.

Finally, we investigate the implicit Euler method, also known as the proximal
point method, which for this quadratic example results in the iteration scheme

${\mathbf{x}}^{k+1}=(I+\alpha Q)^{-1}{\mathbf{x}}^{k}$

and thus

$\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|_{2}\leq\rho_{{\textsc{PPM}}}^{k}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|_{2}$

with

$\rho_{{\textsc{PPM}}}=\|(I+\alpha Q)^{-1}\|_{2}$

which is less than 1 for _any_ $\alpha>0$. At first glance, this seems like an
unnaturally good deal, but we should realize that the proximal point method is
rarely used in its exact form in practice, and is usually solved using an
inexact version, e.g.

${\mathbf{x}}^{k+1}\approx\underset{{\mathbf{x}}}{\mathbf{argmin}}\;\frac{1}{2}{\mathbf{x}}^{T}Q{\mathbf{x}}+\frac{1}{2\alpha}\|{\mathbf{x}}-{\mathbf{x}}^{k}\|_{2}^{2}$

using perhaps $m$ steps of gradient descent with a step size of $\beta$. This
approximate method may be summarized with the iteration scheme

${\mathbf{x}}^{k+1}=\left(A^{m}+\frac{\beta}{\alpha}\sum_{i=1}^{m}A^{i-1}\right){\mathbf{x}}^{k}$

where

$A=\left(1-\frac{\beta}{\alpha}\right)I-\beta Q.$ (3)

We give the upper limits of $\beta$ such that

$\rho_{{\textsc{APPM}}}:=\|A^{m}+\frac{\beta}{\alpha}\sum_{i=1}^{m}A^{i-1}\|$

is smaller than 1 in Table 2), and the optimal radius of convergence given
best choice of $\beta$ in Table 3.

Figure 1 plots this radius of convergence for the different methods. If we
look at the values of $\alpha$ for which $\rho<1$, we see a widening of values
from (GD) to (RK), and then (APPM). The amount of contraction in the
approximate PPM also (logically) depends on the choice of $\beta$ and $m$,
which decides how close to convergence (ProxOp) is computed.

Figure 1: Radius of convergence $\rho$ for various methods as a function of
$\alpha$. Left: explicit methods. Right: implicit PPM and approximations.

#### Multistep implicit method

We now consider an approximate multistep implicit method, e.g. after $m$
gradient steps with inner step size $\beta$,

$\displaystyle{\mathbf{x}}^{k+\tau+1}$ $\displaystyle=$ $\displaystyle
A^{m}{\mathbf{x}}^{k+\tau}+\frac{\beta}{\alpha}\sum_{j=1}^{m}A^{j-1}\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}$
$\displaystyle\approx$
$\displaystyle\underset{{\mathbf{x}}}{\mathbf{argmin}}\;\frac{1}{2}{\mathbf{x}}^{T}Q{\mathbf{x}}+\frac{1}{2\alpha}\|{\mathbf{x}}-\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}\|_{2}^{2}$

where $A$ is as in (3). Defining also

$B:=\frac{\beta}{\alpha}\sum_{j=1}^{m}A^{j-1}$

allows a first-order view of the system

$\underbrace{\begin{bmatrix}{\mathbf{x}}^{k+\tau+1}\\\
{\mathbf{x}}^{k+\tau}\\\ {\mathbf{x}}^{k+\tau-1}\\\ \vdots\\\
{\mathbf{x}}^{k+2}\\\
\end{bmatrix}}_{{\mathbf{z}}^{k+1}}=\underbrace{\begin{bmatrix}A^{m}+\xi_{\tau}B&\xi_{\tau-1}B&\cdots&\xi_{1}B\\\
I&0&\cdots&0\\\ 0&I&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&0\\\
\end{bmatrix}}_{=:M}\underbrace{\begin{bmatrix}{\mathbf{x}}^{k+\tau}\\\
{\mathbf{x}}^{k+\tau-1}\\\ {\mathbf{x}}^{k+\tau-2}\\\ \vdots\\\
{\mathbf{x}}^{k+1}\\\ \end{bmatrix}}_{{\mathbf{z}}^{k}}$

and defining
${\mathbf{z}}^{*}=\mathbf{vec}({\mathbf{x}}^{*},{\mathbf{x}}^{*},...,{\mathbf{x}}^{*})$
we have our usual contraction

$({\mathbf{z}}^{k}-{\mathbf{z}}^{*})=M^{k}({\mathbf{z}}^{0}-{\mathbf{z}}^{*}).$

Unfortunately, the 2-norm $\|M\|_{2}$ is bounded below by 1, making it
difficult to prove contraction. However, invoking Gelfand’s rule (Gelfand,,
1941), we may show that

$\|{\mathbf{z}}^{k}-{\mathbf{z}}^{*}\|_{2}\lesssim(\lambda_{\max}(M))^{k}\|{\mathbf{z}}^{0}-{\mathbf{z}}^{*}\|_{2}.$
(4)

That is, although $M$ is not symmetric, (4) holds asymptotically.

In Tables 2 and 3, we give the limits on step size $\beta$ and on optimal
$\rho$ by finding values in which $\lambda_{\max}(M)<1$, guaranteeing
stability. Alongside, figure 2 shows $\rho$ over a continuation of $\beta$,
over changing $m$ and $\alpha$. There is a clear correlation with greater
stability and faster convergence with increased $\alpha$ (implicit step size),
as expected. The reliance on $m$ is interestingly unexpected; more $m$ leads
to better solving of the inner prox step, but gives little effect on stability
and can even give negative effect on convergence. This suggests that though
implicit methods are thought as too expensive, very approximate versions are
not only practical, they are very close to optimal.

| $L=2$ | $L=10$ | $L=100$
---|---|---|---
GD | 1 | 0.2 | 0.02
RK-44 | 1.3925 | 0.2785 | 0.02785
PPM(1) | 0.667 | 0.182 | 0.0198
PPM(10) | 0.952 | 0.198 | 0.0200
BDM2 (1) | 0.665 | 0.181 | 0.0200
BDM2 (10) | 0.940 | 0.197 | 0.0200
BDM3 (1) | 0.608 | 0.178 | 0.0200
BDM3 (10) | 0.940 | 0.197 | 0.0200
Table 2: Sensitivity of parameters. Upper bound on best step size to avoid divergence ($\alpha$ for GD and RK-44), $\beta$ for PPM and BDFs). PPM$(\alpha)$ means proximal point, with step size $\alpha$. (Bigger is better.) Here, we pick $m=4$, and notice very little deviation for $m>4$. Here, $\mu=1$. | $L=2$ | $L=10$ | $L=100$
---|---|---|---
GD | 0.333 | 0.818 | 0.980
RK-44 | 0.362 | 0.771 | 0.973
PPM(4,1) | 0.500 | 0.596 | 0.926
PPM(20,1) | 0.500 | 0.500 | 0.724
PPM(4, 10) | 0.0935 | 0.466 | 0.923
PPM(20,10) | 0.0909 | 0.100 | 0.676
BDM2 (4,1) | 0.326 | 0.282 | 0.905
BDM2 (20,1) | 0.303 | 0.211 | 0.457
BDM2 (4,10) | 0.059 | 0.423 | 0.941
BDM2 (20,10) | 0.024 | 0.024 | 0.737
BDM3 (4,1) | 0.377 | 0.451 | 0.923
BDM3 (20,1) | 0.377 | 0.306 | 0.470
BDM3 (4,10) | 0.197 | 0.459 | 0.943
BDM3 (20,10) | 0.197 | 0.165 | 0.739
Table 3: Optimal $\rho$. Minimum $\rho$ given best best step size ($\alpha$
for GD and RK-44), $\beta$ for PPM and BDFs). PPM$(m,\alpha)$ means proximal
point, with step size $\alpha$ over $m$ iterations. (Smaller is better.) Here,
$\mu=1$. Figure 2: Multistep radius of convergence. $\lambda_{\max}(M)$ over
changing values of $\beta$, for $\alpha=1$ and varying $L=\lambda_{\max}(Q)$
and $m$ = number of inner gradient steps.

##  4 NUMERICAL EXPERIMENTS

###  4.1 Random ill-conditioned quadratic problem

We now consider the minimization of the quadratic function
$f({\mathbf{x}})=\tfrac{1}{2}{\mathbf{x}}^{T}Q{\mathbf{x}}$ where
$Q=UDU^{T}\in{\mathbb{R}}^{100\times 100}$ for a randomly generated
orthonormal $U$ and $D={\textbf{diag}}(\\{d_{i}\\}_{i=1}^{m})$ generated in
three ways:

  * •

$d_{i}$ are i.i.d. generated and uniformly distributed across the interval
$[1,2]$, leading to condition number $\kappa\lesssim 2$ (well-conditioned);

  * •

$d_{i}=\frac{i^{2}}{n^{2}}$ (badly conditioned); and

  * •

$d_{i}=\exp(-i)$ (very badly conditioned).

In figure 3, we give the parameter sensitivity, numerically computed, by
plotting the objective loss after $K$ gradient calls, using an initialization
of ${\mathbf{x}}^{0}_{i}\sim{\mathcal{N}}(0,1)$, and where all (ProxOp) are
solved via $m=8$ inner iterations of gradient descent. (This means (GD) runs
for $K$ iterations, (RK44) $K/4$ iterations, and (PPM), (BDM2), and (BDM3) for
$K/8$ iterations.)

Figure 3: Sensitivity of step size $\alpha$ over quadratic example. Top:
$K=24$ gradient calls. Bottom: $K=500$ gradient calls.

While for the most part these numerical experiments show similar observations
to that derived in the previous section, one interesting difference is a clear
benefit in convergence behavior by using $\tau$ larger. Another important
observation is between the poor conditioning example (middle column) and the
very poor conditioning example (right column), which actually has the better
convergence rate. This is because poorly conditioned quadratic problems simply
diminish the importance of lesser eigenvalue directions, and do not really
suffer despite ${\mathbf{x}}^{k}$ convergence much more slowly to
${\mathbf{x}}^{*}$. Finally, note that despite normalizing for gradient calls
(which penalizes implicit methods more), they still offer competitive
performance over (GD) and (RK44), especially in very ill-conditioned problems.

###  4.2 Deep neural network

Figure 5 shows the sensitivity of step sizes used to complete a simple binary
classification task over a deep neural network described by

$\begin{array}[]{ll}\underset{W_{1},...,W_{L}}{\mathbf{min}}&\displaystyle\sum_{k=1}^{p}\log(1+\exp(-s_{L,k}y_{k}))\\\
\mathbf{s.t.}&s_{i,k}=\text{ReLU}(W_{i}s_{i-1,k}),\;i=2,...,L,\\\ &\quad
k=1,...,p\\\ &s_{1,k}=x_{k}.\end{array}$ (5)

The motivation behind this experiment is the observation that deepening a
neural network makes training more sensitive to step sizes, perhaps due to
worsening of conditioning. Therefore, these approximate implicit methods,
which are more robust to bad conditioning, may have better training
properties. Here, we use a small neural network of width 10, but comparatively
deeper depth of $L=7$. (The goal here is not to train a large model, but to
investigate the difficulty caused by depth, which is often mollified by
increased width.)

One thing to note is that direct application of (GD) on (5) will not converge
well, due to the notorious gradient vanishing problem, but in practice this is
easily overcome using gradient normalization. Therefore, the experiments in
this section are applied to the gradient _norm_ flow, e.g.

$\dot{x}(t)=-\frac{\nabla f(x(t))}{\|\nabla f(x(t))\|}.$

In fact, we do see significant widening of the eligible step sizes when an
approximate implicit method is used. This suggests that indeed, in
applications like deep learning where severely bad conditioning is the modus
operandi, approximate implicit methods offer a promising family of easily
implementable improvements.

Figure 4: Binary classification problem for experiment in Fig. 5.  Figure 5:
Accuracy after 100,000 iterations, for 7 layer fully connected neural network
with ReLU nonlinearlity. 5 steps of GD is used in PPM and BDF methods.
Experiments were performed over 10 random trials, with solid line for the
mean, and transparent fills for the min and max values.

##  5 CONCLUSION

The goal of this work is to investigate the use of approximate implicit
discretizations of (GF) in badly conditioned problem settings, as a candidate
for training ill-conditioned deep neural networks. Convergence rates are given
for the pure implicit methods, and their rates under $m$ inner iterations of
gradient descent for approximate implicit methods. We give these rates under
strongly convex, convex, and nonconvex settings, and importantly observe that
their key benefit is _sensitivity to step size choice_ , rather than
convergence speed (though the latter can be improved if larger step sizes are
allowed).

Overall, it is clear that these multistep explicit and approximate implicit
methods may offer significant benefits, both in avoiding fine-tuning step
sizes and in improved convergence behavior, but the exact choice of method and
parameters is not easy to optimize. This is shown theoretically, with general
rates and specific rates over quadratic problems, and numerically, over a
badly conditioned quadratic problem and a deepening neural network.

## References

  * Allen-Zhu and Hazan, (2016) Allen-Zhu, Z. and Hazan, E. (2016).  Optimal black-box reductions between optimization objectives.  Advances in Neural Information Processing Systems, 29.
  * Arnold, (1992) Arnold, V. I. (1992).  Ordinary differential equations.  Springer Science & Business Media.
  * Asi and Duchi, (2019) Asi, H. and Duchi, J. C. (2019).  Stochastic (approximate) proximal point methods: Convergence, optimality, and adaptivity.  SIAM Journal on Optimization, 29(3):2257–2290.
  * (4) Buckwar, E. and Winkler, R. (2007a).  Improved linear multi-step methods for stochastic ordinary differential equations.  Journal of computational and applied mathematics, 205(2):912–922.
  * (5) Buckwar, E. and Winkler, R. (2007b).  Multi-step maruyama methods for stochastic delay differential equations.  Stochastic Analysis and Applications, 25(5):933–959.
  * Butcher, (2016) Butcher, J. C. (2016).  Numerical methods for ordinary differential equations.  John Wiley & Sons.
  * Chen et al., (2021) Chen, Z., Lee, M., and Sun, Y. (2021).  Continuous time frank-wolfe does not zig-zag.  arXiv preprint arXiv:2106.05753.
  * Du and Pardalos, (1995) Du, D.-Z. and Pardalos, P. M. (1995).  Minimax and applications, volume 4.  Springer Science & Business Media.
  * Facchinei and Pang, (2003) Facchinei, F. and Pang, J.-S. (2003).  Finite-dimensional variational inequalities and complementarity problems.  Springer.
  * Gelfand, (1941) Gelfand, I. (1941).  Normierte ringe.  Matematicheskii Sbornik, 9(1):3–24.
  * Hirsch et al., (2012) Hirsch, M. W., Smale, S., and Devaney, R. L. (2012).  Differential equations, dynamical systems, and an introduction to chaos.  Academic press.
  * Iusem, (1999) Iusem, A. N. (1999).  Augmented lagrangian methods and proximal point methods for convex optimization.  Investigación Operativa, 8(11-49):7.
  * Kloeden and Platen, (1992) Kloeden, P. E. and Platen, E. (1992).  Higher-order implicit strong numerical schemes for stochastic differential equations.  Journal of statistical physics, 66(1):283–314.
  * Korpelevich, (1976) Korpelevich, G. M. (1976).  The extragradient method for finding saddle points and other problems.  Matecon, 12:747–756.
  * Li et al., (2019) Li, X., Wu, Y., Mackey, L., and Erdogdu, M. A. (2019).  Stochastic Runge-Kutta accelerates Langevin monte carlo and beyond.  Advances in neural information processing systems, 32.
  * Moreau, (1965) Moreau, J.-J. (1965).  Proximité et dualité dans un espace hilbertien.  Bulletin de la Société mathématique de France, 93:273–299.
  * Nesterov, (2020) Nesterov, Y. (2020).  Inexact accelerated high-order proximal-point methods.  Technical report.
  * Nesterov, (1983) Nesterov, Y. E. (1983).  A method of solving a convex programming problem with convergence rate o(k^2).  In Doklady Akademii Nauk, volume 269, pages 543–547. Russian Academy of Sciences.
  * Polyak, (1987) Polyak, B. T. (1987).  Introduction to optimization.  Optimization Software, Inc, New York.
  * Qiu et al., (2005) Qiu, J., Dumbser, M., and Shu, C.-W. (2005).  The discontinuous galerkin method with lax–wendroff type time discretizations.  Computer methods in applied mechanics and engineering, 194(42-44):4528–4543.
  * Rasetarinera and Hussaini, (2001) Rasetarinera, P. and Hussaini, M. (2001).  An efficient implicit discontinuous spectral galerkin method.  Journal of Computational Physics, 172(2):718–738.
  * Romero and Benosman, (2020) Romero, O. and Benosman, M. (2020).  Finite-time convergence in continuous-time optimization.  In International Conference on Machine Learning, pages 8200–8209. PMLR.
  * Shi et al., (2021) Shi, B., Du, S. S., Jordan, M. I., and Su, W. J. (2021).  Understanding the acceleration phenomenon via high-resolution differential equations.  Mathematical Programming, pages 1–70.
  * Shi et al., (2019) Shi, B., Du, S. S., Su, W., and Jordan, M. I. (2019).  Acceleration via symplectic discretization of high-resolution differential equations.  Advances in Neural Information Processing Systems, 32.
  * Strogatz, (2018) Strogatz, S. H. (2018).  Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.  CRC press.
  * Su et al., (2014) Su, W., Boyd, S., and Candes, E. (2014).  A differential equation for modeling nesterov’s accelerated gradient method: theory and insights.  Advances in neural information processing systems, 27.
  * Su et al., (2020) Su, W., Wang, P., Zhang, Y., and Wu, L. (2020).  Implicit discontinuous galerkin method for the boltzmann equation.  Journal of Scientific Computing, 82(2):1–35.
  * Tian and Burrage, (2001) Tian, T. and Burrage, K. (2001).  Implicit taylor methods for stiff stochastic differential equations.  Applied Numerical Mathematics, 38(1-2):167–185.
  * Trahan and Dawson, (2012) Trahan, C. J. and Dawson, C. (2012).  Local time-stepping in runge–kutta discontinuous galerkin finite element methods applied to the shallow-water equations.  Computer Methods in Applied Mechanics and Engineering, 217:139–152.
  * Wang et al., (2015) Wang, H., Shu, C.-W., and Zhang, Q. (2015).  Stability and error estimates of local discontinuous galerkin methods with implicit-explicit time-marching for advection-diffusion problems.  SIAM Journal on Numerical Analysis, 53(1):206–227.
  * Zhang et al., (2018) Zhang, J., Mokhtari, A., Sra, S., and Jadbabaie, A. (2018).  Direct runge-kutta discretization achieves acceleration.  Advances in neural information processing systems, 31.

##  Appendix A COMPLETE BDF PARAMETERS LIST

| $\bar{\xi}$ | $\xi_{1}$ | $\xi_{2}$ | $\xi_{3}$ | $\xi_{4}$ | $\xi_{5}$ | $\xi_{6}$
---|---|---|---|---|---|---|---
BDF1 | 1 | 1 |  |  |  |  |
BDF2 | 2/3 | 4/3 | -1/3 |  |  |  |
BDF3 | 6/11 | 18/11 | -9/11 | 2/11 |  |  |
BDF4 | 12/25 | 48/25 | -36/25 | 16/25 | -3/25 |  |
BDF5 | 60/137 | 300/137 | -300/137 | 200/137 | -75/137 | 12/137 |
BDF6 | 60/147 | 360/147 | -450/147 | 400/147 | -225/147 | 72/142 | -10/147

Note that for BDF 4,5, and 6, the stability region will actually _shrink_ ,
making them less useful.

##  Appendix B CONVERGENCE PROOFS

######  Assumption B.1.

$f({\mathbf{x}})$ is convex and continuously differentiable:

$\displaystyle f({\mathbf{y}})\geq f({\mathbf{x}})+\langle\nabla
f({\mathbf{x}}),{\mathbf{y}}-{\mathbf{x}}\rangle$

######  Assumption B.2.

$f({\mathbf{x}})$ is $\mu$\- strongly convex:

$\displaystyle f({\mathbf{y}})\geq f({\mathbf{x}})+\langle\nabla
f({\mathbf{x}}),{\mathbf{y}}-{\mathbf{x}}\rangle+\frac{\mu}{2}\|{\mathbf{y}}-{\mathbf{x}}\|^{2}$

This assumption is needed when for linear convergence.

######  Assumption B.3.

$f({\mathbf{x}})$ is a $L$\- smooth function:

$\|\nabla f({\mathbf{y}})-\nabla f({\mathbf{x}})\|\leq
L\|{\mathbf{y}}-{\mathbf{x}}\|$

In this section, we are giving the main result we have for the multistep
algorithm. We first need a technical lemma:

######  Lemma B.1.

Let $\\{a_{n}\\}_{n=0}^{\infty}$ be a $\tau$ step sublinear sequence:

$a_{n}\leq\sum_{i=1}^{\tau}\xi_{i}a_{n-i}$ (6)

where $\xi_{i}\geq 0$ and $n\geq\tau$. Then there exists $n>i_{1}>i_{2}>\cdots
i_{[\frac{n}{\tau}]}$ such that

$a_{n}\leq\left(\sum_{i=1}^{n}\xi_{i}\right)a_{i_{1}}\leq\left(\sum_{i=1}^{n}\xi_{i}\right)^{2}a_{i_{2}}\leq\cdots\leq\left(\sum_{i=1}^{n}\xi_{i}\right)^{[\frac{n}{\tau}]}a_{i_{[\frac{n}{\tau}}]}$
(7)

for all $k\in\\{1,2,\cdots,[\frac{n}{\tau}]\\}$

###### Proof.

We prove this result by mathematical induction. When $n=\tau$, by using
Holder’s inequality for $p=1,q=\infty$, we will have:

$\displaystyle
a_{\tau}\leq\sum_{i=1}^{\tau}\xi_{i}a_{n-i}\leq\left(\sum_{i=1}^{\tau}\xi_{i}\right)\left(\max_{i\in\\{1,2,\cdots,\tau\\}}a_{\tau-i}\right).$

In this first step, we pick
$i_{1}=\underset{i\in\\{1,2,\cdots,\tau\\}}{\mathbf{argmax}}a_{\tau-i}[a_{n-i}?]$.
This gives

$a_{n}\leq\left(\sum_{i=1}^{\tau}\xi_{i}\right)a_{i_{1}},\qquad
i_{1}=\underset{i\in\\{1,2,\cdots,\tau\\}}{\mathbf{argmax}}a_{\tau-i}.$

Now applying (6) to $a_{i_{1}}$, we can further conclude

$a_{i_{1}}\leq\left(\sum_{i=1}^{\tau}\xi_{i}\right)a_{i_{2}},\qquad
i_{2}=i_{1}-j,\quad j=\underset{i\in
1,2,...,\tau}{\mathbf{argmax}}a_{i_{1}-i}.$

This we can apply recursively, to get

$a_{i_{k}}\leq\left(\sum_{i=1}^{\tau}\xi_{i}\right)a_{i_{k+1}},\qquad
i_{k+1}=i_{k}-j,\quad j=\underset{i\in
1,2,...,\tau}{\mathbf{argmax}}a_{i_{k}-i}.$

This inequality can be chained for as long as there exists $i_{k}\geq 0$. Note
that at each step, the smallest value $i_{k}$ can take is $n-\tau k$. So, the
largest value of $k$ such that $i_{k}\geq 0$ is guaranteed is
$k=\lceil\frac{n}{\tau}\rceil$. ∎

######  Remark B.4.

Lemma B.1 holds for some $i_{[\frac{n}{\tau}]}\in\\{0,1,2,\cdots,\tau-1\\}$.

Now we give the theorem which shows the linear convergence of the algorithm in
$\mu$ strongly convex case:

######  Theorem B.2.

Let $f({\mathbf{x}})$ be a $\mu$-strongly convex and continuously
differentiable function. Suppose that we minimize $f({\mathbf{x}})$ using a
$\tau$steps backward differentiation formula on the gradient flow $\tau$-BDM
where we assume $\sum_{i=1}^{n}\xi_{i}=1$ and $\xi_{i}\geq 0$. Then the
backward gradient descent method will converge linearly

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq\left(\frac{1}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|.$

###### Proof.

By ($\tau$-BDM), we will have

$\displaystyle\|\sum_{i=1}^{\tau}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{*})\|^{2}$
$\displaystyle\overset{\eqref{bdf}}{=}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}+\alpha\nabla
f({\mathbf{x}}^{k+\tau+1})\|^{2}$
$\displaystyle=\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2}+\alpha^{2}\|\nabla
f({\mathbf{x}}^{k+\tau+1})\|^{2}+2\alpha\langle\nabla
f({\mathbf{x}}^{k+\tau+1}),{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\rangle$
$\displaystyle\overset{(a)}{\geq}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2}+\alpha^{2}\mu^{2}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2}+2\alpha\langle\nabla
f({\mathbf{x}}^{k+\tau+1}),{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\rangle$
$\displaystyle\overset{(b)}{\geq}(1+\alpha^{2}\mu^{2})\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2}+2\alpha(\mu\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2})$
$\displaystyle=(1+2\alpha\mu+\alpha^{2}\mu^{2})\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|^{2}$

where (a) and (b) are consequences of $\mu$-strong convexity. As a result, we
will have:

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq\frac{1}{1+\alpha\mu}\|\sum_{i=1}^{\tau}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{*})\|$
(8)
$\displaystyle\leq\frac{1}{1+\alpha\mu}\sum_{i=1}^{\tau}\xi_{i}\|{\mathbf{x}}^{k+i}-{\mathbf{x}}^{*}\|$
(9)

by convexity of the 2-norm. Now by Lemma B.1, we will have
$i_{1},i_{2},\cdots,i_{[\frac{k+\tau+1}{\tau}]}$ such that

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq\left(\frac{1}{1+\alpha\mu}\right)\|{\mathbf{x}}^{i_{1}}-{\mathbf{x}}^{*}\|\leq\left(\frac{1}{1+\alpha\mu}\right)^{2}\|{\mathbf{x}}^{i_{2}}-{\mathbf{x}}^{*}\|$
$\displaystyle\cdots\left(\frac{1}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\|{\mathbf{x}}^{i_{[\frac{k+\tau+1}{\tau}]}}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq\left(\frac{1}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|$

for all $k\in\\{1,2,\cdots,[\frac{n}{\tau}]\\}$ where
$m\in\\{i_{1},i_{2},\cdots,i_{[\frac{k+\tau+1}{\tau}]}\\}$ Hence, we prove the
desired result. ∎

According to Theorem B.4 we can have the following corollary:

######  Corollary B.3.

Let $f({\mathbf{x}})$ be a $\mu$-strongly convex and continuously
differentiable function. Then, the proximal point method

$\displaystyle{\mathbf{x}}^{k+1}=\underset{{\mathbf{x}}}{\mathbf{argmin}}\left(f({\mathbf{x}})+\frac{1}{2\alpha}\|{\mathbf{x}}-{\mathbf{x}}^{k}\|^{2}\right)$

converges linearly

$\|{\mathbf{x}}^{k}-{\mathbf{x}}^{*}\|\leq\left(\frac{1}{1+\alpha\mu}\right)^{k}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$

with any step size $\alpha>0$.

######  Remark B.5.

Although the above theorem only shows the case that $\xi_{i}\geq 0$ we can
still have linear convergence when $\xi_{i}<0$ by letting
$\alpha\geq\frac{\sum_{i=1}^{\tau}\xi_{i}-1}{\mu}$. See the following result.

######  Theorem B.4.

Let $f({\mathbf{x}})$ be a $\mu$-strongly convex and continuously
differentiable function. Suppose that we minimize $f({\mathbf{x}})$ using a
$\tau$steps backward differentiation formula on the gradient flow $\tau$-BDM
where we assume $\sum_{i=1}^{n}\xi_{i}=1$. Then the backward gradient descent
method will converge linearly

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq\left(\frac{\sum_{i=1}^{\tau}|\xi_{i}|}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|.$

###### Proof.

Similar to B.2, we will have:

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq\frac{1}{1+\alpha\mu}\|\sum_{i=1}^{\tau}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{*})\|$
(10)
$\displaystyle\leq\frac{\sum_{i=1}^{\tau}|\xi_{i}|}{1+\alpha\mu}\sum_{i=1}^{\tau}\frac{|\xi_{i}|}{\sum_{i=1}^{\tau}|\xi_{i}|}\|{\mathbf{x}}^{k+i}-{\mathbf{x}}^{*}\|$
(11)

Now by Lemma B.1 and similar to B.2, we will have the desired result.

∎

######  Remark B.6.

Although the above result gives a linear convergence rate when the step size
is large enough, (LABEL:update_rule) is in general hard to solve to
completion, and is usually approximated in practice. Below, we give the linear
convergence rate when (LABEL:update_rule) is approximated.

######  Theorem B.5.

Let $f({\mathbf{x}})$ be a $\mu$-strongly convex and continuously
differentiable function. Suppose that the exact solution to $\tau$ steps of
the backward differential formula is

$\displaystyle\tilde{{\mathbf{x}}}^{k+\tau+1}=\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\alpha\nabla
f(\tilde{{\mathbf{x}}}^{k+\tau+1})$ (12)

where $\sum_{i=1}^{n}\xi_{i}=1$ and $\xi\geq 0$. Suppose that an approximate
solution ${\mathbf{x}}^{k+\tau+1}$ of $\tilde{{\mathbf{x}}}^{k+\tau+1}$
satisfies a contraction :

$\|{\mathbf{x}}^{k+\tau+1}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|\leq\gamma\|\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|$
(13)

where $\gamma<1$.Then ${\mathbf{x}}^{k+\tau+1}$ converges to
${\mathbf{x}}^{*}$ linaerly

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq\left(\gamma+\frac{1+\gamma}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|$

when $\gamma<\frac{\alpha\mu}{2+\alpha\mu}$.

###### Proof.

order seems a bit illogical, how about

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq$
$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|+\|\tilde{{\mathbf{x}}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\overset{{\color[rgb]{1,0,0}\eqref{eq:approx_contract_assp}}}{\leq}$
$\displaystyle\gamma\|\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|+\|\tilde{{\mathbf{x}}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\leq$
$\displaystyle\gamma\left(\|\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-{\mathbf{x}}^{*}\|+\|\tilde{{\mathbf{x}}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\right)+\|\tilde{{\mathbf{x}}}^{k+\tau+1}-{\mathbf{x}}^{*}\|$
$\displaystyle\overset{(a)}{\leq}$
$\displaystyle(\gamma+\frac{1+\gamma}{1+\alpha\mu})\sum_{i=1}^{\tau}\xi_{i}\|{\mathbf{x}}^{k+i}-{\mathbf{x}}^{*}\|$

where (a) follows from (8)

Now by LemmaB.1, we have

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq\left(\gamma+\frac{1+\gamma}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|$

for all $k\in\\{1,2,\cdots,[\frac{n}{\tau}]\\}$. Hence, we prove the desired
result. ∎

######  Corollary B.6.

Let $f({\mathbf{x}})$ be a $\mu$-strongly convex and $L$-smooth function.
Suppose that we use a $\tau$-step backward differential method with update
($\tau$-BDM),where ($\tau$-BDM) is computed approximately by using $n$ steps
of gradient descent with step size $\beta\leq\frac{\alpha}{\alpha L+1}$. This
method will converge linearly when $\alpha$ is large enough with rate

$\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{*}\|\leq\left((1-\beta\mu-\frac{\beta}{\alpha})^{n}+\frac{1+(1-\beta\mu-\frac{\beta}{\alpha})^{n}}{1+\alpha\mu}\right)^{[\frac{k+\tau+1}{\tau}]}\max_{j\in\\{0,1,\cdots\tau-1\\}}\|{\mathbf{x}}^{j}-{\mathbf{x}}^{*}\|$

###### Proof.

Recall that for a $\rho$-strongly convex function $g$, then $n$ steps of
gradient descent with appropriate step size $\beta$ gives convergence
guarantee

$\|{\mathbf{x}}^{n}-{\mathbf{x}}^{*}\|\leq(1-\rho\beta)^{n}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|.$

Now, note that

$g({\mathbf{x}})=\frac{1}{2\alpha}\|{\mathbf{x}}-\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}\|^{2}+f({\mathbf{x}})$

is $\rho=\mu+\frac{1}{\alpha}$ strongly convex function. Then, using $n$ steps
of gradient descent with starting point
$\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}$ gives

$\displaystyle\|{\mathbf{x}}^{k+\tau+1}_{n}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|\leq(1-\beta\mu-\frac{\beta}{\alpha})^{n}\|\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|$

where $\tilde{{\mathbf{x}}}^{k+\tau+1}$ fully minimizes $g$. As a result, (15)
is satisfied with $\gamma=(1-\beta\mu-\frac{\beta}{\alpha})^{n}$, and by
invoking Theorem B.5, we will have the desired result. ∎

The following 2 theorems are for the case when $f({\mathbf{x}})$ is $L-$
weakly convex for both exact and inexact case. We will also need the following
simple fact:

######  Fact B.7.

$\|\sum_{i=1}^{n}a_{i}\|^{2}\leq n\sum_{i=1}^{n}\|a_{i}\|^{2}$

######  Theorem B.8.

Suppose that $f({\mathbf{x}})$ is an $L-$ weakly convex function. Then for any
fixed step size $\alpha>0$, the update rule from ($\tau$-BDM), given
$\xi_{i}$, will have:

$\displaystyle\min_{0\leq s\leq k}\|\nabla
f({\mathbf{x}}^{k+\tau+1})\|=O(\frac{1}{\sqrt{k}})$

if $\alpha<\frac{2-2(\tau-1)\delta}{L}$ and $\delta<\frac{1}{\tau-1}$ where

$\delta=(\tau-1)\sum_{j=1}^{\tau-1}\left(\sum_{i=1}^{j}(\tau-i)\xi_{i}^{2}\right).$

###### Proof.

By $L$-smoothness, we will have:

$\displaystyle f({\mathbf{x}}^{k+\tau})$ $\displaystyle\geq
f({\mathbf{x}}^{k+\tau+1})+\langle\nabla
f({\mathbf{x}}^{k+\tau+1}),{\mathbf{x}}^{k+\tau}-{\mathbf{x}}^{k+\tau+1}\rangle-\frac{L}{2}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle=f({\mathbf{x}}^{k+\tau+1})+\frac{1}{\alpha}\langle\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau+1},{\mathbf{x}}^{k+\tau}-{\mathbf{x}}^{k+\tau+1}\rangle-\frac{L}{2}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle\overset{(*)}{=}f({\mathbf{x}}^{k+\tau+1})+\frac{1}{\alpha}\underbrace{\|{\mathbf{x}}^{k+\tau+1}-\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}\|^{2}}_{=\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}-\frac{1}{\alpha}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2})$

where (*) is $2a^{T}b=\|a\|^{2}+\|b\|^{2}-\|a-b\|^{2}$. Telescoping,

$f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*})\geq\alpha\sum_{k=0}^{T}\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\sum_{k=0}^{T}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}-\frac{1}{\alpha}\sum_{k=0}^{T}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2}$

Now from Fact B.7, notice that:

$\displaystyle\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2}$
$\displaystyle\leq(\tau-1)\sum_{i=1}^{\tau-1}\xi_{i}^{2}\|{\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle\overset{(*)}{\leq}(\tau-1)\sum_{i=1}^{\tau-1}\xi_{i}^{2}\|\sum_{j=i}^{\tau-1}{\mathbf{x}}^{(k+j)}-{\mathbf{x}}^{(k+j-1)}\|^{2}$
$\displaystyle\overset{\text{Fact
\ref{fact:variance}}}{\leq}(\tau-1)\sum_{i=1}^{\tau-1}(\tau-i)\sum_{j=i}^{\tau-1}\xi_{i}^{2}\|{\mathbf{x}}^{k+j}-{\mathbf{x}}^{k+j+1}\|^{2}$
$\displaystyle=(\tau-1)\sum_{j=1}^{\tau-1}\left(\sum_{i=1}^{j}(\tau-i)\xi_{i}^{2}\right)\|{\mathbf{x}}^{k+j}-{\mathbf{x}}^{k+j+1}\|^{2}$

where (*) is from reverse telescoping

${\mathbf{x}}^{(k+i)}-{\mathbf{x}}^{(k+\tau)}=\sum_{j=i}^{\tau-1}({\mathbf{x}}^{(k+j)}-{\mathbf{x}}^{(k+j-1)})$

Hence for $c_{j}=\sum_{i=1}^{j}(\tau-i)\xi_{i}^{2}$ we have:

$\displaystyle\sum_{s=0}^{k}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{s+i}-{\mathbf{x}}^{s+\tau})\|^{2}$
$\displaystyle\leq(\tau-1)\sum_{r=0}^{k}\sum_{j=1}^{\tau-1}c_{j}\|{\mathbf{x}}^{r+j}-{\mathbf{x}}^{r+j+1}\|^{2}$
$\displaystyle=(\tau-1)\sum_{j=1}^{\tau-1}c_{j}\sum_{r=0}^{k}\|{\mathbf{x}}^{r+j}-{\mathbf{x}}^{r+j+1}\|^{2}$
$\displaystyle\overset{s=r+j}{=}(\tau-1)\sum_{j=1}^{\tau-1}c_{j}\sum_{s=j}^{j+k}\|{\mathbf{x}}^{s}-{\mathbf{x}}^{s+1}\|^{2}$
$\displaystyle\leq\underbrace{(\tau-1)\sum_{j=1}^{\tau-1}c_{j}}_{=\delta}\sum_{s=0}^{k+\tau-1}\|{\mathbf{x}}^{s+1}-{\mathbf{x}}^{s}\|^{2}$

As a result, we will have:

$f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*})\geq\alpha\sum_{k=0}^{T}\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}-\frac{\delta}{\alpha}\right)\sum_{k=0}^{T}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}-\frac{\delta}{\alpha}\sum_{s=0}^{\tau-1}\|{\mathbf{x}}^{s+1}-{\mathbf{x}}^{s}\|^{2}$

Now suppose we pick $\alpha,\delta$ such that
$\frac{1}{\alpha}-\frac{L}{2}-\frac{\delta}{\alpha}\geq 0$ (e.g.
$\alpha<\frac{2-2(\tau-1)\delta}{L}$) then

$\frac{1}{T}\sum_{k=0}^{T}\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}\leq\frac{1}{T\alpha}(f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*}))+\frac{\delta}{T\alpha^{2}}\sum_{s=0}^{\tau-1}\|{\mathbf{x}}^{s+1}-{\mathbf{x}}^{s}\|^{2}=O(1/T).$

and we have the desired result. ∎

######  Theorem B.9.

Suppose that $f({\mathbf{x}})$ is an $L$-weakly convex function. Then for any
fixed step size $\alpha>0$,suppose that the exact solution to $\tau$ steps of
backward differential formula is

$\displaystyle\tilde{{\mathbf{x}}}^{k+\tau+1}=\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}-\alpha\nabla
f(\tilde{{\mathbf{x}}}^{k+\tau+1})$ (14)

where $\sum_{i=1}^{n}\xi_{i}=1$ and $\xi\geq 0$. Suppose that an approximate
solution ${\mathbf{x}}^{k+\tau+1}$ of $\tilde{{\mathbf{x}}}^{k+\tau+1}$
satisfies the following contraction with respect to ${\mathbf{x}}^{k+\tau}$:

$\|{\mathbf{x}}^{k+\tau+1}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|\leq\gamma\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|$
(15)

Then, we will have:

$\displaystyle\min_{0\leq s\leq k}\|\nabla
f({\mathbf{x}}^{k+\tau+1})\|=O(\frac{1}{\sqrt{k}})$

if $\alpha<\min({\frac{2-2(\tau-1)\delta-16\gamma^{2}}{L},\frac{1}{L}})$ and
$\delta<\frac{1}{2\tau-2}$ where

$\delta=(\tau-1)\sum_{j=1}^{\tau-1}\left(\sum_{i=1}^{j}(\tau-1-i)\xi_{i}^{2}\right).$

###### Proof.

By $L$-smoothness, we will have:

$\displaystyle f({\mathbf{x}}^{k+\tau})$ $\displaystyle\geq
f({\mathbf{x}}^{k+\tau+1})+\frac{1}{\alpha}\langle\ \alpha\nabla
f({\mathbf{x}}^{k+\tau+1}),{\mathbf{x}}^{k+\tau}-{\mathbf{x}}^{k+\tau+1}\rangle-\frac{L}{2}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle\overset{(*)}{=}f({\mathbf{x}}^{k+\tau+1})+\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle-\frac{1}{\alpha}\|{\mathbf{x}}^{k+\tau}-{\mathbf{x}}^{k+\tau+1}-\alpha\nabla
f({\mathbf{x}}^{k+\tau+1})\|^{2})$
$\displaystyle=f({\mathbf{x}}^{k+\tau+1})+\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle-\frac{1}{\alpha}\|{\mathbf{x}}^{k+\tau}-{\mathbf{x}}^{k+\tau+1}-\alpha\nabla
f({\mathbf{x}}^{k+\tau+1})+\tilde{{\mathbf{x}}}^{k+\tau+1}-\sum_{i=1}^{\tau}\xi_{i}{\mathbf{x}}^{k+i}+\alpha\nabla
f(\tilde{{\mathbf{x}}}^{k+\tau+1})\|^{2}$ $\displaystyle\geq
f({\mathbf{x}}^{k+\tau+1})+\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle-\frac{2}{\alpha}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2}-\frac{2}{\alpha}\|{\mathbf{x}}^{k+\tau+1}+\alpha\nabla
f({\mathbf{x}}^{k+\tau+1})-\tilde{{\mathbf{x}}}^{k+\tau+1}-\alpha\nabla
f(\tilde{{\mathbf{x}}}^{k+\tau+1})\|^{2}$ $\displaystyle\geq
f({\mathbf{x}}^{k+\tau+1})+\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle-\frac{2}{\alpha}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2}-\frac{2(1+\alpha
L)^{2}}{\alpha}\|{\mathbf{x}}^{k+\tau+1}-\tilde{{\mathbf{x}}}^{k+\tau+1}\|^{2}$
$\displaystyle\geq f({\mathbf{x}}^{k+\tau+1})+\alpha\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}\right)\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$
$\displaystyle-\frac{2}{\alpha}\|\sum_{i=1}^{\tau-1}\xi_{i}({\mathbf{x}}^{k+i}-{\mathbf{x}}^{k+\tau})\|^{2}-\frac{2\gamma^{2}(1+\alpha
L)^{2}}{\alpha}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}$

Now we repeat the step from B.8, we will have:

$f({\mathbf{x}}^{\tau})-f({\mathbf{x}}^{*})\geq\alpha\sum_{k=0}^{T}\|\nabla
f({\mathbf{x}}^{(k+\tau+1)})\|^{2}+\left(\frac{1}{\alpha}-\frac{L}{2}-\frac{2\gamma^{2}(1+\alpha
L)^{2}}{\alpha}-\frac{2\delta}{\alpha}\right)\sum_{k=0}^{T}\|{\mathbf{x}}^{k+\tau+1}-{\mathbf{x}}^{k+\tau}\|^{2}-\frac{2\delta}{\alpha}\sum_{s=0}^{\tau-1}\|{\mathbf{x}}^{s+1}-{\mathbf{x}}^{s}\|^{2}$

Hence, we prove the desired result. ∎