# Stochastic Gradient Descent with Nonlinear Conjugate Gradient-Style Adaptive
Momentum

Bao Wang
Department of Mathematics
Scientific Computing and Imaging Institute
University of Utah
<EMAIL_ADDRESS>
Qiang Ye
Department of Mathematics
University of Kentucky
<EMAIL_ADDRESS>
###### Abstract

Momentum plays a crucial role in stochastic gradient-based optimization
algorithms for accelerating or improving training deep neural networks (DNNs).
In deep learning practice, the momentum is usually weighted by a well-
calibrated constant. However, tuning hyperparameters for momentum can be a
significant computational burden. In this paper, we propose a novel _adaptive
momentum_ for improving DNNs training; this adaptive momentum, with no
momentum related hyperparameter required, is motivated by the nonlinear
conjugate gradient (NCG) method. Stochastic gradient descent (SGD) with this
new adaptive momentum eliminates the need for the momentum hyperparameter
calibration, allows a significantly larger learning rate, accelerates DNN
training, and improves final accuracy and robustness of the trained DNNs. For
instance, SGD with this adaptive momentum reduces classification errors for
training ResNet110 for CIFAR10 and CIFAR100 from $5.25\%$ to $4.64\%$ and
$23.75\%$ to $20.03\%$, respectively. Furthermore, SGD with the new adaptive
momentum also benefits adversarial training and improves adversarial
robustness of the trained DNNs.

##  1 Introduction

Given a training dataset $\Omega_{N}:=\\{{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ with
${\bm{x}}_{i}$ and $y_{i}$ being the data-label pair of the $i$th instance.
Natural training, i.e., training a machine learning classifier
$y=g({\bm{x}},{\bm{w}})$ for clean image classification, can be formulated as
solving the following empirical risk minimization (ERM) problem [42]:

${\min_{{\bm{w}}\in{\mathbb{R}}^{d}}f({\bm{w}}):=\frac{1}{N}\sum_{i=1}^{N}f_{i}({\bm{w}}):=\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}(g({\bm{x}}_{i},{\bm{w}}),y_{i}),}$
(1)

where $\mathcal{L}$ is typically the cross-entropy loss between the predicted
label $g({\bm{x}}_{i},{\bm{w}})$ and the ground truth $y_{i}$. To train a
robust classifier under adversarial attacks [15, 27], we often solve the
following empirical adversarial risk minimization (EARM) problem:

${\min_{{\bm{w}}\in{\mathbb{R}}^{d}}\frac{1}{N}\sum_{i=1}^{N}\max_{\|{\bm{x}}_{i}-{\bm{x}}_{i}^{\prime}\|_{2}\leq\epsilon}\mathcal{L}(g({\bm{x}}_{i}^{\prime},{\bm{w}}),y_{i}),}$
(2)

with $\epsilon>0$ being a constant. Training deep neural networks (DNNs) by
solving (1) or (2) is a very difficult task: 1) the objective function is
highly nonconvex [25]; 2) $N$ is very large, e.g., in ImageNet classification
$N\sim 10^{6}$ [37], which makes computing the gradient of the loss function
difficult and inefficient; 3) the dimension of ${\bm{w}}$ is very high; for
instance, in training ResNet200 for ImageNet classification, ${\bm{w}}$ is of
dimension $\sim 65M$ [18]. Due to the above challenges, stochastic gradient
descent (SGD) becomes the method of choice for training DNNs for image
classification [4]; momentum scaled by a well-calibrated weight is usually
integrated with SGD to accelerate or improve training DNNs [34, 40, 3, 31,
44].

Starting from ${\bm{w}}_{0}\in{\mathbb{R}}^{d}$ and ${\bm{p}}_{0}=0$, in the
$n(n\geq 1)$th iteration of SGD with momentum (scaled by a constant $\beta\geq
0$), we randomly sample a mini-batch $\\{i_{k}\\}_{k=1}^{m}\subset[N]\ (m\ll
N)$; update ${\bm{w}}^{n}$ as follows [31]:

$\displaystyle{\begin{aligned}
{\bm{p}}_{n}&=\beta{\bm{p}}_{n-1}+\frac{1}{m}\sum_{k=1}^{m}\nabla
f_{i_{k}}({\bm{w}}_{n}),\\\
{\bm{w}}_{n+1}&={\bm{w}}_{n}-\alpha{\bm{p}}_{n},\end{aligned}}$ (3)

where $\alpha>0$ is the step size. $\nabla f_{i_{k}}({\bm{w}}_{n})$ in (3) can
be replaced with $\nabla f_{i_{k}}({\bm{w}}_{n}-\alpha{\bm{p}}_{n-1})$ to get
the Nesterov momentum [40]. In training DNNs, the tuning of $\alpha$ and
$\beta$ is time consuming. As a result, several adaptive learning rate
algorithms have been developed and are widely used: Adagrad [10] adapts
learning rate to the parameters based on the sum of the squares of the
gradients; Adadelta [45] and RMSprop [21] modify Adagrad by restricting the
window of accumulated past gradients to some fixed size; Adam integrates
momentum with adaptive step size and achieves remarkable performances in many
applications [22, 9, 26, 36]. However, there have not been many works on the
varying or adaptive momentum ($\beta$).

Nonlinear conjugate gradient (NCG) method [11] can be considered as an
adaptive momentum method combined with steepest descent along the search
direction. It significantly accelerates convergence of the gradient descent
method and it has some nice theoretical convergence guarantees [2, 12, 7, 16,
35, 47]. A major obstacle to using NCG is the need for a line search at each
iteration; even an inexact search requires several function/derivative
evaluations. As a result, it is not used very often in DNN applications [24].
However, the success of adaptive momentum in NCG and the success of SGD for
DNN training motivate us to use the adaptive momentum from NCG to improve SGD
based DNN training.

###  1.1 Our Contributions

In this paper, we leverage the celebrated NCG-style adaptive momentum to
accelerate GD/SGD; in particular, to improve training DNNs for image
classification. We propose a GD/SGD method with a fixed learning rate but
leverage an adaptive momentum as used in NCG. We will present some convergence
results to show global convergence under certain conditions on the learning
rate. For the case of quadratic functions, we will show the accelerated
convergence rate under some quite general conditions. We summarize the major
advantages of SGD with adaptive momentum below:

  * •

It converges faster and allows us to use significantly larger step sizes to
train DNNs.

  * •

It improves the accuracy and adversarial robustness of the trained DNNs for
image classification. For instance, it reduces test errors of training
ResNet110 for CIFAR10 and CIFAR100 classification from $5.25\%$ to $4.64\%$
and $23.75\%$ to $20.03\%$, respectively. Furthermore, SGD with the new
adaptive momentum also benefits adversarial training and hence improves
adversarial robustness of the trained DNNs.

  * •

It eliminates the work for momentum-related hyperparameter tuning with almost
no computational overhead.

###  1.2 Related Works

Non-constant momentum has been used to accelerate GD. One of the most exciting
results is the Nesterov accelerated gradient (NAG) [30, 29], which replaces
the constant momentum with an iteration dependent momentum and achieves a
convergence rate of $O(1/k^{2})$ for convex optimization (vs. GD with a
convergence rate $O(1/k)$). However, directly applying NAG to SGD suffers from
error accumulation [8, 44], which can be alleviated by using NAG with
scheduled restart [44] at the cost of hyperparameter calibration.

NCG [11] is a popular optimization method that has been studied extensively.
For various formulations, it has been proved for a general function to have a
descent property and global convergence under some assumptions on the step
size known as Wolf conditions; see [2, 7, 12, 16, 35, 47]. NCG has been
applied to deep learning; [24] empirically compares NCG, L-BFGS, and the
momentum methods and found each to be superior in some problems. There are
some related works in avoiding the line search in NCG. For example, [28] uses
some estimate of the Hessian to approximate optimal step size $\alpha$, while
[1] replaces the line search by starting with some initial $\alpha$ and then
increasing or decreasing its value at each iteration by checking whether the
loss decreases or increases.

###  1.3 Organization

We organize this paper as follows: In Section Algorithm: GD/SGD with Adaptive
Momentum, we give a brief review of NCG, and then we leverage the adaptive
momentum in NCG to accelerate GD and SGD. In Section 3, we give the
convergence guarantees of the proposed algorithms. We then present experiments
to demonstrate the performance of the SGD with adaptive momentum for training
DNNs in Section 4. We end with some concluding remarks in the last section.
Technical proofs and more experimental details are provided in the appendix.

###  1.4 Notations

We denote scalars by lower or upper case letters; vectors/ matrices by
lower/upper case bold face letters. For a vector
${\bm{x}}=(x_{1},\cdots,x_{d})^{T}\in\mathbb{R}^{d}$, we use
$\|{\bm{x}}\|={(\sum_{i=1}^{d}|x_{i}|^{2})^{1/2}}$ to denote its $\ell_{2}$
norm, and the $\ell_{\infty}$ norm of ${\bm{x}}$ by
$\|{\bm{x}}\|_{\infty}=\max_{i=1}^{d}|x_{i}|$. For a matrix ${\mathbf{A}}$, we
use $\|{\mathbf{A}}\|_{2/\infty}$ to denote its induced norm by the vector
$\ell_{2/\infty}$ norm. We denote the set $\\{1,2,\cdots,N\\}$ as $[N]$. For a
function $f({\bm{x}}):\mathbb{R}^{d}\rightarrow\mathbb{R}$, we denote $\nabla
f({\bm{x}})$ and $\nabla^{2}f({\bm{x}})$ the gradient and Hessian of
$f({\bm{x}})$, respectively.

##  2 Algorithm: GD/SGD with Adaptive Momentum

###  2.1 Nonlinear Conjugate Gradient Methods

The classical GD with an optimal learning rate (step size) has a local
convergence rate that depends on the condition number $\kappa$ of the Hessian
matrix at a local minimum. The conjugate gradient (CG) method augments the
gradient with a suitable momentum term as the search direction. In the
quadratic case
($\min_{\bm{w}}f({\bm{w}}\in{\mathbb{R}}^{d}):=1/2{\bm{w}}^{T}{\mathbf{A}}{\bm{w}}+{\bm{b}}^{T}{\bm{w}}$
with ${\mathbf{A}}\in{\mathbb{R}}^{d\times d}$ and
${\bm{b}}\in{\mathbb{R}}^{d}$ being a known matrix and vector, respectively.),
the modified directions maintain orthogonality in the ${\mathbf{A}}$-inner
product (${\bm{w}}^{T}{\mathbf{A}}{\bm{w}}$), and it allows a significantly
accelerated convergence rate that depends on $\sqrt{\kappa}$. It has been
generalized to general nonlinear functions as follows.

  * •

In the first iteration, perform a line search along ${\bm{p}}_{0}:=\nabla
f({\bm{w}}_{0})$ to get the initial step size, i.e.,
$\alpha_{0}:=\arg\min_{\alpha}f({\bm{w}}_{0}-\alpha{\bm{p}}_{0})$, and update
${\bm{w}}$ by ${\bm{w}}_{1}={\bm{w}}_{0}-\alpha_{0}{\bm{p}}_{0}$.

  * •

For $n\geq 1$, we perform the following updates for the $n$th iteration:

    * –

Compute

$\beta_{n}=\beta_{n}^{FR}:=\frac{(\nabla f({\bm{w}}_{n})^{T}\nabla
f({\bm{w}}_{n}))}{(\nabla f({\bm{w}}_{n-1})^{T}\nabla f({\bm{w}}_{n-1}))}.$

    * –

Update the search direction:

${\bm{p}}_{n}=\nabla f({\bm{w}}_{n})+\beta_{n}{\bm{p}}_{n-1}.$

    * –

Perform a line search:

$\alpha_{n}=\arg\min_{\alpha}f({\bm{w}}_{n}-\alpha{\bm{p}}_{n}).$

    * –

Update the position:

${\bm{w}}_{n+1}={\bm{w}}_{n}-\alpha_{n}{\bm{p}}_{n}.$

There are several possible formulations on $\beta_{n}$ in the literature, and
the one we present $\beta_{n}^{FR}$ is known as the Fletcher-Reeves formula
[11]; see [20, 7, 32, 39, 16] for other formulations and related theoretical
properties. The NCG has been empirically found to have some similar
convergence properties as the classical linear CG method. There have been
several analyses to show a descent property and convergence of NCG for a
general function under some forms of the Wolf conditions for inexact line
search of $\alpha_{n}$; see [2, 12, 7, 16, 35, 47]. However, there appears to
be no result characterizing its CG-like accelerated convergence rate.

In the NCG method, a line search is performed to determine $\alpha_{n}$, while
$\beta_{n}$ can be regarded as the momentum coefficient. Here, the momentum
leverages the past gradient instead of the past position, which is different
from the momentum in (3). This has an advantage over the traditional momentum
method in that the momentum coefficient is adaptively determined and no tuning
is needed. A disadvantage of NCG is that even an inexact line search for
$\alpha_{n}$ requires several function/gradient evaluations and would make the
method less appealing for training DNNs. Therefore, NCG is rarely used for
DNNS.

###  2.2 (Stochastic) Gradient Descent with NCG Momentum

The discussions in the previous parts motivate us to integrate the adaptive
momentum coefficient $\beta_{n}$ into GD/SGD. We propose to consider GD/SGD
with an adaptive momentum, i.e., with a fixed $\alpha$ but with
$\beta=\beta_{n}^{FR}$ as the momentum coefficient at each step. This has two
potential benefits. As a generalization of NCG, this may preserve some
convergence properties of NCG. As a momentum method, there is no need to
determine or tune the hyperparameter for momentum. We have chosen the
Fletcher-Reeves formula $\beta_{n}^{FR}$ for its simplicity and robustness, as
indicated by our preliminary numerical testing. We call the resulting
algorithms FRGD/FRSGD, which we state as follows. Starting with
${\bm{w}}_{0}$, we set ${\bm{p}}_{-1}=\mathbf{0}$ and $\beta_{0}=0$ and
iterate for $n\geq 0$ as follows:

$\displaystyle{\bm{p}}_{n}=\bm{r}_{n}+\beta_{n}{\bm{p}}_{n-1},$ (4)
$\displaystyle{\bm{w}}_{n+1}={\bm{w}}_{n}-\alpha{\bm{p}}_{n},$

where

  * •

FRGD:

$\bm{r}_{n}=\nabla f({\bm{w}}_{n})$

and

$\beta_{n}=(\bm{r}_{n}^{T}\bm{r}_{n})/(\bm{r}_{n-1}^{T}\bm{r}_{n-1});$

  * •

FRSGD:

$\bm{r}_{n}=\frac{1}{m}{\sum_{j=1}^{m}}\nabla f_{i_{j}}({\bm{w}}_{n})$

and

$\beta_{n}=(\bm{r}_{n}^{T}\bm{r}_{n})/(\bm{r}_{n-1}^{T}\bm{r}_{n-1})$

.

We have found this adaptive momentum method significantly accelerates the
convergence of GD with momentum and outperforms NAG as well. The only extra
computational cost over the momentum method is in computing an inner product
$\bm{r}_{n}^{T}\bm{r}_{n}$ at each step, which is negligible.

Before testing its performance in training DNNs, we first present an academic
example to illustrate its potential advantage.

Example 2.1: We consider the following quadratic optimization problem [17]:

${\small\min_{\bm{w}}f({\bm{w}})=\frac{1}{2}{\bm{w}}^{T}\mathbf{L}{\bm{w}}-{\bm{w}}^{T}{\bm{b}},}$
(5)

where $\mathbf{L}\in\mathbb{R}^{500\times 500}$ is the Laplacian of a cycle
graph, and ${\bm{b}}$ is a $500$-dimensional vector whose first entry is $1$
and all the other entries are $0$. It is easy to see that $f({\bm{w}})$ is
convex (not strongly convex) with Lipschitz constant $4$. We run GD, GD with
momentum scaled by $0.9$ (GD $+$ Momentum), NAG, and FRGD with step size $1/4$
(the same hyperparameters as that used in [17]). As shown in Fig. 1, GD $+$
Momentum converges faster than GD, while NAG speeds up GD $+$ Momentum
dramatically and converges to the minimum in an oscillatory fashion. More
interestingly, FRGD converges exponentially fast and significantly outperforms
all other methods in this case.

---
Figure 1: Comparison between a few optimization algorithms for optimizing the
quadratic function
$f({\bm{w}})=1/2{\bm{w}}^{T}{\mathbf{L}}{\bm{w}}+{\bm{b}}^{T}{\bm{w}},\
{\bm{w}},{\bm{b}}\in{\mathbb{R}}^{500},\ \mbox{and}\
{\mathbf{L}}\in{\mathbb{R}}^{500\times 500}$, where ${\mathbf{L}}$ is the
Laplacian of a cycle graph, and ${\bm{b}}$ is vector with the first entry
being $1$ and all the others are $0$s [17]. Momentum accelerates GD slightly;
NAG oscillates to the minimum, ${\bm{w}}^{*}$, and converges much faster than
GD (with momentum); FRSGD converge almost exponentially fast to
${\bm{w}}^{*}$.

This example demonstrates that FRGD converges at a rate much faster than GD.
We will present in the next section some theoretical results to demonstrate
this property.

##  3 Main Theory

In this section, we present some convergence results to demonstrate descent
property and convergence of FRGD (4) under some conditions on $\alpha$. We
shall focus on strongly convex functions and quadratic functions. Our results
are applicable only locally for a general function. For the case of quadratic
functions, we present convergence bounds to demonstrate accelerated
convergence rate. The proofs of these results will be given in Appendix A. We
first consider a convex function $f({\bm{w}})$ with Lipschitz continuous
Hessian matrix.

At each FRGD iteration, there exists
${\bm{y}}_{i}=(1-\xi_{i})\bm{w}_{n}+\xi_{i}\bm{w}_{n+1}$ for some
$\xi_{i}\in[0,1]$ ($1\leq i\leq d$) s.t.

${\nabla f(\bm{w}_{n+1})-\nabla
f(\bm{w}_{n})={\mathbf{H}}_{n}(\bm{w}_{n+1}-\bm{w}_{n}),}$ (6)

where ${\mathbf{H}}_{n}=\left[\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}({\bm{y}}_{i})\right]_{i,j=1}^{d}$; see Lemma 1 in the Appendix A. This
${\mathbf{H}}_{n}$ is approximately a Hessian matrix.

######  Theorem 1.

Consider the adaptive momentum method FRGD (4) for
$f(\bm{w}):{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$. Assume that
${\mathbf{H}}_{ij}(\bm{w}):=\frac{\partial^{2}f(\bm{w})}{\partial
w_{i}\partial w_{j}}$ is Lipschitz continuous with the Lipschiz constant $C$,
i.e. $\|{\mathbf{H}}_{ij}(\bm{w})-{\mathbf{H}}_{ij}(\tilde{\bm{w}})\|\leq
C\|\bm{w}-\tilde{\bm{w}}\|$ for all $1\leq i,j\leq d$. Let ${\mathbf{H}}_{n}$
be as defined in (6). Assume for some $K>0$ that the eigenvalues of
${\mathbf{H}}_{n}$ are on $[\lambda_{\min},\lambda_{\max}]$ with
$\lambda_{\min}>0$ for $n\leq K$ and
$\alpha\leq\frac{\lambda_{\min}}{(\lambda_{\max}^{2}+2Cd\|\bm{r}_{0}\|)K^{2}}$.
Then

${{\bm{p}}_{n}^{T}\bm{r}_{n}>0}$

and

$\|\bm{r}_{n}\|\leq\sqrt{1-\alpha\lambda_{\min}}\|\bm{r}_{n-1}\|,$

where $\bm{r}_{n}=\nabla f(\bm{w}_{n})$ and $n\leq K$.

The theorem shows that ${\bm{p}}_{n}$ is a descent direction and $\bm{r}_{n}$
converges monotonically with a rate of at least
$\sqrt{1-\alpha\lambda_{\min}}$. Although such properties are expected for GD,
it is important that with the adaptive momentum, FRGD maintains these
properties. However, NAG does not converge monotonically to the minimum;
instead it oscillates. Moreover, it is worth noting that: 1) our assumption on
the Hessian is different from the convergence theory of GD, we need Lipschitz
Hessian while GD requires bounded Hessian; 2) both the step size constraint
and the convergence rate, $\sqrt{1-\alpha\lambda_{min}}$, are also different
from that of GD.

###  3.1 Quadratic Functions

To further study the convergence rate, we consider a quadratic function
$f({\bm{w}})=\frac{1}{2}{\bm{w}}^{T}{\mathbf{A}}{\bm{w}}-{\bm{b}}^{T}{\bm{w}}$,
where ${\mathbf{A}}$ is a symmetric positive definite matrix. In this case,
the problem reduces to a linear system ${\mathbf{A}}{\bm{w}}={\bm{b}}$. Then
the classical CG method has been studied extensively and strong convergence
results exist. With a fixed $\alpha$, most properties that the analysis of CG
relies on no longer hold. Fortunately, a techniques used for the analysis of
inexact CG due to round off errors or inexact preconditioning [13, 41] can be
adapted to our method. Using coupled two-term recurrences, our adaptive
momentum becomes a Krylov subspace method, and we can derive the following
convergence bound.

######  Theorem 2.

Consider the adaptive momentum method FRGD (4) for
$f({\bm{w}})=\frac{1}{2}{\bm{w}}^{T}{\mathbf{A}}{\bm{w}}-{\bm{b}}^{T}{\bm{w}}$.
Let ${\mathbf{Z}}_{n}=[{\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n-1}]$ where
${\bm{z}}_{i}=\bm{r}_{i}/\|\bm{r}_{i}\|$. If
${\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n}$ are linearly independent, then

${\|\bm{r}_{n}\|\leq
2(1+K_{n})\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{n}\|\bm{r}_{0}\|,}$
(7)

where $K_{n}\leq n(1+n\rho/2)\|{\mathbf{A}}\|\kappa({\mathbf{Z}}_{n+1})$,
$\rho=\max_{0\leq j<i\leq n-1}\|\bm{r}_{i}\|^{2}/\|\bm{r}_{j}\|^{2}$, $\kappa$
is the spectral condition number of ${\mathbf{A}}$, and
$\kappa({\mathbf{Z}}_{n+1})$ is the spectral condition number of
${\mathbf{Z}}_{n+1}$.

The bound in (7) contains a linearly converging term with the rate
$(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$, but it also depends on the term
$K_{n}$, which may grow with $n$. The key term in $K_{n}$ is
$\kappa({\mathbf{Z}}_{n+1})=\|{\mathbf{Z}}_{n+1}\|\|{\mathbf{Z}}_{n+1}^{+}\|\approx\sqrt{n+1}\|{\mathbf{Z}}_{n+1}^{+}\|$,
which measures linear independence among
${\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n}$. Thus, as long as
${\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n}$ does not completely lose
linear independence, $K_{n}$ may be a modestly increasing term so that
$\|\bm{r}_{n}\|$ converges at a rate close to
$(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$. Note that $\rho$ may be expected to be
bounded. In particular, if $\|\bm{r}_{n}\|$ is monotonic, which holds under
the condition of Theorem 1, then $\rho\leq 1$.

Note that the classical CG method converges at the rate of
$(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$ and so does the momentum method with the
following optimal $\alpha$ and $\beta$:

${\alpha=\frac{4}{(\sqrt{\lambda_{\max}}+\sqrt{\lambda_{\min}})^{2}},}$
$\beta=\frac{(\sqrt{\lambda_{\max}}-\sqrt{\lambda_{\min}})^{2}}{(\sqrt{\lambda_{\max}}+\sqrt{\lambda_{\min}})^{2}},$

where $\lambda_{\max}$ and $\lambda_{\min}$ are the largest and the smallest
eigenvalues of $A$ respectively; see [33]. However, they require strong
conditions with the former requiring variable $\alpha_{n}$ (see its definition
in the previous section), while the latter requiring the optimal $\alpha$ and
$\beta$. All these methods significantly accelerates GD with constant step
size $\alpha$, which has at best a convergence rate of
$(\kappa-1)/(\kappa+1)$.

We remark that a slightly more general bound
$\|\bm{r}_{n}\|\leq(1+K_{n})\min_{p\in{\cal
P}_{n},p(0)=1}\|p({\mathbf{A}})\bm{r}_{0}\|$ holds in place of (7) without the
assumption that ${\mathbf{A}}$ is symmetric positive definite; see the proof
in the Appendix. This bound and thus our method are applicable to the more
general situations of positive semi-definite ${\mathbf{A}}$ (Example 2.1) or
the trust region problem with indefinite ${\mathbf{A}}$ [6, 46].

We also note that the theorem also _does not require any explicit condition on
the learning rate $\alpha$_. Although $\alpha$ may affect the quality of the
basis ${\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n}$ generated, as long as
$K_{n}$ increases gradually at a rate slower than
$(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$, we have convergence of $\bm{r}_{n}$.
This may explain the success of our method with quite large learning rates
(see our numerical results in the Experiment section). Of course, this is only
true to the extent that $\alpha$ is not so large that the condition number of
the basis generated grows unbounded.

##  4 Experiments

In this section, we present numerical results to illustrate the advantage of
FRSGD over the baseline SGD with constant momentum in training DNNs for image
classification. We run all experiments with five independent random seeds and
report the means and standard deviations of the results.

|
---|---
Figure 2: Plots of epochs vs. training loss and accuracy of PreResNet56
trained by FRSGD, SGD (with momentum), and SGD $+$ NM with different learning
rate. FRSGD is the most robust optimizer to large learning rate; where the
performance of SGD and SGD $+$ NM deteriorate when a large learning rate is
used.

#### Objective.

Our experimental results will demonstrate the following advantages of FRSGD
over the baseline methods of SGD with momentum or Nesterov momentum: 1) FRSGD
converges significantly faster; 2) FRSGD is significantly more robust to large
step sizes; 3) DNN trained by FRSGD is more accurate and more adversarially
robust than that trained by the baseline methods.

#### Datasets.

We consider CIFAR10 & CIFAR100 datasets [23] due to limited computing
resources. Both datasets consist of 60K 32 by 32 color images with 50K/10K
training/test split. CIFAR10/CIFAR100 contains 10/100 different classes, with
each class having the same number of images.

#### Tasks, Experimental Settings, and Baselines.

We consider both natural and adversarial training, by solving (1) and (2),
respectively. We use pre-activated ResNets (PreResNets) models of different
depths [19]. As baseline methods, we consider SGD with the standard momentum
and with the Nesterov momentum scaled by $0.9$, and we denote them as SGD and
SGD $+$ NM, respectively, in the following context. We note that SGD is the
optimizer used in the original ResNet implementations [18, 19]. For the SGD
and SGD $+$ NM baselines, we follow the standard-setting of ResNets by running
it for $200$ epochs with an initial learning rate of $0.1$ and decay it by a
factor of $10$ at the $80$th, $120$th, and $160$th epoch, respectively. For
FRSGD, we run $240$ epochs with an initial learning rate of $0.5$ and reduce
it by a factor of $10$ at the $180$th, $220$th, and $230$th epoch,
respectively 111Here, we are able to use a larger learning rate under which
FRSGD is still stable, and we decay the learning rate when the decrease of the
loss function becomes slower. . For adversarial training, we use the same SGD
and FRSGD solvers described above to solve the outer minimization problem, and
we run $10$ iterations of the iterative fast gradient sign method (IFGSM10)
attack with $\alpha=2/255$ and $\epsilon=8/255$ to approximate the solution of
the inner maximization problem. We provide the details of IFGSM [14] and a few
other attacks in the Appendix.

###  4.1 FRSGD is Robust Under Large Learning Rates

In this subsection, we compare the performance of SGD, SGD $+$ NM, and FRSGD
in training PreResNet56 for CIFAR10 classification using different learning
rates. We set the learning rate to be 0.1, 0.3, and 0.5, with all the other
parameters the same as before. We plot epochs vs. training loss and training
accuracy in Fig. 2. These results show that: 1) under the same small learning
rate, e.g., 0.1, FRSGD converges remarkably faster than both SGD and SGD $+$
NM; 2) the convergence of SGD and SGD $+$ NM gets deteriorated severely when a
larger learning rate is used; in particular, the training loss will not
converge when 0.5 is used as the learning rate. However, FRSGD maintains
convergence even when a very large learning rate is used; 3) the training loss
and accuracy curves of SGD and SGD $+$ NM get plateau very quickly at a large
loss, while FRSGD continues to decay.

###  4.2 FRSGD Improves Accuracy of DNNs for Image Classification

#### CIFAR10.

We consider training PreResNets with different depths using the settings
mentioned before for the CIFAR10 classification. We list the test errors of
different ResNets trained by different stochastic optimization algorithms in
Table 1. In general, SGD performs on par with SGD $+$ NM; SGD $+$ NM has small
advantages over SGD for training shallow DNNs. FRSGD outperforms SGD by
$0.5\sim 0.7$% for ResNets with the depth ranging from 56 to 470. These
improvements over already small error rates are significant in the relative
sense, e.g., for PreResNet470, the relative error reduction is $\sim 13$%
($4.92$% vs. $4.27$%).

Table 1: Test error (%) on CIFAR10 using the SGD (with momentum), SGD $+$ NM, and FRSGD. We also include the reported results from [19] (in parentheses) in addition to our reproduced results. ResNets trained by FRSGD is consistently more accurate than those trained by SGD and SGD $+$ NM. Network | SGD (baseline) | SGD$+$NM | FRSGD | Improve over SGD | Improve over SGD$+$NM
---|---|---|---|---|---
PreResNet56 | $6.12\pm 0.24$ | $5.90\pm 0.17$ | $\boldsymbol{5.39\pm 0.13}$ | $0.73$ | $0.51$
PreResNet110 |  $5.25\pm 0.14$ ($6.37$) | $5.24\pm 0.16$ | $4.73\pm 0.12$ | $0.52$ | $0.51$
PreResNet164 |  $5.10\pm 0.19$ ($5.46$) | $5.08\pm 0.21$ | $4.50\pm 0.16$ | $0.60$ | $0.58$
PreResNet290 | $5.05\pm 0.23$ | $5.04\pm 0.12$ | $4.44\pm 0.10$ | $0.61$ | $0.60$
PreResNet470 | $4.92\pm 0.10$ | $4.97\pm 0.15$ | $4.27\pm 0.09$ | $0.65$ | $0.70$

#### CIFAR100.

Here, we consider CIFAR100 classification with the same DNNs and the same
settings as those used for CIFAR10 classification. We report the test errors
in Table 2. In this case, FRSGD improves the test accuracy over both SGD and
SGD $+$ NM by $\sim 1.0$ to $1.6$%.

Table 2: Test error (%) on CIFAR100 using the SGD (with momentum), SGD $+$ NM, and FRSGD. We also include the reported results from [19] (in parentheses) in addition to our reproduced results. ResNets trained by FRSGD is uniformly more accurate than those trained by SGD and SGD $+$ NM. Network | SGD (baseline) | SGD$+$NM | FRSGD | Improve over SGD | Improve over SGD$+$NM
---|---|---|---|---|---
PreResNet56 | $26.60\pm 0.33$ | $26.14\pm 0.38$ | $\boldsymbol{25.00\pm 0.32}$ | $1.60$ | $1.14$
PreResNet110 | $23.75\pm 0.20$ | $23.65\pm 0.36$ | $22.52\pm 0.35$ | $1.23$ | $1.13$
PreResNet164 |  $22.76\pm 0.37$ ($24.33$) | $22.79\pm 0.29$ | $21.38\pm 0.34$ | $1.38$ | $1.41$
PreResNet290 | $21.78\pm 0.21$ | $21.68\pm 0.21$ | $20.66\pm 0.31$ | $1.12$ | $1.02$
PreResNet470 | $21.43\pm 0.30$ | $21.21\pm 0.30$ | $19.92\pm 0.29$ | $1.51$ | $1.29$
Table 3: Lists of the optimal training/test loss and accuracy of PreResNet110 trained by FRSGD, SGD (with momentum), SGD $+$ NM, and Adam with $240$ epochs. Adam has smaller training loss than the others, but the PreResNet110 trained by FRSGD has the smallest test loss/error. Optimizer |  Training Loss |  Training Error Rate (%) |  Test Loss |  Test Error rate (%)
---|---|---|---|---
SGD | $0.00529\pm 0.00043$ | $0.042\pm 0.006$ | $0.1950\pm 0.00091$ | $5.23\pm 0.15$
SGD $+$ NM | $0.00462\pm 0.00047$ | $0.032\pm 0.005$ | $0.1846\pm 0.00101$ | $5.19\pm 0.16$
Adam | $\boldsymbol{0.00033\pm 0.00003}$ | $\boldsymbol{0.002\pm 0.002}$ | $0.3237\pm 0.00125$ | $6.47\pm 0.31$
FRSGD | $0.00680\pm 0.00021$ | $0.090\pm 0.002$ | $\boldsymbol{0.1611\pm 0.00086}$ | $\boldsymbol{4.73\pm 0.12}$

#### FRSGD vs. Adam and SGD with More Epochs.

In the original ResNet experiments, we have run SGD and SGD $+$ NM for 200
epochs, after which no training loss decay is observed. To compare over longer
iterations, we train PreResNet110 by running SGD, SGD $+$ NM, and FRSGD for
240 epochs 222Based on trial and error, we found that adding 20 epochs each in
the first and second learning rate stages gives the best performance. All the
reported results are based on using this setting.. Moreover, we also compare
them with running Adam [22] for 240 epochs using a well-calibrated initial
learning rate 0.003 and the same decaying schedule as SGD with 240 epochs, and
we use the default value for all Adam’s other hyperparameters. Table 3 lists
the training and testing losses as well as errors of PreResNet110 trained by
different optimizers on CIFAR10. The best test error of SGD, SGD $+$ NM, and
Adam for CIFAR10 classification are $5.23\pm 0.15$%, $5.19\pm 0.16$%, and
$6.47\pm 0.31$% respectively, compared with $4.73\pm 0.12$ for FRSGD. We see
that adding 40 more epochs to SGD and SGD $+$ NM does not improve
classification accuracy much, and this is because the training loss has
reached the plateau at each stage with a budget of 200 epochs. Adam converges
faster with a smaller final training loss, but the testing loss and accuracy
are far behind those of SGD, SGD $+$ NM, and FRSGD.

Table 4:  Test accuracy (%) of PreResNet110 on CIFAR10 using PGD adversarial training with SGD $+$ NM and FRSGD as the outer solver. FRSGD improves accuracies for classifying both clean and adversarial images. Optimizer | Natural | FGSM | IFGSM10 | IFGSM20 | IFGSM40 | IFGSM100 | C&W
---|---|---|---|---|---|---|---
SGD $+$NM | $82.19\pm 0.29$ | $57.61\pm 0.33$ | $55.35\pm 0.42$ | $52.02\pm 0.34$ | $51.45\pm 0.33$ | $51.08\pm 0.35$ | $62.92\pm 0.50$
FRSGD | $82.36\pm 0.27$ | $58.27\pm 0.29$ | $55.83\pm 0.31$ | $53.07\pm 0.28$ | $52.39\pm 0.25$ | $52.15\pm 0.19$ | $63.05\pm 0.33$
Table 5:  Test accuracy (%) of PreResNet110 on CIFAR100 using PGD adversarial training with SGD $+$ NM and FRSGD as the outer solver. FRSGD improves accuracies for classifying both clean and adversarial images. Optimizer | Clean | FGSM | IFGSM10 | IFGSM20 | IFGSM40 | IFGSM100 | C&W
---|---|---|---|---|---|---|---
SGD | $54.75\pm 0.52$ | $30.75\pm 0.41$ | $29.61\pm 0.45$ | $27.87\pm 0.44$ | $27.51\pm 0.42$ | $27.40\pm 0.48$ | $38.97\pm 0.66$
FRSGD | $54.95\pm 0.49$ | $31.77\pm 0.43$ | $30.79\pm 0.33$ | $29.32\pm 0.37$ | $29.09\pm 0.40$ | $29.01\pm 0.39$ | $39.01\pm 0.50$

#### Training FRSGD with a Large Number of Epochs.

In the previous experiments, considering the training efficiency, we limited
the budget for training epochs of FRSGD by dropping the learning rate when the
training loss convergence slows down but before reaching plateaus. This
learning rate reduction may be premature and the result may not be the best
accuracy our method can achieve. In this experiment, we relax this budget and
use a much larger number of epochs for FRSGD and see if we can get more
improvement in classification accuracy. In particular, we train
PreResNet110/PreResNet290 by running 400 epochs of FRSGD with an initial
learning rate of 0.5 and reduce the learning rate by a factor of 10 at the
200th, 300th, and 350th epochs, respectively.

In this setting, we get the best test error rates of $4.64\pm 0.12/4.26\pm
0.09$% for CIFAR10 and $20.03\pm 0.29/19.89\pm 0.19$% for CIFAR100, which
remarkably improves what we get by using 240 epochs ($4.73\pm 0.12/4.44\pm
0.10$ for CIFAR10 and $22.52\pm 0.35/20.66\pm 0.31$ for CIFAR100).
Furthermore, the training loss of FRSGD at the last epoch in training the
PreResNet290 for CIFAR10 classification also becomes significantly smaller
than that of SGD or SGD $+$ NM ($0.00138\pm 0.00012$ (FRSGD), vs. $0.0048\pm
0.0003$ (SGD), and $0.0052\pm 0.0003$ (SGD $+$ NM)).

###  4.3 FRSGD Improves Adversarial Training

Finally, we numerically demonstrate that FRSGD can also improve the
adversarial robustness of the trained DNNs through adversarial training. We
train the PreResNet110 by applying the adversarial training using the settings
listed before. Then we apply the well-trained PreResNet110 to classify the
test set under three kinds of benchmark adversarial attacks: fast gradient
sign method (FGSM), $m$ steps IFGSM (IFGSMm with $m=10,20,40,100$) [14], and
C&W attacks [5]. We apply the same set of hyperparameters for these attacks as
that used in [43, 27] in the following experiments. A brief introduction of
these attacks and the used hyperparameters are available in Appendix.

Tables 4 and 5 list the accuracy of the adversarially trained PreResNet110 for
classifying CIFAR10 and CIFAR100 images with or without adversarial attacks
333We only compare FRSGD with SGD $+$ NM since SGD performance is weaker than
SGD $+$ NM in this case.. First, we see that the robust PreResNet110 trained
by FRSGD is slightly more accurate than that trained by SGD $+$ NM for
classifying the clean CIFAR10 and CIFAR100 images without any attack, e.g.,
the accuracy of FRSGD is $82.36\pm 0.27$% and $54.95\pm 0.49$% for CIFAR10 and
CIFAR100 classification, while the corresponding accuracy of the model trained
by SGD $+$ NM is $82.19\pm 0.29$% and $54.75\pm 0.52$%, respectively. Second,
the model trained by FRSGD is more robust than that trained by SGD $+$ NM
under all the adversarial attacks mentioned before, e.g., under the IFGSM100
attack, the robust accuracy of these two models are $52.15\pm 0.19$% vs.
$51.08\pm 0.35$% for CIFAR10 classification, and are $27.40\pm 0.48$% vs.
$29.01\pm 0.39$% for CIFAR100 classification. Although the improvements are
minor, the good performance of FRSGD in this difficult setting illustrates its
robustness in different problem types.

##  5 Concluding Remarks

In this paper, we leveraged adaptive momentum from the NCG to improve SGD, and
the resulting algorithm performs surprisingly well in the following sense: 1)
It can accelerate GD significantly; in particular, we observed that it
achieves exponential convergence for optimizing a specific convex function; 2)
It allows us to use much larger step sizes and converges faster than SGD with
(Nesterov) momentum in training DNNs; 3) DNNs trained by FRSGD have remarkably
higher classification accuracy and are more robust to adversarial attacks for
image classification. The method is as simple as SGD and is easy to implement.
It is well suited for DNN training.

There are several interesting open problems that are worth further
investigations. First, can we integrate the adaptive momentum with adaptive
step size to further improve stochastic optimization algorithms? Second, can
we prove stronger convergence results for FRGD/FRSGD under more general
conditions? Third, can we leverage adaptive momentum to improve training DNNs
for other deep learning tasks beyond image classification? These will
potentially be our future works.

## References

  * [1] Saurabh Adya, Vinay Palakkode, and Oncel Tuzel.  Nonlinear conjugate gradients for scaling synchronous distributed dnn training, 2018.
  * [2] M. AL-BAALI.  Descent Property and Global Convergence of the Fletcher—Reeves Method with Inexact Line Search.  IMA Journal of Numerical Analysis, 5(1):121–124, 01 1985.
  * [3] Yoshua Bengio, Nicolas Boulanger-Lewandowski, and Razvan Pascanu.  Advances in optimizing recurrent networks.  In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 8624–8628. IEEE, 2013.
  * [4] Léon Bottou, Frank E Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  Siam Review, 60(2):223–311, 2018.
  * [5] N. Carlini and D.A. Wagner.  Towards evaluating the robustness of neural networks.  IEEE European Symposium on Security and Privacy, pages 39–57, 2016\.
  * [6] Yair Carmon and John C Duchi.  Analysis of krylov subspace solutions of regularized non-convex quadratic problems.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 10705–10715. Curran Associates, Inc., 2018.
  * [7] Yu-Hong Dai and Yaxiang Yuan.  A nonlinear conjugate gradient method with a strong global convergence property.  SIAM Journal on optimization, 10(1):177–182, 1999.
  * [8] Olivier Devolder, François Glineur, and Yurii Nesterov.  First-order methods of smooth convex optimization with inexact oracle.  Mathematical Programming, 146(1-2):37–75, 2014.
  * [9] Timothy Dozat.  Incorporating nesterov momentum into adam.  2016\.
  * [10] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of machine learning research, 12(Jul):2121–2159, 2011.
  * [11] Reeves Fletcher and Colin M Reeves.  Function minimization by conjugate gradients.  The computer journal, 7(2):149–154, 1964.
  * [12] Jean Charles Gilbert and Jorge Nocedal.  Global convergence properties of conjugate gradient methods for optimization.  SIAM Journal on Optimization, 2(1):21–42, 1992.
  * [13] Gene H. Golub and Qiang Ye.  Inexact preconditioned conjugate gradient method with inner-outer iteration.  SIAM Journal on Scientific Computing, 21(4):1305–1320, 1999.
  * [14] I. J. Goodfellow, J. Shlens, and C. Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6275, 2014.
  * [15] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572, 2014.
  * [16] William W Hager and Hongchao Zhang.  A new conjugate gradient method with guaranteed descent and an efficient line search.  SIAM Journal on optimization, 16(1):170–192, 2005.
  * [17] Moritz Hardt.  Robustness versus acceleration.  http://blog.mrtz.org/2014/08/18/robustness-versus-acceleration.html, 2014\.
  * [18] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [19] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  In European conference on computer vision, pages 630–645. Springer, 2016.
  * [20] Magnus R Hestenes et al.  Methods of conjugate gradients for solving linear systems.  Journal of research of the National Bureau of Standards, 49(6):409–436, 1952.
  * [21] Geoffrey Hinton, Nitish Srivastava, and Kevin Swersky.  Neural networks for machine learning lecture 6a overview of mini-batch gradient descent.
  * [22] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [23] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [24] Quoc V. Le, Jiquan Ngiam, Adam Coates, Abhik Lahiri, Bobby Prochnow, and Andrew Y. Ng.  On optimization methods for deep learning, 2011.
  * [25] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [26] Ilya Loshchilov and Frank Hutter.  Fixing weight decay regularization in adam.  2018\.
  * [27] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  arXiv preprint arXiv:1706.06083, 2017.
  * [28] Martin Moller.  Moller, m.f.: A scaled conjugate gradient algorithm for fast supervised learning. neural networks 6, 525-533.  Neural Networks, 6:525–533, 12 1993.
  * [29] Yurii Nesterov.  Introductory lectures on convex programming volume i: Basic course.  1998\.
  * [30] Yurii E Nesterov.  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In Dokl. akad. nauk Sssr, volume 269, pages 543–547, 1983.
  * [31] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems, pages 8024–8035, 2019.
  * [32] Elijah Polak and Gerard Ribiere.  Note sur la convergence de méthodes de directions conjuguées.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 3(R1):35–43, 1969.
  * [33] Boris Polyak.  Some methods of speeding up the convergence of iteration methods.  Ussr Computational Mathematics and Mathematical Physics, 4:1–17, 12 1964.
  * [34] Boris T Polyak.  Some methods of speeding up the convergence of iteration methods.  USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
  * [35] Michael James David Powell.  Some convergence properties of the conjugate gradient method.  Mathematical Programming, 11(1):42–49, 1976.
  * [36] Sashank J Reddi, Satyen Kale, and Sanjiv Kumar.  On the convergence of adam and beyond.  arXiv preprint arXiv:1904.09237, 2019.
  * [37] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  International journal of computer vision, 115(3):211–252, 2015\.
  * [38] Yousef Saad.  Iterative Methods for Sparse Linear Systems.  Society for Industrial and Applied Mathematics, second edition, 2003.
  * [39] Jonathan Richard Shewchuk et al.  An introduction to the conjugate gradient method without the agonizing pain, 1994.
  * [40] Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton.  On the importance of initialization and momentum in deep learning.  In International conference on machine learning, pages 1139–1147, 2013.
  * [41] Charles H. Tong and Qiang Ye.  Analysis of the finite precision bi-conjugate gradient algorithm for nonsymmetric linear systems.  Math. Comput., 69:1559–1575, 2000.
  * [42] Vladimir Vapnik.  Principles of risk minimization for learning theory.  In Advances in neural information processing systems, pages 831–838, 1992.
  * [43] B. Wang, B. Yuan, Z. Shi, and S. Osher.  ResNet ensemble via the Feynman-Kac formalism to improve natural and robust acurcies.  In Advances in Neural Information Processing Systems, 2019.
  * [44] Bao Wang, Tan M Nguyen, Andrea L Bertozzi, Richard G Baraniuk, and Stanley J Osher.  Scheduled restart momentum for accelerated stochastic gradient descent.  arXiv preprint arXiv:2002.10583, 2020.
  * [45] Matthew D Zeiler.  Adadelta: an adaptive learning rate method.  arXiv preprint arXiv:1212.5701, 2012.
  * [46] Lei-Hong Zhang, Chungen Shen, and Ren-Cang Li.  On the generalized lanczos trust-region method.  SIAM Journal on Optimization, 27(3):2110–2142, 2017.
  * [47] G. Zoutendijk.  Nonlinear programming, computational methods, in integer and nonlinear programming.  Integer and Nonlinear Programming, J. Abadie, ed., 1970, pp. 37-86.

##  Appendix A Proof of Main Theorems

We first give a lemma and then prove Theorem 1. For convenience, we restate
all the theorems.

######  Lemma 1.

Assume $\frac{\partial^{2}f}{\partial w_{i}\partial w_{j}}$ is Lipschitz
continuous with the Lipschiz constant $C$ for all $1\leq i,j\leq d$, i.e.
$|\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}(\bm{w})-\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}(\tilde{\bm{w}})|\leq C\|\bm{w}-\tilde{\bm{w}}\|$. For each fixed $n\geq
0$, there exists ${\bm{y}}_{i}=(1-\xi_{i})\bm{w}_{n}+\xi_{i}\bm{w}_{n+1}$ for
some $\xi_{i}\in[0,1]$ ($1\leq i\leq d$) s.t.

$\nabla f(\bm{w}_{n+1})-\nabla
f(\bm{w}_{n})={\mathbf{H}}_{n}(\bm{w}_{n+1}-\bm{w}_{n})\;\;\mbox{ where
}\;\;H_{n}=\left[\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}({\bm{y}}_{i})\right]_{i,j=1}^{d}.$

Furthermore,

$\|{\mathbf{H}}_{n+1}-{\mathbf{H}}_{n}\|\leq
Cd(\|\bm{w}_{n+2}-\bm{w}_{n+1}\|+\|\bm{w}_{n+1}-\bm{w}_{n}\|).$ (8)

###### Proof.

Write $\nabla f(\bm{w})=[f_{1}(\bm{w}),\cdots,f_{d}(\bm{w})]^{T}$, where
$f_{i}(\bm{w}):=\frac{\partial f}{\partial w_{i}}$. Then, applying Taylor’s
theorem to each $f_{i}$, we have

$f_{i}(\bm{w}_{n+1})-f_{i}(\bm{w}_{n})=\nabla
f_{i}({\bm{y}}_{i})^{T}(\bm{w}_{n+1}-\bm{w}_{n}),$

where ${\bm{y}}_{i}=(1-\xi_{i})\bm{w}_{n}+\xi_{i}\bm{w}_{n+1}$ for some
$\xi_{i}\in[0,1]$, which is $n$ dependent. Since

${\mathbf{H}}_{n}=\left[\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}({\bm{y}}_{i})\right]_{i,j=1}^{d}=\left(\begin{array}[]{c}\nabla\frac{\partial
f}{\partial w_{1}}({\bm{y}}_{1})^{T}\\\ \nabla\frac{\partial f}{\partial
w_{2}}({\bm{y}}_{2})^{T}\\\ \vdots\\\ \nabla\frac{\partial f}{\partial
w_{d}}({\bm{y}}_{d})^{T}\\\ \end{array}\right)=\left(\begin{array}[]{c}\nabla
f_{1}({\bm{y}}_{1})^{T}\\\ \nabla f_{2}({\bm{y}}_{2})^{T}\\\ \vdots\\\ \nabla
f_{d}({\bm{y}}_{d})^{T}\\\ \end{array}\right),$

we have

$\nabla f(\bm{w}_{n+1})-\nabla
f(\bm{w}_{n})=\left(\begin{array}[]{c}f_{1}(\bm{w}_{n+1})-f_{1}(\bm{w}_{n})\\\
\vdots\\\ f_{d}(\bm{w}_{n+1})-f_{d}(\bm{w}_{n})\\\
\end{array}\right)={\mathbf{H}}_{n}(\bm{w}_{n+1}-\bm{w}_{n}).$

We have ${\mathbf{H}}_{n+1}$ defined similarly and we can write

${\mathbf{H}}_{n+1}=\left[\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}({\bm{z}}_{i})\right]_{i,j=1}^{d}=\left[\nabla
f_{1}({\bm{z}}_{1}),\nabla f_{2}({\bm{z}}_{2}),\cdots,\nabla
f_{d}({\bm{z}}_{d})\right]^{T},$

for some ${\bm{z}}_{i}=(1-\eta_{i})\bm{w}_{n+1}+\eta_{i}\bm{w}_{n+2}$ with
$\eta_{i}\in[0,1]$. Now, by the Lipschitz continuity,

$\displaystyle\|\nabla f_{i}({\bm{z}}_{i})-\nabla
f_{i}({\bm{y}}_{i})\|_{\infty}$ $\displaystyle\leq$ $\displaystyle
C\|{\bm{z}}_{i}-{\bm{y}}_{i}\|$ $\displaystyle=$ $\displaystyle
C\|(1-\eta_{i})\bm{w}_{n+1}+\eta_{i}\bm{w}_{n+2}-(1-\xi_{i})\bm{w}_{n}-\xi_{i}\bm{w}_{n+1}\|$
$\displaystyle=$ $\displaystyle
C\|\eta_{i}(\bm{w}_{n+2}-\bm{w}_{n+1})+(1-\xi_{i})(\bm{w}_{n+1}-\bm{w}_{n})\|$
$\displaystyle\leq$ $\displaystyle
C(\eta_{i}\|\bm{w}_{n+2}-\bm{w}_{n+1}\|+(1-\xi_{i})\|\bm{w}_{n+1}-\bm{w}_{n}\|)$
$\displaystyle\leq$ $\displaystyle
C(\|\bm{w}_{n+2}-\bm{w}_{n+1}\|+\|\bm{w}_{n+1}-\bm{w}_{n}\|).$

Thus
$\|{\mathbf{H}}_{n+1}-{\mathbf{H}}_{n}\|\leq\|{\mathbf{H}}_{n+1}-{\mathbf{H}}_{n}\|_{F}\leq
d\max_{i}\|\nabla f_{i}({\bm{z}}_{i})-\nabla
f_{i}({\bm{y}}_{i})\|_{\infty}\leq
Cd(\|\bm{w}_{n+2}-\bm{w}_{n+1}\|+\|\bm{w}_{n+1}-\bm{w}_{n}\|)$, where
$\|\cdot\|_{F}$ is the Frobenius norm of matrices. ∎

We now prove Theorem 1.

######  Theorem 3 (Theorem 1 Restate).

Consider the adaptive momentum method FRGD (4) for $f({\bm{w}})$. Assume that
${\mathbf{H}}_{ij}({\bm{w}}):=\frac{\partial^{2}f}{\partial w_{i}\partial
w_{j}}$ is Lipschitz continuous with the Lipschiz constant $C$, i.e.
$\|{\mathbf{H}}_{ij}({\bm{w}})-{\mathbf{H}}_{ij}(\tilde{{\bm{w}}})\|\leq
C\|{\bm{w}}-\tilde{{\bm{w}}}\|$ for all $1\leq i,j\leq d$. Let
${\mathbf{H}}_{n}$ be as defined in (6). Assume for some $K>0$ that the
eigenvalues of ${\mathbf{H}}_{n}$ are on $[\lambda_{\min},\lambda_{\max}]$
with $\lambda_{\min}>0$ for $n\leq K$ and
$\alpha\leq\frac{\lambda_{\min}}{(\lambda_{\max}^{2}+2Cd\|\bm{r}_{0}\|)K^{2}}$.
Then

$\bm{r}_{n}^{T}\bm{p}_{n}>0\;\;\mbox{ and
}\;\;\|\bm{r}_{n}\|\leq\sqrt{1-\alpha\lambda_{\min}}\|\bm{r}_{n-1}\|,$

where $\bm{r}_{n}=\nabla f(\bm{w}_{n})$ and $k\leq K$.

###### Proof.

Consider $\bm{w}_{n-1},\bm{w}_{n},\bm{w}_{n+1}$. By Lemma 1, there is
${\mathbf{H}}_{n}$ such that

$\bm{r}_{n+1}-\bm{r}_{n}=\nabla f(\bm{w}_{n+1})-\nabla
f(\bm{w}_{n})={\mathbf{H}}_{n}(\bm{w}_{n+1}-\bm{w}_{n})=-\alpha{\mathbf{H}}_{n}\bm{p}_{n}.$
(9)

Furthermore

$\|{\mathbf{H}}_{n}-{\mathbf{H}}_{n-1}\|\leq
Cd(\|\bm{w}_{n+1}-\bm{w}_{n}\|+\|\bm{w}_{n}-\bm{w}_{n-1}\|)=Cd\alpha(\|\bm{p}_{n}\|+\|\bm{p}_{n-1}\|).$
(10)

Now we prove by induction in $n$ that, for $1\leq n\leq K$,

$\frac{\bm{r}_{n-1}^{T}{\mathbf{H}}_{n-1}\bm{p}_{n-1}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}\geq\lambda_{\min},\;\|\bm{p}_{n-1}\|\leq
n\|\bm{r}_{n-1}\|,\;\mbox{ and
}\;\|\bm{r}_{n}\|\leq\sqrt{1-\alpha\lambda_{\min}}\|\bm{r}_{n-1}\|.$ (11)

First, since $\bm{p}_{0}=\bm{r}_{0}$ and
$\bm{r}_{1}=\bm{r}_{0}-\alpha{\mathbf{H}}_{0}\bm{p}_{0}=\bm{r}_{0}-\alpha{\mathbf{H}}_{0}\bm{r}_{0}$,
we have
$\bm{r}_{0}^{T}{\mathbf{H}}_{0}\bm{p}_{0}=\bm{r}_{0}^{T}{\mathbf{H}}_{0}\bm{r}_{0}\geq\lambda_{\min}\bm{r}_{0}^{T}\bm{r}_{0}$,
$\|\bm{p}_{0}\|=\|\bm{r}_{0}\|$, and

$\displaystyle\bm{r}_{1}^{T}\bm{r}_{1}$ $\displaystyle=$
$\displaystyle\bm{r}_{0}^{T}\bm{r}_{0}-2\alpha\bm{r}_{0}^{T}{\mathbf{H}}_{0}\bm{r}_{0}+\alpha^{2}\bm{r}_{0}^{T}{\mathbf{H}}_{0}^{2}\bm{r}_{0}$
$\displaystyle\leq$
$\displaystyle\bm{r}_{0}^{T}\bm{r}_{0}-2\alpha\lambda_{\min}\bm{r}_{0}^{T}\bm{r}_{0}+\alpha^{2}\lambda_{\max}^{2}\bm{r}_{0}^{T}\bm{r}_{0}$
$\displaystyle\leq$
$\displaystyle(1-\alpha\lambda_{\min})\bm{r}_{0}^{T}\bm{r}_{0},$

where we have used $-\lambda_{\min}+\alpha\lambda_{\max}^{2}\leq 0$. So, (11)
holds for $n=1$.

Assume that (11) holds for some $n\leq K-1$. We prove it for $n+1$. Using (4)
and the induction assumption, we have
$\beta_{n}=\frac{\|\bm{r}_{n}\|^{2}}{\|\bm{r}_{n-1}\|^{2}}\leq 1$ and then

$\displaystyle\bm{p}_{n}^{T}\bm{p}_{n}$ $\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}+2\beta_{n}\bm{r}_{n}^{T}\bm{p}_{n-1}+\beta_{n}^{2}\bm{p}_{n-1}^{T}\bm{p}_{n-1}$
$\displaystyle\leq$
$\displaystyle\|\bm{r}_{n}\|^{2}+2\frac{\|\bm{r}_{n}\|^{2}}{\|\bm{r}_{n-1}\|^{2}}\|\bm{r}_{n}\|\|\bm{p}_{n-1}\|+\frac{\|\bm{r}_{n}\|^{2}}{\|\bm{r}_{n-1}\|^{2}}\|\bm{p}_{n-1}\|^{2}$
$\displaystyle\leq$
$\displaystyle\|\bm{r}_{n}\|^{2}+2\|\bm{r}_{n}\|^{2}n+\|\bm{r}_{n}\|^{2}n^{2}$
$\displaystyle=$ $\displaystyle(1+n)^{2}\|\bm{r}_{n}\|^{2}.$

where we have used $\|\bm{p}_{n-1}\|\leq n\|\bm{r}_{n-1}\|$.

Next, using (4) and (9), we have

$\displaystyle\bm{r}_{n}^{T}{\mathbf{H}}_{n}\bm{p}_{n}$ $\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}{\mathbf{H}}_{n}\bm{r}_{n}+\beta_{n}\bm{r}_{n}^{T}{\mathbf{H}}_{n}\bm{p}_{n-1}$
(12) $\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}{\mathbf{H}}_{n}\bm{r}_{n}+\beta_{n}\bm{r}_{n-1}^{T}{\mathbf{H}}_{n}\bm{p}_{n-1}-\alpha\beta_{n}\bm{p}_{n-1}^{T}{\mathbf{H}}_{n-1}{\mathbf{H}}_{n}\bm{p}_{n-1}$
$\displaystyle\geq$
$\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\frac{\bm{r}_{n-1}^{T}{\mathbf{H}}_{n-1}\bm{p}_{n-1}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}+\bm{r}_{n}^{T}\bm{r}_{n}\frac{\bm{r}_{n-1}^{T}({\mathbf{H}}_{n}-{\mathbf{H}}_{n-1})\bm{p}_{n-1}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}$
$\displaystyle-\alpha\frac{\bm{r}_{n}^{T}\bm{r}_{n}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}\lambda_{\max}^{2}\bm{p}_{n-1}^{T}\bm{p}_{n-1}$
$\displaystyle\geq$
$\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\lambda_{\min}-\bm{r}_{n}^{T}\bm{r}_{n}\frac{\|\bm{p}_{n-1}\|}{\|\bm{r}_{n-1}\|}\|{\mathbf{H}}_{n}-{\mathbf{H}}_{n-1}\|-\alpha\lambda_{\max}^{2}\bm{r}_{n}^{T}\bm{r}_{n}\frac{\|\bm{p}_{n-1}\|^{2}}{\|\bm{r}_{n-1}\|^{2}}$
$\displaystyle\geq$
$\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\lambda_{\min}-\bm{r}_{n}^{T}\bm{r}_{n}nCd\alpha(\|\bm{p}_{n}\|+\|\bm{p}_{n-1}\|)-\alpha\lambda_{\max}^{2}\bm{r}_{n}^{T}\bm{r}_{n}n^{2}$
$\displaystyle\geq$
$\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\lambda_{\min}-\alpha\bm{r}_{n}^{T}\bm{r}_{n}Cdn((n+1)\|\bm{r}_{n}\|+n\|\bm{r}_{n-1}\|)-\alpha\lambda_{\max}^{2}\bm{r}_{n}^{T}\bm{r}_{n}n^{2}$
$\displaystyle\geq$
$\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\lambda_{\min}-\alpha\bm{r}_{n}^{T}\bm{r}_{n}Cd2K^{2}\|\bm{r}_{0}\|-\alpha\lambda_{\max}^{2}\bm{r}_{n}^{T}\bm{r}_{n}K^{2}$
$\displaystyle\geq$ $\displaystyle\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n},$

where the last inequality follows from the condition on $\alpha$.

Finally, using the two inequalities above, we have

$\displaystyle\bm{r}_{n+1}^{T}\bm{r}_{n+1}$ $\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}-2\alpha\bm{r}_{n}^{T}{\mathbf{H}}_{n}\bm{p}_{n}+\alpha^{2}\bm{p}_{n}^{T}{\mathbf{H}}_{n}^{2}\bm{p}_{n}$
$\displaystyle\leq$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}-2\alpha\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\alpha^{2}\lambda_{\max}^{2}\bm{p}_{n}^{T}\bm{p}_{n}$
$\displaystyle\leq$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}-2\alpha\lambda_{\min}\bm{r}_{n}^{T}\bm{r}_{n}+\alpha^{2}\lambda_{\max}^{2}(n+1)^{2}\bm{r}_{n}^{T}\bm{r}_{n}$
$\displaystyle\leq$
$\displaystyle(1-\alpha\lambda_{\min})\bm{r}_{n}^{T}\bm{r}_{n}-\alpha(\lambda_{\min}-\alpha\lambda_{\max}^{2}(n+1)^{2})\bm{r}_{n}^{T}\bm{r}_{n}$
$\displaystyle\leq$
$\displaystyle(1-\alpha\lambda_{\min})\bm{r}_{n}^{T}\bm{r}_{n},$

where we note that $n+1\leq K$ and hence
$\lambda_{\min}-\alpha\lambda_{\max}^{2}(n+1)^{2}\geq 0$. This completes the
proof of (11).

We now prove $\bm{r}_{n}^{T}\bm{p}_{n}\geq\bm{r}_{n}^{T}\bm{r}_{n}>0$ by
induction. The case $n=0$ is trivial and assume it hold for $n-1$. Then

$\displaystyle\bm{r}_{n}^{T}\bm{p}_{n}$ $\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}+\beta_{n}\bm{r}_{n}^{T}\bm{p}_{n-1}$
$\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}+\beta_{n}\bm{r}_{n-1}^{T}\bm{p}_{n-1}-\alpha\beta_{n}\bm{p}_{n-1}^{T}H_{n-1}\bm{p}_{n-1}$
$\displaystyle=$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}\frac{\bm{r}_{n-1}^{T}\bm{p}_{n-1}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}-\alpha\bm{r}_{n}^{T}\bm{r}_{n}\frac{\bm{p}_{n-1}^{T}H_{n-1}\bm{p}_{n-1}}{\bm{r}_{n-1}^{T}\bm{r}_{n-1}}$
$\displaystyle\geq$
$\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}+\bm{r}_{n}^{T}\bm{r}_{n}-\alpha\lambda_{\max}\bm{r}_{n}^{T}\bm{r}_{n}\frac{\|\bm{p}_{n-1}\|^{2}}{\|\bm{r}_{n-1}\|^{2}}$
$\displaystyle\geq$ $\displaystyle
2\bm{r}_{n}^{T}\bm{r}_{n}-\alpha\lambda_{\max}n^{2}\bm{r}_{n}^{T}\bm{r}_{n}$
$\displaystyle\geq$ $\displaystyle\bm{r}_{n}^{T}\bm{r}_{n}$

where the last inequality follows from
$\alpha\leq\frac{\lambda_{\min}}{(\lambda_{\max}^{2}+2Cd\|\bm{r}_{0}\|)K^{2}}<\frac{\lambda_{\min}}{\lambda_{\max}^{2}K^{2}}\leq\frac{1}{\lambda_{\max}K^{2}}$.
This completes the proof the theorem. ∎

Next, we present the proof of Theorem 2.

######  Theorem 4 (Theorem 2 Restate).

Consider the CG-momentum method (4) for
$f({\bm{w}})=\frac{1}{2}{\bm{w}}^{T}{\mathbf{A}}{\bm{w}}-{\bm{b}}^{T}{\bm{w}}$.
Let ${\mathbf{Z}}_{n}=[{\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n-1}]$ where
${\bm{z}}_{i}=\bm{r}_{i}/\|\bm{r}_{i}\|$. If
${\bm{z}}_{0},{\bm{z}}_{1},\cdots,{\bm{z}}_{n}$ are linearly independent, then

$\|\bm{r}_{n}\|\leq
2(1+K_{n})\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{n}\|\bm{r}_{0}\|,$
(13)

where $K_{n}\leq n(1+\frac{n}{2}\rho)\|A\|\kappa({\mathbf{Z}}_{n+1})$,
$\rho=\max_{0\leq j<i\leq n-1}(\|\bm{r}_{i}\|^{2})/(\|\bm{r}_{j}\|^{2})$, and
$\kappa$ is the spectral condition number of $A$ and
$\kappa({\mathbf{Z}}_{n+1})$ is the spectral condition number of
${\mathbf{Z}}_{n+1}$.

###### Proof.

Let ${\mathbf{R}}_{n}=[\bm{r}_{0},\bm{r}_{1},\cdots,\bm{r}_{n-1}]$,
${\mathbf{D}}_{n}=\mbox{diag}\\{\|\bm{r}_{0}\|,\|\bm{r}_{1}\|,\cdots,\|\bm{r}_{n-1}\|\\}$
and ${\mathbf{P}}_{n}=[\bm{p}_{0},\cdots,\bm{p}_{n-1}]$. Then
${\mathbf{Z}}_{n}={\mathbf{R}}_{n}{\mathbf{D}}_{n}^{-1}$. Using
$\bm{r}_{k+1}=\bm{r}_{k}-\alpha{\mathbf{A}}\bm{p}_{k}$, we have

$\alpha{\mathbf{A}}{\mathbf{P}}_{n}=[\bm{r}_{0}-\bm{r}_{1},\bm{r}_{1}-\bm{r}_{2},\cdots,\bm{r}_{n-1}-\bm{r}_{n}]={\mathbf{R}}_{n}{\mathbf{L}}_{n}-\bm{r}_{n}e_{n}^{T}={\mathbf{Z}}_{n}{\mathbf{D}}_{n}{\mathbf{L}}_{n}-\bm{r}_{n}{\bm{e}}_{n}^{T},$

and using $\bm{p}_{k}=\bm{r}_{k}+\beta_{k}\bm{p}_{k-1}$

${\mathbf{Z}}_{n}={\mathbf{R}}_{n}{\mathbf{D}}_{n}^{-1}={\mathbf{P}}_{n}{\mathbf{U}}_{n}{\mathbf{D}}_{n}^{-1},$

where ${\bm{e}}_{n}=[0,\cdots,0,1]^{T}$ and ${\mathbf{L}}_{n}$ is the $n\times
n$ lower bidiagonal matrix with $1$ on the diagonal and $-1$ on the
subdiagonal, and ${\mathbf{U}}_{n}$ is the upper bidiagonal matrix with 1 on
the diagonal and $-\beta_{1},\cdots,-\beta_{n-1}$ on the superdiagonal.
Combining the two equations, we obtain

${\mathbf{A}}{\mathbf{Z}}_{n}={\mathbf{Z}}_{n}{{\mathbf{T}}}_{n}-\frac{1}{\alpha^{\prime}}\frac{\bm{r}_{n}}{\|\bm{r}_{0}\|}{\bm{e}}_{n}^{T},$

where
${\mathbf{T}}_{n}=\frac{1}{\alpha}{\mathbf{D}}_{n}{\mathbf{L}}_{n}{\mathbf{U}}_{n}{\mathbf{D}}_{n}^{-1}$
and $\alpha^{\prime}=\alpha\|\bm{r}_{n-1}\|/\|\bm{r}_{0}\|.$ Note that
$\alpha^{\prime}={\bm{e}}_{n}^{T}{\mathbf{T}}_{n}^{-1}{\bm{e}}_{1}$. Apply
Theorem 3.5 of [41] (with $\hat{\Delta}_{n}$ there equal to 0 and the indexes
shifted by 1) to the above equation, we have

$\|\bm{r}_{n}\|\leq(1+K_{n})\min_{p\in{\cal
P}_{n},p(0)=1}\|p({\mathbf{A}})\bm{r}_{0}\|,$ (14)

where
$K_{n}=\|{\mathbf{A}}{\mathbf{Z}}_{n}{\mathbf{T}}_{n}^{-1}[{\mathbf{I}}_{n}\;\;0]{\mathbf{Z}}_{n+1}^{+}\|\leq\|{\mathbf{A}}\|\|{\mathbf{T}}_{n}^{-1}\|\|{\mathbf{Z}}_{n}\|\|{\mathbf{Z}}_{n+1}^{+}\|$
and ${\cal P}_{n}$ is the set of polynomials of degree $n$.

Note that $\beta_{k}=\frac{\|\bm{r}_{k}\|^{2}}{\|\bm{r}_{k-1}\|^{2}}$. Write
${\mathbf{T}}_{n}=\frac{1}{\alpha}{\mathbf{D}}_{n}{\mathbf{L}}_{n}{\mathbf{D}}_{n}^{-1}{\mathbf{D}}_{n}{\mathbf{U}}_{n}{\mathbf{D}}_{n}^{-1}=\frac{1}{\alpha}\hat{{\mathbf{L}}}_{n}\hat{{\mathbf{U}}}_{n}$,
where

$\displaystyle\hat{{\mathbf{L}}}_{n}:={\mathbf{D}}_{n}{\mathbf{L}}_{n}{\mathbf{D}}_{n}^{-1}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccc}1&{}\hfil&{}\hfil&{}\hfil&{}\hfil\\\
-\frac{\|\bm{r}_{1}\|}{\|\bm{r}_{0}\|}&1&{}\hfil&{}\hfil&{}\hfil\\\
{}\hfil&-\frac{\|\bm{r}_{2}\|}{\|\bm{r}_{1}\|}&\ddots&{}\hfil&{}\hfil\\\
&{}\hfil&\ddots&1&{}\hfil\\\
&{}\hfil&{}\hfil&-\frac{\|\bm{r}_{n-1}\|}{\|\bm{r}_{n-2}\|}&1\end{array}\right)$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}1&{}\hfil&{}\hfil&{}\hfil&{}\hfil\\\
-\sqrt{\beta_{1}}&1&{}\hfil&{}\hfil&{}\hfil\\\
{}\hfil&-\sqrt{\beta_{2}}&\ddots&{}\hfil&{}\hfil\\\
{}\hfil&{}\hfil&\ddots&1&{}\hfil\\\
{}\hfil&{}\hfil&{}\hfil&-\sqrt{\beta_{n-1}}&1\end{array}\right),$

and

$\displaystyle\hat{{\mathbf{U}}}_{n}:={\mathbf{D}}_{n}{\mathbf{U}}_{n}{\mathbf{D}}_{n}^{-1}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}1&-\beta_{1}\frac{\|\bm{r}_{0}\|}{\|\bm{r}_{1}\|}&{}\hfil&{}\hfil&{}\hfil\\\
{}\hfil&1&-\beta_{2}\frac{\|\bm{r}_{1}\|}{\|\bm{r}_{2}\|}&{}\hfil&{}\hfil\\\
{}\hfil&{}\hfil&\ddots&\ddots&{}\hfil\\\
{}\hfil&{}\hfil&{}\hfil&1&-\beta_{n-1}\frac{\|\bm{r}_{n-2}\|}{\|\bm{r}_{n-1}\|}\\\
{}\hfil&{}\hfil&{}\hfil&{}\hfil&1\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}1&-\sqrt{\beta_{1}}&{}\hfil&{}\hfil&{}\hfil\\\
{}\hfil&1&-\sqrt{\beta_{2}}&{}\hfil&{}\hfil\\\
{}\hfil&{}\hfil&\ddots&\ddots&{}\hfil\\\
{}\hfil&{}\hfil&{}\hfil&1&-\sqrt{\beta_{n-1}}\\\
{}\hfil&{}\hfil&{}\hfil&{}\hfil&1\end{array}\right)=\hat{{\mathbf{L}}}_{n}^{T}.$

Then $\hat{{\mathbf{L}}}_{n}^{-1}$ is a lower triangular matrix with the
diagonals being 1 and with the $(i,j)$ entry being
$\sqrt{\beta_{j}\beta_{j+1}\cdots\beta_{i-1}}=\|\bm{r}_{i-1}\|/\|\bm{r}_{j-1}\|$
for $i>j$. Then bounding $\|\bm{r}_{i-1}\|^{2}/\|\bm{r}_{j-1}\|^{2}$ by
$\rho$, we have $\|\hat{{\mathbf{L}}}_{n}^{-1}\|_{F}^{2}\leq n+n(n-1)\rho/2$.
So,
$\|{\mathbf{T}}_{n}^{-1}\|=\alpha\|\hat{{\mathbf{L}}}_{n}^{-1}\|^{2}\leq\alpha\|\hat{{\mathbf{L}}}_{n}^{-1}\|_{F}^{2}\leq\alpha
n(1+n\rho/2)$. Combining this with
$\|{\mathbf{Z}}_{n}\|\|{\mathbf{Z}}_{n+1}^{+}\|\leq\|{\mathbf{Z}}_{n+1}\|\|{\mathbf{Z}}_{n+1}\|=\kappa({\mathbf{Z}}_{n+1})$
results in $K_{n}\leq
n\alpha(1+\frac{n}{2}\rho)\|{\mathbf{A}}\|\kappa({\mathbf{Z}}_{n+1})$.
Finally, the bound follows from the standard CG convergence bound [38, p.215]
that shows

$\min_{p\in{\cal
P}_{n},p(0)=1}\|p({\mathbf{A}})\bm{r}_{0}\|\leq\min_{p\in{\cal
P}_{n},p(0)=1}\max_{i}|p(\lambda_{i})|\,\|\bm{r}_{0}\|\leq
2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{n}\|\bm{r}_{0}\|,$

where $\lambda_{i}$ (for $1\leq i\leq d$) are eigenvalues of ${\mathbf{A}}$. ∎

##  Appendix B Adversarial Attacks

We focus on the $\ell_{\infty}$ norm-based FGSM, IFGSM, and C&W white-box
attacks. For a given image-label pair $\\{{\bm{x}},y\\}$, a given ML model
$g({\bm{x}},{\bm{w}})$, and the associated loss
$f({\bm{x}},y):=f(g({\bm{x}},{\bm{w}}),y)$:

  * •

Fast gradient sign method (FGSM) searches an adversarial, ${\bm{x}}^{\prime}$,
within an $\ell_{\infty}$-ball as

${\footnotesize{\bm{x}}^{\prime}={\bm{x}}+\epsilon\cdot{\rm
sign}\left(\nabla_{\bm{x}}f({\bm{x}},y)\right),}$

and we set $\epsilon=8/255$ in all of our experiments.

  * •

Iterative FGSM (IFGSMM) [14] iterates FGSM and clip the range as

${\footnotesize{\bm{x}}^{(m)}={\rm
Clip}_{{\bm{x}},\epsilon}\left\\{{\bm{x}}^{(m-1)}+\alpha\cdot{\rm
sign}\left(\nabla_{{\bm{x}}^{(m-1)}}f({\bm{x}}^{(m-1)},y)\right)\right\\},\
\mbox{w/}\ {\bm{x}}^{(0)}={\bm{x}},\ m=1,\cdots,M,}$

and we set $\epsilon=8/255$ and $\alpha=1/255$ in IFGSM attacks with different
number of iterations.

  * •

C&W attack [5] searches the minimal perturbation ($\delta$) attack as

${\footnotesize\min_{\bm{\delta}}||\bm{\delta}||_{\infty},\ \ \mbox{subject
to}\ \ g({\bm{w}},{\bm{x}}+\bm{\delta})=t,\;{\bm{x}}+\bm{\delta}\in[0,1]^{d},\
\mbox{for}\ \forall t\neq y,}$

we use the same setting as that used in [43] for C&W attack.