# An Efficient Nonlinear Acceleration method that Exploits Symmetry of the
Hessian

Huan He
Harvard University
<EMAIL_ADDRESS>
&Shifan Zhao
Emory University
<EMAIL_ADDRESS>
&Ziyuan Tang
University of Minnesota
<EMAIL_ADDRESS>
Joyce C Ho
Emory University
<EMAIL_ADDRESS>
&Yousef Saad
University of Minnesota
<EMAIL_ADDRESS>
&Yuanzhe Xi
Emory University
<EMAIL_ADDRESS>

###### Abstract

Nonlinear acceleration methods are powerful techniques to speed up fixed-point
iterations. However, many acceleration methods require storing a large number
of previous iterates and this can become impractical if computational
resources are limited. In this paper, we propose a nonlinear Truncated
Generalized Conjugate Residual method (nlTGCR) whose goal is to exploit the
symmetry of the Hessian to reduce memory usage. The proposed method can be
interpreted as either an inexact Newton or a quasi-Newton method. We show
that, with the help of global strategies like residual check techniques,
nlTGCR can converge globally for general nonlinear problems and that under
mild conditions, nlTGCR is able to achieve superlinear convergence. We further
analyze the convergence of nlTGCR in a stochastic setting. Numerical results
demonstrate the superiority of nlTGCR when compared with several other
competitive baseline approaches on a few problems. Our code will be available
in the future.

##  1 Introduction

In this paper, we consider solving the fixed-point problem:

$\text{ Find }x\in\mathbb{R}^{n}\text{ such that }x=H(x)\text{. }$ (1)

This problem has received a surge of interest due to its wide range of
applications in mathematics, computational science and engineering. Most
optimization algorithms are iterative, and their goal is to find a related
fixed-point of the form (1), where $H:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$
is the iteration mapping which can potentially be nonsmooth or noncontractive.
When the optimization problem is convex, $H$ is typically nonexpansive, and
the solution set of the fixed-point problem is the same as that of the
original optimization problem, or closely related to it. Consider the simple
fixed-point iteration $x_{k+1}=H(x_{k})$ which produces a sequence of iterates
$\\{x_{0},x_{1},\cdots,x_{K}\\}$. When the iteration converges, its limit is a
fixed-point, i.e., $x^{\ast}=H(x^{\ast})$. However, an issue with fixed-point
iteration is that it does not always converge, and when it does, it might
reach the limit very slowly.

To address this issue, a number of acceleration methods have been proposed and
studied over the years, such as the reduced-rank extrapolation (RRE) [58],
minimal-polynomial extrapolation (MPE) [12], modified MPE (MMPE) [33], and the
vector $\epsilon$-algorithms [8]. Besides these algorithms, Anderson
Acceleration (AA) [1] has received enormous recent attention due to its nice
properties and its success in machine learnin applications [56, 22, 57, 14,
60, 44, 25, 63]. In practice, since computing the Hessian of the objective
function is commonly difficult or even unavailable, AA can be seen as a
practical alternative to Newton’s method [34]. Also, compared with the
classical iterative methods such as the nonlinear conjugate gradient (CG)
method [24], no line-search or trust-region technique is performed in AA, and
this is a big advantage in large-scale unconstrained optimization.
Empirically, it is observed that AA is quite successful in accelerating
convergence. We refer readers to [9] for a recent survey of acceleration
methods.

However, classical AA has one undesirable disadvantage in that it is expensive
in terms of memory as well as computational cost, especially in a nonconvex
stochastic setting, where only sublinear convergence can be expected when only
stochastic gradients can be accessed [3]. In light of this, a number of
variants of AA have been proposed which aim at improving its performance and
robustness (e.g., [39, 63, 62, 56, 67]). The above-cited works focus on
improving the convergence behavior of AA, but they do not consider reducing
the memory cost. In machine learning, we often encounter practical situations
where the number of parameters is quite large and for this reason, it is not
practical to use a large number of vectors in the acceleration methods. It is
not clear whether or not the symmetric structure of the Hessian can be
exploited in a scheme like AA to reduce the memory cost while still
maintaining the convergence guarantees. In this paper, we will demonstrate how
this can be accomplished with a new algorithm that is superior to AA in
practice.

Our contributions. This paper develops a nonlinear acceleration method,
nonlinear Truncated Generalized Conjugate Residual method (nlTGCR), that takes
advantage of symmetry. This work is motivated by the observation that the
Hessian of a nonlinear function, or the Jacobian of a gradient of a mapping,
$f$ is symmetric and therefore more effective, conjugate gradient-like schemes
can be exploited.

We demonstrate that nonlinear acceleration methods can benefit from the
symmetry property of the Hessian. In particular, we study both linear and
nonlinear problems and give a systematic analysis of TGCR and nlTGCR. We show
that TGCR is efficient and optimal for linear problems. By viewing the method
from the angle of an inexact Newton approach, we also show that adding a few
global convergence strategies ensures that nlTGCR can achieve global
convergence guarantees.

We complement our theoretical results with numerical simulations on several
different problems. The experimental results demonstrate advantages of our
methods. To the best of our knowledge, this is still the first work to
investigate and improve the AA dynamics by exploiting symmetry of the Hessian.

Related work. Designing efficient optimization methods has received much
attention. Several recent works [65, 4, 40, 48, 18] consider second order
optimization methods that employ sketching or approximation techniques.
Different from these approaches, our method is a first-order method that
utilizes symmetry of the Hessian instead of constructing it. A variant of
inexact Newton method was proposed in [47] where the least-squares sub-
problems are solved approximately using Minimum Residual method. Similarly, a
new type of quasi Newton symmetric update [54] uses several secant equations
in a least-squares sense. These approaches have the same goal as ours.
However, they are more closely related to a secant or a multi-secant
technique, and as will be argued it does a better job of capturing the
nonlinearity of the problem. [63] proposed a short-term AA algorithm that is
different from ours because it is still based on the parameter sequence
instead of the gradient sequence and does not exploit symmetry of the Hessian.

##  2 Background

###  2.1 Extrapolation, acceleration, and the Anderson Acceleration procedure

Consider a general fixed-point problem and the associated fixed-point
iteration as shown in (1). Denote by $r_{j}=H(x_{j})-x_{j}$ the _residual_
vector at the $j$th iteration. Classical extrapolation methods including RRE,
MPE and the vector $\epsilon$-Algorithm, have been designed to accelerate the
convergence of the original sequence by generating a new and independent
sequence of the form: $t_{j}^{(k)}=\sum_{i=0}^{k}\alpha_{i}x_{j+i}$. An
important characteristic of these classical extrapolation methods is that the
two sequences are not mixed in the sense that no accelerated item
$t_{j}^{(k)}$, is used to produce the iterate $x_{j}$. These _extrapolation_
methods must be distinguished from _acceleration_ methods such as the AA
procedure which aim at generating their own sequences to find a fixed point of
a certain mapping $H$.

AA was originally designed to solve a system of nonlinear equations written in
the form $F(x)=H(x)-x=0$ [1, 61, 41, 30]. Denote $F_{i}=F(x_{i})$. AA starts
with an initial $x_{0}$ and sets $x_{1}=H(x_{0})=x_{0}+\beta F_{0}$, where
$\beta>0$ is a parameter. At step $j>1$ we define
${X}_{j}=[x_{j-m},\ldots,x_{j-1}],$ and $\bar{F}_{j}=[F_{j-m},\ldots,F_{j-1}]$
along with the differences:

$\displaystyle\mathcal{X}_{j}$ $\displaystyle=[\Delta x_{j-m}\;\cdots\;\Delta
x_{j-1}]\in\mathbb{R}^{n\times m},$ (2) $\displaystyle\mathcal{F}_{j}$
$\displaystyle=[\Delta F_{j-m}\;\cdots\;\Delta F_{j-1}]\in\mathbb{R}^{n\times
m}.$

We then define the next AA iterate as follows:

$\displaystyle x_{j+1}$ $\displaystyle=x_{j}+\beta
F_{j}-(\mathcal{X}_{j}+\beta\mathcal{F}_{j})\ \theta^{(j)}\quad\mbox{where:}$
(3) $\displaystyle\theta^{(j)}$
$\displaystyle=\text{argmin}_{\theta\in\mathbb{R}^{m}}\|F_{j}-\mathcal{F}_{j}\
\theta\|_{2}.$ (4)

To define the next iterate in (3) the algorithm uses the term
$F_{j+1}=F(x_{j+1})$ where $x_{j+1}$ is the current accelerated iterate. AA
belongs to the class of _multi-secant methods_. Indeed, the approximation (3)
can be written as:

$\begin{split}x_{j+1}&=x_{j}-[-\beta
I+(\mathcal{X}_{j}+\beta\mathcal{F}_{j})(\mathcal{F}_{j}^{T}\mathcal{F}_{j})^{-1}\mathcal{F}_{j}^{T}]F_{j}\\\
&\equiv x_{j}-G_{j}F_{j}.\end{split}$ (5)

Thus, $G_{j}$ can be seen as an update to the (approximate) inverse Jacobian
$G_{j-m}=-\beta I$

$G_{j}=G_{j-m}+(\mathcal{X}_{j}-G_{j-m}\mathcal{F}_{j})(\mathcal{F}_{j}^{T}\mathcal{F}_{j})^{-1}\mathcal{F}_{j}^{T},$
(6)

and is the minimizer of $\|G_{j}+\beta I\|_{F}$ under the _multi-secant
condition_ of type II 111Type I Broyden conditions involve approximations to
the Jacobian, while type II conditions deal with the inverse Jacobian.

$G_{j}\mathcal{F}_{j}=\mathcal{X}_{j}.$ (7)

This link between AA and Broyden multi-secant type updates was first unraveled
by Eyert [21] and expanded upon in [46].

###  2.2 Inexact and quasi-Newton methods

Given a nonlinear system of equations $F(x)=0$. Inexact Newton methods [15,
10], start with an initial guess $x_{0}$ and compute a sequence of iterates as
follows

Solve $\displaystyle J(x_{j})\delta_{j}$ $\displaystyle\approx-F(x_{j})$ (8)
Set $\displaystyle x_{j+1}$ $\displaystyle=x_{j}+\delta_{j}$ (9)

Here, $J(x_{j})$ is the Jacobian of $F$ at the current iterate $x_{j}$. In (8)
the system is solved inexactly, typically by some iterative method. In quasi-
Newton methods [17, 17, 54], the inverse of the Jacobian is approximated
progressively. Because it is the inverse Jacobian that is approximated, the
method is akin to Broyden’s second (or type-II) update method. This method
replaces Newtons’s iteration:
$x_{j+1}=x_{j}-DF(x_{j})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}F(x_{j})$
with $x_{j+1}=x_{j}-G_{j}F(x_{j})$ where $G_{j}$ approximates the inverse of
the Jacobian $DF(x_{j})$ at $x_{j}$ by the update formula
$G_{j+1}=G_{j}+(\Delta x_{j}-G_{j}\Delta F(x_{j}))v_{j}^{T}$ in which $v_{j}$
is defined in different ways see [46] for details.

##  3 Exploiting symmetry

In the following, we specifically consider the case where the nonlinear
mapping $F$ is the gradient of some objective function
$\phi:\mathbb{R}^{n}\rightarrow\mathbb{R}$ to be minimized, i.e.,

$F(x)=\nabla\phi(x).$

In this situation, the Jacobian of $F$ becomes $\nabla^{2}\phi$ the Hessian of
$\phi$. An obvious observation here is that the symmetry of the Hessian is not
taken into account in the approximate inverse Hessian update formula (6). This
has only been considered in the literature (very) recently (e.g., [6, 55, 7]).
In a 1983 report, [53] showed that the matrix $G_{j}$ obtained by a multi-
secant method that satisfies the secant condition (7) is symmetric iff the
matrix $\mathcal{X}_{j}^{T}\mathcal{F}_{j}$ is symmetric. It is possible to
explicitly force symmetry by employing generalizations of the symmetric
versions of Broyden-type methods. Thus, the authors of [6, 7] developed a
multisecant version of the Powell Symmetric Broyden (PSB) update due to Powell
[45] while the article [55] proposed a symmetric multisecant method based on
the popular Broyden-Fletcher-Goldfarb-Shanno (BFGS) approach as well as the
Davidon-Fletcher-Powell (DFP) update. However, there are a number of issues
with the symmetric versions of multisecant updates, some of which are
discussed in [55].

We observe that when we are close to the limit, the condition
$\mathcal{X}_{j}^{T}\mathcal{F}_{j}=\mathcal{F}_{j}^{T}\mathcal{X}_{j}$ is
nearly satisfied. This is because if $x^{\ast}$ is the limit with
$F(x^{\ast})=0$ we can write

$\displaystyle F(x_{k})-F(x_{k-1})$ $\displaystyle=[F(x_{k})-F(x^{\ast})]$
(10) $\displaystyle\quad-[F(x_{k-1})-F(x^{\ast})]$
$\displaystyle\approx\nabla^{2}\phi(x^{\ast})(x_{k}\ -\ x_{k-1}).$

This translates to
$\mathcal{F}_{j}\approx\nabla^{2}\phi(x^{\ast})\mathcal{X}_{j}$ from which it
follows that
$\mathcal{X}_{j}^{T}\mathcal{F}_{j}\approx\mathcal{X}_{j}^{T}\nabla^{2}\phi(x^{\ast})\mathcal{X}_{j}$
which is a symmetric matrix under mild smoothness conditions on $\phi$.
Therefore, the issue of symmetry can be mitigated if we are able to develop
nonlinear acceleration methods that take advantage of near-symmetry.

###  3.1 The linear case: Truncated GCR (TGCR)

We first consider solving the linear system $Ax=b$ with a general matrix $A$.
The Generalized Conjugate Residual (GCR) algorithm, see, e.g., [19, 51],
solves this linear system by building a sequence of search directions $p_{i}$,
for $i=0,\cdots,j$ at step $j$ so that the vectors $Ap_{i}$ are orthogonal to
each other. With this property it is easy to generate iterates that minimize
the residual at each step, and this leads to GCR, see [51, pp 195-196] for
details.

Next we will make two changes to GCR. First, we will develop a truncated
version in which any given $Ap_{j}$ is orthogonal to the previous $m$
$Ap_{i}$’s only. This is dictated by practical considerations, because keeping
all $Ap_{i}$ vectors may otherwise require too much memory. Second, we will
keep a set of vectors for the $p_{i}$’s and another set for the vectors
$v_{i}\equiv Ap_{i}$, for $i=1,\cdots,j$ at step $j$ in order to avoid
unnecessary additional products of the matrix $A$ with vectors. The Truncated
GCR (TGCR) is summarized in Algorithm 1.

Algorithm 1 TGCR (m)

1:Input: Matrix $A$, RHS $b$, initial $x_{0}$.

2:Set $r_{0}\equiv b-Ax_{0}$; $v=Ar_{0}$;

3:$v_{0}=v/\|v\|$; $p_{0}=r_{0}/\|v\|$;

4:for $j=0,1,2,\cdots,$ Until convergence do

5: $\alpha_{j}=(r_{j},v_{j})$

6: $x_{j+1}=x_{j}+\alpha_{j}p_{j}$

7: $r_{j+1}=r_{j}-\alpha_{j}v_{j}$

8: $p=r_{j+1}$; $v=Ap$;

9: $i_{0}=\max(1,j-m+1)$

10: for $i=i_{0}:j$ do

11: $\beta_{ij}:=(v,Ap_{i})$

12: $p:=p-\beta_{ij}p_{i}$;

13: $v:=v-\beta_{ij}v_{i}$;

14: end for

15: $p_{j+1}:=p/\|v\|$ ; $v_{j+1}:=v/\|v\|$ ;

16:end for

With $m=\infty$ we obtain the non-restarted GCR method, which is equivalent to
the non-restarted (i.e., full) GMRES. However, when $A$ is symmetric, but not
necessarily symmetric positive definite, then TGCR (1) is identical with TGCR
(m) in exact arithmetic. This leads to big savings both in terms of memory and
in computational costs.

######  Theorem 3.1.

When the coefficient matrix $A$ is symmetric, TGCR (m) generates the same
iterates as TGCR (1) for any $m>0$. In addition, when $A$ is positive
definite, the $k$-th residual vector $r_{k}=b-Ax_{k}$ satisfies the following
inequality where $\kappa$ is the spectral condition number of $A$:

$\|r_{k}\|\leq
2\left[\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right]^{k}\|r_{0}\|.$ (11)

###  3.2 The nonlinear case: nlTGCR

Assume now that we want to solve the nonlinear problem $F(x)=0.$ We need to
make three major changes to Algorithm 1. First, any residual is now the
negative of $F(x)$ so Line 2 and Line 7 must be replaced by $r_{0}=-F(x_{0})$
and $r_{j+1}=-F(x_{j+1})$, respectively. In addition, we originally need to
calculate the products $Ar_{0}$ and $Ap$ in Line 2 and Line 8 respectively.
Here $A$ needs to be replaced by the Jacobian $J(x_{j})$ of $F$ at the current
iterate. We also use the notation $P_{j}\equiv[p_{i_{0}},\cdots,p_{j}]$, and
$V_{j}\equiv[J(x_{i_{0}})p_{i_{0}},\cdots,J(x_{j})p_{j}]$. The most important
change is in lines 5-6 where $\alpha_{j}$ of Algorithm 1 needs to be replaced
by a vector $y_{j}$. This is because when we write the linear model used in
the form of an inexact Newton method:

$\displaystyle F(x_{j}+P_{j}y)$ $\displaystyle\approx
F(x_{j})+[J]P_{j}y\quad\mbox{where}$ (12) $\displaystyle[J]P_{j}$
$\displaystyle\equiv[J(x_{i_{0}})p_{i_{0}},\cdots,J(x_{j})p_{j}]=V_{j}.$

The projection method that minimizes the norm
$\|F(x_{j})+[J]P_{j}y\|=\|F(x_{j})+V_{j}y\|$ of the right-hand side determines
$y$ in such a way that

$\displaystyle F(x_{j})+V_{j}y\perp\mbox{Span}\\{V_{j}\\}$
$\displaystyle\rightarrow(V_{j})^{T}[F(x_{j})+V_{j}y]=0$ (13)
$\displaystyle\rightarrow y=V_{j}^{T}r_{j}$

where it is assumed the $v_{i}$’s are fully orthogonal. Note that in the
linear case, it can be shown that $V_{j}^{T}r_{j}$ has only one nonzero
component when one assumes that the vectors $J(x_{i})p_{i}$ are fully
orthogonal, i.e., that $i_{0}=1$ always. The nonlinear version of TGCR (m) is
summarized in Algorithm 2 where the indication ‘Use Frechet’ means that the
vector $v=J(x)u$ is to be computed as $v=(F(x+\epsilon u)-F(x))/\epsilon$ for
some small $\epsilon$.

Algorithm 2 nlTGCR (m)

1:Input: $F(x)$, initial $x_{0}$.

2:Set $r_{0}=-F(x_{0})$.

3:Compute $v=J(x_{0})r_{0}$; (Use Frechet)

4:$v_{0}=v/\|v\|$, $p_{0}=r_{0}/\|v\|$;

5:for $j=0,1,2,\cdots,$ Until convergence do

6: $y_{j}=V_{j}^{T}r_{j}$

7: $x_{j+1}=x_{j}+P_{j}y_{j}$

8: $r_{j+1}=-F(x_{j+1})$

9: Set: $p:=r_{j+1}$;

10: $i_{0}=\max(1,j-m+1)$

11: Compute $v=J(x_{j+1})p$ (Use Frechet)

12: for $i=i_{0}:j$ do

13: $\beta_{ij}:=(v,v_{i})$

14: $p:=p-\beta_{ij}p_{i}$

15: $v:=v-\beta_{ij}v_{i}$

16: end for

17: $p_{j+1}:=p/\|v\|$ ; $v_{j+1}:=v/\|v\|$ ;

18:end for

######  Remark.

nlTGCR (m) requires 2 function evaluations per step: one in Line 8 and the
other in Line 11. In the situation when computing the Jacobian is inexpensive,
then one can compute $Jp$ in Line 11 as a matrix-vector product and this will
reduce the number of function evaluations per step from 2 to 1. The inner loop
in Line 12-16 corresponds to exploiting symmetry of Hessian. At a given step,
nlTGCR (m) attempts to approximate a Newton step: $x_{j+1}=x_{j}+\delta$ where
$\delta$ is an approximate solution to $J(x_{j})\delta+F(x_{j})=0$.

High-Level Clarification. At this point, one might ask the question: why not
just use an inexact Newton method whereby the Jacobian system is solved with
the _linear_ GCR or TGCR method? This is where AA provides an interesting
insight on some weaknesses of Newton-Krylov method. A Newton-Krylov method
generates a Krylov subspace $\mbox{Span}\\{r_{0},Jr_{0},\cdots,J^{k}r_{0}\\}$
at a current iterate – say $K=x_{0}$ – (so $J\equiv J(x_{0})\equiv DF(x_{0})$)
and tries to minimize $F(x_{0}+\delta)$ where $\delta\ \in\ K$, by exploiting
the linear model: $F(x_{0}+\delta)\approx F(x_{0})+J\delta$. If we generate a
basis $V=[v_{1},\cdots v_{k}]$ of $K$ and express $\delta$ as $\delta=Vy$ then
we would need to minimize $\|F(x_{0})+JVy\|$ which is a small least-squares
problem. One usually adds to this a global convergence strategy, e.g., a
linesearch or a trust-region technique to produce the next iterate $x_{1}$.
The problem with this approach is this: _the approximate solution obtained
after $k$ steps of a Krylov subspace approach is based on the Jacobian at the
initial point $x_{0}$._ The intermediate calculation is entirely linear and
based on $J(x_{0})$. It is not exploited in any way to produce intermediate
(nonlinear) iterates which in turn could be used to produce more accurate
information on some local Jacobian. In contrast, a method like AA (or in fact
any of the secant or multisecant methods) will do just this, i.e., it will use
information on the nonlinear mapping near the most recent approximation to
produce the new iterate. This distinction is rather important although if the
problem is nearly linear, then it could make little difference.

In nlTGCR , we try to improve on the Newton-Krylov approach, since our
starting point is TGCR which is generalized to nonlinear problems. We also
take the viewpoint of improving on AA or multisecant methods by not relying on
the approximation $F(x_{j+1})-F(x_{j})\approx J(x_{j+1}-x_{j})$ mentioned
above. This is achieved by adopting the projection viewpoint. Instead of
minimizing $\|F(x_{0})+JPy\|$ as in the inexact Newton mode, we would like to
now minimize $\|F(x_{k})+JPy\|$ where $x_{k}$ is the most recent iterate. This
initial idea leads to difficulty since there is not one but several $J$ at
previous points and a single one of them will not be satisfactory. Thus, we
have a few directions $p_{i}$ just like the differences $\Delta x_{i}$ in
Anderson, _but now each $p_{i}$ will lead to a $J(x_{i})p_{i}$ which - unlike
in AA - is accurately computed and then saved._ This feature is what we
believe makes a difference in the performance of the algorithm, although this
is something that is rather difficult to prove theoretically. We leave it as
future work. Overall, the method described in this paper mixes a number of
ideas coming from different horizons. A further high-level discussion and
detailed complexity analysis are provided in Appendix A.

Next, we analyze two possible versions of nlTGCR in the next two sections. In
what follows we assume that all the $Jp_{i}$’s are computed exactly.

####  3.2.1 Linearized update version

First, we consider a variant of Algorithm 2 which we call the “linearized
update version” – whereby in Line 8 we update the residual by using the linear
model, namely, we replace Line 8 by its linear analogue: 8a:
$r_{j+1}=r_{j}-V_{j}y_{j}$. In addition, the matrix-vector product in Line 11
is performed with $J(x_{0})$ instead of $J(x_{j+1})$. When $F$ is linear, it
turns out that $y_{j}$ has only one nonzero component, namely the last one and
this will yield the standard truncated GCR algorithm. Assume that we perform
$k$ steps of the algorithm to produce $x_{k}$, i.e., that Line 5 is replaced
by ‘‘for $j=0,1,2,\cdots,k$ do’’. Then the algorithm is exactly equivalent to
an _inexact Newton method_ in which GMRES (or GCR) is invoked to solve the
Jacobian linear system [10]. Indeed, in this situation Lines 4-15 of Algorithm
1 and Lines 5-17 of Algorithm 2 are identical. In other words, in Lines 5-17,
Algorithm 2 performs $k$ steps of the GCR algorithm for approximately solving
the linear systems $J(x_{0})\delta=-F(x_{0})$. Note that while the update is
written in progressive form as $x_{j+1}=x_{j}+\alpha_{j}p_{j}$, _the right-
hand side does not change during the algorithm and it is equal to
$r_{0}=-F(x_{0})$. _ In effect $x_{k}$ is updated from $x_{0}$ by adding a
vector from the span of $P_{k}$. See the related global convergence result
shown in Theorem B.7 in the Appendix, for a version of this algorithm that
includes a line-search. A weakness of this linear update version is that the
Jacobian is not evaluated at the most recent update but at $x_{0}$, which in
practice is the iterate at each restart.

####  3.2.2 Non-linear update version with residual check

Next we consider the ‘nonlinear-update version’ as described in Algorithm 2.
This version explicitly enforces the linear optimality condition of GCR, as
represented by the Equation (13). In this section, we will analyze the
convergence of nlTGCR through the function $\phi(x)=\frac{1}{2}\|F(x)\|^{2}$.

In order to prove the global convergence of nlTGCR, we need to make a small
modification to Algorithm 2 because as implemented in Algorithm 2 $P_{j}$ is
fixed and the solution obtained at this step may not necessarily satisfy the
following _residual check_ condition which is often used in inexact Newton
methods [15, 11, 20] to guarantee the global convergence:

$\|F(x_{j})+[J]P_{j}y\|\leq\eta\|F(x_{j})\|,$ (14)

where $\eta<1$ is a parameter.

The residual norm on the left-hand side of (14) is readily available at no
additional cost and this can help devise globally converging strategies, by
monitoring to what extent (14) is satisfied. If (14) is not satisfied, we can
either use a line-search technique 4 or restart the process and take the next
iterate as the output of the fixed point iteration mapping $H$. When the
residual check condition is implemented after Line 8 in Algorithm 2, we can
prove the global convergence of nlTGCR in the next theorem. Similar global
strategies have also been proposed in [68, 49, 23, 64, 59, 43].

######  Theorem 3.2 (Global convergence of nlTGCR with residual check).

Assume $\phi$ is twice differentiable and $F(x)$ is L-lipschitz. If the
residual check is satisfied
$\|J(x_{n})P_{n}y_{n}+F(x_{n})\|\leq\eta_{n}\|F(x_{n})\|$ where
$0\leq\eta_{n}\leq\eta<1$ and $J(x_{n})$ is non-singular and its norm is
bounded from above for all n, then $P_{n}y_{n}$ produced in line 7 of
Algorithm 2 is a descent direction and the iterates $x_{n}$ produced by
Algorithm 2will converge to the minimizer $x^{*}$:

$\displaystyle\lim_{n\rightarrow\infty}\phi(x_{n})=\phi(x^{*})=0.$

In the next theorem, we prove that nlTGCR can achieve superlinear and
quadratic convergence under mild conditions.

######  Theorem 3.3 (Superlinear and quadratic convergence of nlTGCR ).

With the same setting as Theorem 3.2. Assume $\nabla^{2}\phi$ is L-lipschitz.
Consider a sequence generated by Algorithm 2 such that residual check is
satisfied $\|J(x_{n})P_{n}y_{n}+F(x_{n})\|\leq\eta_{n}\|F(x_{n})\|$ where
$0\leq\eta_{n}\leq\eta<1$. Moreover, if the following conditions hold

$\displaystyle\phi(x_{n}+P_{n}y_{n})$
$\displaystyle\leq\phi(x_{n})+\alpha\nabla\phi(x_{n})^{T}P_{n}y_{n}$
$\displaystyle\phi(x_{n}+P_{n}y_{n})$
$\displaystyle\geq\phi(x_{n})+\beta\nabla\phi(x_{n})^{T}P_{n}y_{n}$

for $\alpha<\frac{1}{2}$ and $\beta>\frac{1}{2}$. Then there exists a $N_{s}$
such that $x_{n}\rightarrow x^{*}$ superlinearly for $n\geq N_{s}$ if
$\eta_{n}\rightarrow 0$, as $n\rightarrow\infty$. Moreover, if
$\eta_{n}=O(\|F(x_{n})\|^{2})$, the convergence is quadratic.

If the property of the function is bad (non-expansive/non-convex), it will be
more difficult to satisfy the assumptions of Theorem 3.2 and 3.3. For example,
in Theorem 3.2, the non-singularity and boundedness is required for $J(x)$. If
the function does not satisfy the assumption, say, degenerate at a point, then
the algorithm may not converge to a stationary point.

######  Remark.

This superlinear(quadratic) convergence of nlTGCR does not contradict with the
linear convergence of TGCR shown in 11. 11 is obtained from the equivalence
between TGCR and CG in that short-term recurrence holds for symmetric matrix.
Like CG, TGCR can still have a superlinear convergence rate. In practice, the
second stage of convergence of Krylov Space methods is typically well defined
by the theoretical convergence bound with $\sqrt{\kappa(A)}$ but may be super-
linear, depending on a distribution of the spectrum of the matrix $A$ and the
spectral distribution of the error.

Finally, we analyze the convergence of nlTGCR when the gradient $F$ is
subsampled. In the analysis, we make the following five assumptions.

##### Assumptions for stochastic setting

$\mathbf{A_{1}:}$ The variance of subsampled gradients is uniformly bounded by
a constant $C$, $\operatorname{tr}(Cov(F(x)))\leq C^{2},\ \forall x$.

$\mathbf{A_{2}:}$ The eigenvalues of the Hessian matrix for any sample
$|\mathcal{H}|=\beta$ is bounded from below and above in Loewner order
$\mu_{\beta}I\preceq J(x,\mathcal{H})\preceq L_{\beta}I$. Further more, we
require there is uniform lower and upper bound for all subsmaples. That is,
there exists $\hat{\mu}$ and $\hat{L}$ such that
$0\leq\hat{\mu}\leq\mu_{\beta}\quad\text{and}\quad
L_{\beta}\leq\hat{L}<\infty,\quad\forall\beta\in\mathbb{N}.$ And the full
Hessian is bounded below and above $\mu I\preceq J(x)\preceq LI,\quad\forall
x.$

$\mathbf{A_{3}:}$ Hessian is M-Lipschitz, that is $\|J(x)-J(y)\|\leq
M\|x-y\|,\quad\forall x,y$

$\mathbf{A_{4}:}$The variance of subsampled Hessian is bounded by a constant
$\sigma$.

$\displaystyle\|\mathbb{E}_{\mathcal{H}}[(J(x;\mathcal{H})-J(x))]\|\leq\sigma,\quad\forall
x$ (15)

$\mathbf{A_{5}:}$ There exists a constant $\gamma$ such that
$\mathbb{E}[\|x_{n}-x^{*}\|^{2}]\leq\gamma(\mathbb{E}[\|x_{n}-x^{*}\|])^{2}.$

######  Theorem 3.4 (Convergence of stochastic version of nlTGCR ).

Assume $|\mathcal{H}_{n}|=\beta\geq\frac{16\sigma^{2}}{\mu},~{}\forall n$,
residual check is satisfied for $\eta_{n}\leq\eta\leq\frac{1}{4L}$ and
assumptions $A1-A5$ hold. The iterates generated by the stochastic version
Algorithm 2 converge to $x^{*}$ if $\|x_{k}-x^{*}\|\leq\frac{\mu}{2M\gamma}$
and

$\displaystyle\mathbb{E}\|x_{n+1}-x^{*}\|\leq\frac{3}{4}\mathbb{E}\|x_{n}-x^{*}\|.$
(16)

###  3.3 Connections with other methods

This section explores the connection between nlTGCR with inexact Newton and
AA. We provide the connection between nlTGCR and quasi-Newton in A.1.

1) The inexact Newton viewpoint. Inexact Newton methods minimize
$\|F(x_{0})+J(x_{0})P_{j}y\|$ over $y$ by using some iterative method and
enforcing a condition like

$\|F(x_{0})+J(x_{0})P_{j}y\|\leq\eta\|F(x_{0})\|$

where $\eta<1$ is a parameter, see, e.g., [15, 10, 11, 20]. In nlTGCR, we are
trying to solve a similar equation

$F(x_{j})+J(x_{j})\delta=0$

by minimizing $\|F(x_{j})+[J]P_{j}y\|$. We can prove the following properties
of nlTGCR.

######  Proposition 1.

As defined in Algorithm 2, $\delta_{j}=x_{j+1}-x_{j}=Py_{j}$ minimizes
$\|F(x_{j})+\delta\|$ over vectors of the form $\delta=V_{j}y$, where $y\ \in\
\mathbb{R}^{n_{j}}$ and $n_{j}=j-i_{0}+1$.

As noted earlier, in the linear case, the vector $y_{j}$ has only one nonzero
component, namely the top one. In the general case, it is often observed that
the other components are not zero but small. Let us then suppose that we
replace the update in Lines 6-7 by the simpler form: $\mu_{j}=v_{j}^{T}r_{j}$,
and $x_{j+1}=x_{j}+\mu_{j}p_{j}$. Then the direction
$\delta_{j}=x_{j+1}-x_{j}$ is a descent direction for
$\frac{1}{2}\|F(x)\|^{2}$.

######  Proposition 2.

Assume that $J(x_{j})$ is nonsingular and that $\mu_{j}\equiv
v_{j}^{T}r_{j}\neq 0$. Then $\delta_{j}=\mu_{j}p_{j}$ is a descent direction
for the function $\frac{1}{2}\|F(x)\|^{2}$ at $x_{j}$.

2) The quasi-Newton viewpoint. It is also possible to view the algorithm from
the alternative angle of a quasi-Newton approach instead of inexact Newton. In
nlTGCR, the approximate inverse Jacobian $G_{j}$ at step $j$ is equal to

$G_{j}=P_{j}V_{j}^{T}.$ (17)

If we apply this to the vector $v_{j}$ we get
$G_{j}v_{j}=P_{j}V_{j}^{T}v_{j}=p_{j}=J(x_{j})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}v_{j}.$
So $G_{j}$ inverts $J(x_{j})$ exactly when applied to $v_{j}$. It therefore
satisfies the _secant_ equation ([46, sec. 2.3])

$G_{j}v_{j}=p_{j}.$ (18)

This is equivalent to the secant condition $G_{j}\Delta f_{j}=\Delta x_{j}$
used in Broyden’s second update method.

In addition, the update $G_{j}$ satisfies the ‘no-change’ condition:

$G_{j}q=0\quad\forall q\perp v_{j}.$ (19)

The usual no-change condition for secant methods is of the form
$(G_{j}-G_{j-m})q=0$ for $q\perp\Delta f_{j}$ which in our case would be
$(G_{j}-G_{j-m})q=0$ for $q\perp\ v_{j}$. One can therefore consider that we
are updating $G_{j-m}\equiv 0$. In this sense, we can prove the optimality of
nlTGCR(m).

######  Theorem 3.5 (Optimality of nlTGCR).

The matrix $G_{j}$ in (24) is the best approximation to the inverse Jacobian
$J(x_{j})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}$ of $F(x)$ at $x_{j}$
among all the matrices $G$ whose range
$\operatorname{Range}(G)=\operatorname{Span}\\{V_{j}\\}$ and satisfies the
multisecant equation Equation $GV_{j}=P_{j}$. That is,

$\displaystyle G_{j}=\underset{\\{G\in\mathbb{R}^{d\times
d}|\operatorname{Range}(G)=\operatorname{Span}\\{V_{j}\\},GV_{j}=P_{j}\\}}{\operatorname{arg\min}}\|GJ(x_{i})-I\|.$
(20)

3) Comparison with Anderson Acceleration. Let us set $\beta=0$ in Anderson
Acceleration. Without loss of generality and in an effort to simplify notation
we also assume that $i_{0}=1$ each time. According to (3–4), the $j$-th
iterate becomes simply $x_{j+1}=x_{j}-\mathcal{F}_{j}\theta_{j}$ where
$\theta_{j}$ is a vector that minimizes $\|F_{j}-\mathcal{F}_{j}\theta\|$. For
nlTGCR , we have $x_{j+1}=x_{j}+P_{j}y_{j}$ where $y_{j}$ minimizes
$\|F_{j}+V_{j}y\|$. So this is identical with Equation (3) when $\beta=0$ in
which $P_{j}\equiv\mathcal{X}_{j}$, and $\mathcal{F}_{j}$ is replaced by
$V_{j}$.

The most important relation for both cases is the multi-secant relation. For
Anderson, with $G_{j-m}=0$, the multi-secant matrix in (6) becomes

$G_{j}=\mathcal{X}_{j}(\mathcal{F}_{j}^{T}\mathcal{F}_{j})^{-1}\mathcal{F}_{j}^{T}$
(21)

which can be easily seen to minimizes $\|G\|_{F}$ for matrices $G$ that
satisfy the multisecant condition $G\mathcal{F}_{j}=\mathcal{X}_{j}$ and the
no-change condition $G_{j}^{T}(G_{j}-G)=0.$ Therefore the two methods differ
mainly in the way in which the sets $\mathcal{F}_{j}/V_{j}$ , and
$\mathcal{X}_{j}/P_{j}$ are defined. Let us use the more general notation
$V_{j},P_{j}$ for the pair of subspaces.

In both cases, a vector $v_{j}$ is related to the corresponding $p_{j}$ by the
fact that $v_{j}\approx J(x_{j})p_{j}$. In the case of nlTGCR this relation is
explicitly enforced by a Frechet differentiation (Line 10)– before we perform
an orthogonalization - which combines this vector with others – without
changing the span of the new $P_{j}$ (and also $V_{j}$).

In the case of AA, we have $v_{j}=\Delta F_{j-1}=F_{j}-F_{j-1}$ and the
relation exploited is that

$\displaystyle f_{j}$ $\displaystyle\approx
F_{j-1}+J(x_{j-1})(x_{j}-x_{j-1})\to\Delta f_{j-1}$ (22) $\displaystyle\approx
J(x_{j-1})\Delta x_{j-1}$

However, the approximation $v_{j}\approx J(x_{j})p_{j}$ in nlTGCR is _more
accurate_ \- because we use an additional function evaluation to explicitly
obtain a more accurate approximation (ideally exact value) for
$J(x_{j})p_{j}$. In contrast when $x_{j}$ and $x_{j-1}$ are not close, then
(22) can be a very rough approximation. This is a key difference between the
two methods.

##  4 Experimental Results

This section compares our proposed algorithms TGCR and nlTGCR to existing
methods in the literature with the help of a few experiments. We first compare
the convergence for linear problems and then for a softmax classification
problem in the case where the gradients are either deterministic or
stochastic. More experiments and experimental details are available in the
Appendix C.

###  4.1 Linear Problems

We first compare the performance on linear equations
$\mathbf{A}\mathbf{x}=\mathbf{b}$ with Conjugate Gradient [29], generalized
minimal residual method (GMRES) [52] and Anderson Acceleration under different
settings.

##### Linear Systems.

The advantages of TGCR for linear systems are two-folds. 1:) Theorem 11 shows
that TGCR ($1$) is already optimal (equivalent to Conjugate Residual) when
$\mathbf{A}$ is symmetric positive definite. A larger table size is
unnecessary while AA and GMRES require more past iterates to converge fast. It
can be observed from Figure 1(a) and 1(b) that TGCR ($1$) requires much less
memory and computation overhead to converge compared to GMRES and AA. It also
has the same convergence behavior and similar running time as CG. 2:) It is
easy to show that TGCR can converge for indefinite systems while CG fails.
Figure 1(c) verifies our point. This can be helpful when it is not known in
advance if the system is indefinite. The numerical results shown in Figure 1
demonstrate the power of TGCR as a variant of Krylov subspace methods. Figure
1 clearly verifies the correctness of Theorem 11 that TGCR (1) is identical
with TGCR $(m)$ in exact arithmetic, which leads to big savings both in terms
of memory and in computational costs. We include more experimental results in
the Appendix C.2.

(a)  $\mathbf{A}$ is SPD

(b) Time comparison (1(a))

(c) $\mathbf{A}$ is symmetric indefinite

Figure 1: Linear Systems
$\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{A}\in\mathbb{R}^{1000\times 1000}$:
1(a): Comparison in terms of iteration, TGCR [$m$, mv] means table size$=m$
and moving window (no restart). 1(b): Comparison in terms of time for problem
in 1(a). 1(c): Indefinite System. It is well known that CG fails for
indefinite systems. The gap between full GMRES and TGCR is due to the
numerical issue. It can be concluded that TGCR is ideal for solving linear
systems because of its nice convergence property (compared to CG and AA) as
well as the memory-efficient design (compared to GMRES).

(a) Function Value

(b) Time Comparison

(c) Test Accuracy

(d) Stochastic Gradients

Figure 2:  Softmax Classification (MNIST Dataset): 2(a): Function Value vs.
Iterations; 2(b): Function Value vs. Time; 2(c): Test Accuracy vs. Iterations.
2(d): Function Value vs. Iterations Stochastic gradients are calculated using
a batch size of 500.

###  4.2 Nonlinear Problems: Softmax Classification

Next, we consider a softmax multi-class classification problem shown in (79)
without regularization.

$f=-\frac{1}{s}\sum_{i=1}^{s}\log\left(\frac{e^{w_{y_{j}}^{T}x^{(i)}}}{\sum_{j=1}^{k}e^{w_{j}^{T}x^{(i)}}}\right),$
(23)

where $s$ is the total number of sample, $k$ is the total number of classes,
$x^{(i)}$ is vector of all features of sample $i$, $w_{j}$ is the weights for
the $j^{th}$ class, and $y_{j}$ is the correct class for the $i^{th}$ sample.
We compare nlTGCR with Gradient Descent (GD), Nonlinear Conjugate Gradient
(NCG) [13], L-BFGS [38] and Anderson Acceleration using the MNIST dataset [16]
and report results in Figure 2. Figure 2(a) and 2(b) plot the objective value
vs. iteration number and wall-clock time respectively. It can be seen that
nlTGCR converges significantly faster than baselines even without a line-
search strategy. In addition, for this convex and symmetric problem, it is not
surprising to observe that nlTGCR$(1)$ exhibits a similar convergence rate
with nlTGCR$(m)$, which saves even more memory and computation time. Figure
2(c) shows that nlTGCR greatly outperforms baselines by achieving high test
accuracy in the very early stage of training. Figure 2(d) shows the
effectiveness of our method in the stochastic setting. ‘S-’ stands for a
stochastic version. We use a step size of 0.2 for SGD and a batch size ($B$)
of 500 for all stochastic algorithms. It can be observed that nlTGCR(1) with a
small batch size is comparable with the full batch GD with line-search, which
confirms that TGCR takes advantage of symmetry in a very effective way even in
the case of stochastic gradients.

###  4.3 Deep learning applications

We then evaluate nlTCGR on several widely-used deep learning applications
using different frameworks. We run experiments on image classification using
CNN [42] and ResNet [27], time series forecasting using LSTM [31], and node
classification using GCN [35]. Due to space limitation, we provide full
results in Appendix C.5. It shows that nlTGCR(1) outperforms baselines (SGD,
Nesterov, and Adam) for the above DL experiments, highlighting its
effectiveness in large-scale and stochastic non-convex optimization.

##  5 Conclusion

This paper describes an efficient nonlinear acceleration method that takes
advantage of the symmetry of the Hessian. We studied the convergence
properties of the proposed method and established a few connections with
existing methods. The numerical results suggest that nlTGCR can be a
competitive iterative algorithm from both theoretical and practical
perspectives. We plan to conduct a more detailed theoretical and experimental
investigation of the method for a nonconvex stochastic setting.

Social Impact. This work does not present any foreseeable societal
consequence.

## References

  * [1] D. G. Anderson.  Iterative procedures for non-linear integral equations.  Assoc. Comput. Mach., 12(547):547–560, 1965.
  * [2] W. Azizian, I. Mitliagkas, S. Lacoste-Julien, and G. Gidel.  A tight and unified analysis of gradient-based methods for a whole spectrum of games, 2019.
  * [3] C. Blair.  Problem complexity and method efficiency in optimization (a. s. nemirovsky and d. b. yudin).  SIAM Review, 27(2):264–265, 1985.
  * [4] R. Bollapragada, R. Byrd, and J. Nocedal.  Exact and inexact subsampled newton methods for optimization, 2016.
  * [5] R. Bollapragada, R. H. Byrd, and J. Nocedal.  Exact and inexact subsampled newton methods for optimization.  IMA Journal of Numerical Analysis, 39(2):545–578, 2019.
  * [6] N. Boutet, R. Haelterman, and J. Degroote.  Secant update version of quasi-newton psb with weighted multisecant equations.  Computational Optimization and Applications, 75(2):441–466, 2020\.
  * [7] N. Boutet, R. Haelterman, and J. Degroote.  Secant update generalized version of psb: a new approach.  Computational Optimization and Applications, 78(3):953–982, 2021\.
  * [8] C. Brezinski and M. Redivo-Zaglia.  The simplified topological $\epsilon$-algorithms for accelerating sequences in a vector space.  SIAM Journal on Scientific Computing, 36(5):A2227–A2247, 2014.
  * [9] C. Brezinski, M. Redivo-Zaglia, and Y. Saad.  Shanks sequence transformations and anderson acceleration.  SIAM Review, 60(3):646–669, 2018.
  * [10] P. N. Brown and Y. Saad.  Hybrid Krylov methods for nonlinear systems of equations.  SIAM J. Sci. Stat. Comp., 11:450–481, 1990.
  * [11] P. N. Brown and Y. Saad.  Convergence theory of nonlinear Newton-Krylov algorithms.  SIAM Journal on Optimization, 4:297–330, 1994.
  * [12] S. Cabay and L. W. Jackson.  A polynomial extrapolation method for finding limits and antilimits of vector sequences.  SIAM Journal on Numerical Analysis, 13(5):734–752, 1976.
  * [13] Y. H. Dai and Y. Yuan.  A nonlinear conjugate gradient method with a strong global convergence property.  SIAM Journal on Optimization, 10(1):177–182, 1999.
  * [14] A. d’Aspremont, D. Scieur, and A. Taylor.  Acceleration methods.  Foundations and Trends® in Optimization, 5(1-2):1–245, 2021.
  * [15] R. S. Dembo, S. C. Eisenstat, and T. Steihaug.  Inexact Newton methods.  SIAM J. Numer. Anal., 18(2):400–408, 1982.
  * [16] L. Deng.  The mnist database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine, 29(6):141–142, 2012.
  * [17] J. E. Dennis and J. J. Moré.  Quasi-newton methods, motivation and theory.  SIAM Rev., 19:46–89, 1977.
  * [18] M. Derezinski, J. Lacotte, M. Pilanci, and M. W. Mahoney.  Newton-LESS: Sparsification without trade-offs for the sketched newton update.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [19] S. C. Eisenstat, H. C. Elman, and M. H. Schultz.  Variational iterative methods for nonsymmetric systems of linear equations.  SIAM Journal on Numerical Analysis, 20:345–357, 1983.
  * [20] S. C. Eisenstat and H. F. Walker.  Globally convergent inexact newton methods.  SIAM Journal on Optimization, 4:393–422, 1994.
  * [21] V. Eyert.  A comparative study on methods for convergence acceleration of iterative vector sequences.  J. Computational Phys., 124:271–285, 1996.
  * [22] M. Geist and B. Scherrer.  Anderson acceleration for reinforcement learning, 2018.
  * [23] D. Goldfarb, Y. Ren, and A. Bahamou.  Practical quasi-newton methods for training deep neural networks.  In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, Red Hook, NY, USA, 2020. Curran Associates Inc.
  * [24] W. W. Hager and H. Zhang.  A survey of nonlinear conjugate gradient methods.  2005\.
  * [25] H. He, S. Zhao, Y. Xi, J. Ho, and Y. Saad.  GDA-AM: ON THE EFFECTIVENESS OF SOLVING MIN-IMAX OPTIMIZATION VIA ANDERSON MIXING.  In International Conference on Learning Representations, 2022.
  * [26] H. He, S. Zhao, Y. Xi, J. C. Ho, and Y. Saad.  Solve minimax optimization by anderson acceleration, 2021.
  * [27] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  arXiv preprint arXiv:1512.03385, 2015.
  * [28] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition, 2015.
  * [29] M. R. Hestenes and E. Stiefel.  Methods of conjugate gradients for solving linear systems.  Journal of research of the National Bureau of Standards, 49:409–435, 1952.
  * [30] N. J. Higham and N. Strabi.  Anderson acceleration of the alternating projections method for computing the nearest correlation matrix.  to appear.
  * [31] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural Computation, 1997.
  * [32] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural computation, 9:1735–80, 12 1997.
  * [33] K. Jbilou and H. Sadok.  Lu implementation of the modified minimal polynomial extrapolation method for solving linear and nonlinear systems.  IMA Journal of Numerical Analysis, 19(4):549–561, 1999.
  * [34] C. T. Kelley.  Newton’s method in mixed precision.  SIAM Review, 64(1):191–211, 2022.
  * [35] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.  In International Conference on Learning Representations (ICLR), 2017\.
  * [36] T. N. Kipf and M. Welling.  Semi-Supervised Classification with Graph Convolutional Networks.  arXiv:1609.02907 [cs, stat], Feb. 2017.  arXiv: 1609.02907.
  * [37] A. Krizhevsky.  Learning multiple layers of features from tiny images.  2009\.
  * [38] D. C. Liu and J. Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical Programming, 45:503–528, 1989.
  * [39] V. V. Mai and M. Johansson.  Nonlinear acceleration of constrained optimization algorithms.  In ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4903–4907, 2019.
  * [40] J. Martens.  Deep learning via hessian-free optimization.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, page 735–742, Madison, WI, USA, 2010. Omnipress.
  * [41] P. Ni.  Anderson Acceleration of Fixed-point Iteration with Applications to Electronic Structure Computations.  PhD thesis, Worcester Polytechnic Institute, Worcester, Massachusetts, USA, 2009.
  * [42] K. O’Shea and R. Nash.  An introduction to convolutional neural networks, 2015.
  * [43] W. Ouyang, Y. Liu, and A. Milzarek.  Descent properties of an anderson accelerated gradient method with restarting.  06 2022.
  * [44] M. L. Pasini, J. Yin, V. Reshniak, and M. K. Stoyanov.  Anderson acceleration for distributed training of deep learning models.  In SoutheastCon 2022, pages 289–295, 2022.
  * [45] M. POWELL.  A new algorithm for unconstrained optimization.  In J. Rosen, O. Mangasarian, and K. Ritter, editors, Nonlinear Programming, pages 31–65. Academic Press, 1970.
  * [46] H. ren Fang and Y. Saad.  Two classes of multisecant methods for nonlinear acceleration.  Numerical Linear Algebra with Applications, 16(3):197–221, 2009\.
  * [47] F. Roosta, Y. Liu, P. Xu, and M. W. Mahoney.  Newton-mr: Inexact newton method with minimum residual sub-problem solver, 2018.
  * [48] F. Roosta-Khorasani and M. W. Mahoney.  Sub-sampled newton methods i: Globally convergent algorithms, 2016.
  * [49] C. W. Royer, M. O’Neill, and S. J. Wright.  A newton-cg algorithm with complexity guarantees for smooth unconstrained optimization.  Math. Program., 180(1–2):451–488, mar 2020.
  * [50] Y. Saad.  Numerical Methods for Large Eigenvalue Problems.  Halstead Press, New York, 1992.
  * [51] Y. Saad.  Iterative Methods for Sparse Linear Systems, 2nd edition.  SIAM, Philadelpha, PA, 2003.
  * [52] Y. Saad and M. H. Schultz.  GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7:856–869, 1986.
  * [53] R. B. Schnabel.  Quasi-newton methods using multiple secant equations.  Technical Report CU-CS-247-83, Department of Computer Science, University of Colorado at Boulder, Boulder, CO, 1983.
  * [54] D. Scieur, L. Liu, T. Pumir, and N. Boumal.  Generalization of quasi-newton methods: Application to robust symmetric multisecant updates, 2020.
  * [55] D. Scieur, L. Liu, T. Pumir, and N. Boumal.  Generalization of quasi-newton methods: Application to robust symmetric multisecant updates.  In A. Banerjee and K. Fukumizu, editors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 550–558. PMLR, 13–15 Apr 2021.
  * [56] D. Scieur, E. Oyallon, A. d’Aspremont, and F. Bach.  Online regularized nonlinear acceleration, 2018.
  * [57] W. Shi, S. Song, H. Wu, Y. Hsu, C. Wu, and G. Huang.  Regularized anderson acceleration for off-policy deep reinforcement learning.  In NeurIPS, 2019.
  * [58] D. A. Smith, W. F. Ford, and A. Sidi.  Extrapolation methods for vector sequences.  SIAM Review, 29(2):199–233, 1987.
  * [59] H. D. Sterck and Y. He.  On the asymptotic linear convergence speed of anderson acceleration, nesterov acceleration, and nonlinear gmres.  SIAM Journal on Scientific Computing, 43(5):S21–S46, 2021.
  * [60] K. Sun, Y. Wang, Y. Liu, Y. Zhao, B. Pan, S. Jui, B. Jiang, and L. Kong.  Damped anderson mixing for deep reinforcement learning: Acceleration, convergence, and stabilization.  In NeurIPS, 2021.
  * [61] H. F. Walker and P. Ni.  Anderson acceleration for fixed-point iterations.  SIAM J. Numer. Anal., 49(4):1715–1735, 2011.
  * [62] F. Wei, C. Bao, and Y. Liu.  Stochastic anderson mixing for nonconvex stochastic optimization, 2021\.
  * [63] F. Wei, C. Bao, and Y. Liu.  A class of short-term recurrence anderson mixing methods and their applications.  In International Conference on Learning Representations, 2022.
  * [64] Y. Xie, R. H. Byrd, and J. Nocedal.  Analysis of the bfgs method with errors.  SIAM Journal on Optimization, 30(1):182–209, 2020.
  * [65] Z. Yao, A. Gholami, S. Shen, M. Mustafa, K. Keutzer, and M. W. Mahoney.  Adahessian: An adaptive second order optimizer for machine learning.  2020\.
  * [66] G. Zhang and Y. Yu.  Convergence of gradient methods on bilinear zero-sum games.  In ICLR, 2020.
  * [67] J. Zhang, B. O’Donoghue, and S. Boyd.  Globally convergent type-i anderson acceleration for non-smooth fixed-point iterations, 2018.
  * [68] J. Zhang, B. O’Donoghue, and S. Boyd.  Globally convergent type-i anderson acceleration for nonsmooth fixed-point iterations.  SIAM Journal on Optimization, 30(4):3170–3197, 2020.

##  Appendix A Additional Discussion

###  A.1 High-Level Clarification

The method described in this paper mixes a number of ideas coming from
different horizons. Some of the high-level discussion provided below is
expanding further in later sections.

Linear case: TGCR. Our initial idea was motivated by considering the linear
case, in an attempt to exploit Conjugate-gradient like methods for solving a
linear system $Ax=b$. When $A$ is symmetric, it is known that it is possible
to minimize the objective function $f(x)=\|b-Ax\|_{2}^{2}$ on the $k$-th
Krylov subspace $\mbox{Span}\\{r_{0},Ar_{0},\cdots,A^{k-1}r_{0}\\}$ by a nice
algorithm that uses a short-term recurrence. This algorithm, called the
Conjugate Residual algorithm, is quite similar to the Conjugate Gradient but
its residual vectors are conjugate (instead of being orthogonal) and its
search directions are $A^{T}A-$ conjugate (instead of being $A$-conjugate).
Its generalization to the nonsymmetric case, called the Generalized Conjugate
Residual method, is easy to obtain by enforcing these two properties.
Enforcing the $A^{T}A$ conjugacy of the $p_{i}$’s is the same as enforcing the
orthogonality of the vectors $Ap_{i}$ and this is expensive when we have many
vectors. For this reason, practical versions of the algorithm are _truncated_
, i.e., the orthogonalization is enforced against only a few previous
directions. The result is the TGCR(m) algorithm (Algorithm 1) – which has been
known since the 1980s. It is clear that _we expect that when the matrix $A$ is
nearly symmetric_ TGCR(m) will perform nearly as well as the full version GCR
- because when $A$ is symmetric, taking $m=1$ will yield full orthogonality of
the $Ap_{i}$s (Theorem B.2).

Nonlinear case: Newton Krylov. Suppose now that we have to solve the nonlinear
system $F(x)=0$ (in optimization $F$ is just the gradient of the objective
function). At this point, we may ask the question: why not just use an inexact
Newton method whereby the Jacobian system is solved with the _linear_ GCR or
TGCR method? This is where Anderson acceleration provides an interesting
insight on some weaknesses of Newton-Krylov method. A Newton Krylov method
generates a Krylov subspace $\mbox{Span}\\{r_{0},Jr_{0},\cdots,J^{k}r_{0}\\}$
at a current iterate – say $K=x_{0}$ – (so $J\equiv J(x_{0})\equiv DF(x_{0})$)
and tries to minimize $F(x_{0}+\delta)$ where $\delta\ \in\ K$, by exploiting
the linear model: $F(x_{0}+\delta)\approx F(x_{0})+J\delta$. If we generate a
basis $V=[v_{1},\cdots v_{k}]$ of $K$ and express $\delta$ as $\delta=Vy$ then
we would need to minimize $\|F(x_{0})+JVy\|$ which is a small least-squares
problem. One usually adds to this a global convergence strategies, e.g., a
linesearch or a trust-region technique to produce the next iterate $x_{1}$.
The problem with this approach is this: _the approximate solution obtained
after $k$ steps of a Krylov subspace approach is based on the Jacobian at the
initial point $x_{0}$._ The intermediate calculation is entirely linear and
based on $J(x_{0})$. It is not exploited in any way to produce intermediate
(nonlinear) iterates which in turn could be used to produce more accurate
information on some local Jacobian. In contrast, a method like Anderson
acceleration (or in fact any of the secant or multisecant methods) will do
just this, i.e., it will tend to use information on the nonlinear mapping near
the most recent approximation to produce the new iterate. This distinction is
rather important although if the problem is nearly linear, then it could make
little difference.

Nonlinear case: Anderson and nlTGCR. Anderson acceleration can be viewed as a
form of Quasi-Newton method whereby the approximate inverse Jacobian is
updated at each step by using the collection of the previous iterates
$x_{k},x_{k-1},\cdots x_{k-m+1}$ and the corresponding function values
$F_{k},F_{k-1},\cdots F_{k-m+1}$. To be more accurate it uses the differences
$\Delta x_{j}=x_{j+1}-x_{j}$ and the corresponding $\Delta F_{j}$ defined in
the same way. Similarly to Newton-Krylov, it generates an approximation of the
form $x_{k}+Py$ where $P$ is a basis of the subspace spanned by the $\Delta
x_{j}$’s. Notice how the update now is on $x_{k}$ the latest point generated.
The previous iterates are used to essentially provide information on the
nonlinear mapping and its differential. This information is constantly updated
using the most recent iterate. Note that this is informally stated: Anderson
does not formally get an approximation to the Jacobian. It is based implicitly
on exploiting the relation $F(x_{j+1})-F(x_{j})\approx J(x_{j+1}-x_{j})$.
Herein lies a problem that nlTGCR aims at correcting: this relation is only
vaguely verified. For example, _if we take $J$ to be $J(x_{j})$, the Jacobian
at $x_{j}$, the resulting linear model is bound to be extremely inaccurate at
the beginning of the iteration._

In nlTGCR, we try to improve on the Newton-Krylov approach, since our starting
point is TGCR which is generalized to nonlinear problems. We also take the
viewpoint of improving on Anderson Acceleration or multisecant methods by not
relying on the approximation $F(x_{j+1})-F(x_{j})\approx J(x_{j+1}-x_{j})$
mentioned above. This is achieved by adopting the projection viewpoint.
Instead of minimizing $\|F(x_{0})+JPy\|$ as in the inexact Newton mode, we
would like to now minimize $\|F(x_{k})+JPy\|$ where $x_{k}$ is the most recent
iterate. This initial idea leads to a difficulty since there is not one $J$
but several ones at previous points and a single one of them will not be
satisfactory. Thus, we have a few directions $p_{i}$ just like the differences
$\Delta x_{i}$ in Anderson, _but now each $p_{i}$ will lead to a
$J(x_{i})p_{i}$ which - unlike in AA - is accurately computed and then saved._
This feature is what we believe makes a difference in the performance of the
algorithm – although this is something that would be rather difficult to prove
theoretically.

The Quasi-Newton viewpoint.. It is also possible to view the algorithm from
the alternative angle of a Quasi-Newton approach instead of Inexact Newton. In
this viewpoint, the inverse of the Jacobian is approximated progressively.
Because it is the inverse Jacobian that is approximated, the method is akin to
Broyden’s second update method.

In our case, the approximate inverse Jacobian $G_{j}$ at step $j$ is equal to

$G_{j}=P_{j}V_{j}^{T}.$ (24)

If we apply this to the vector $v_{j}$ we get
$G_{j}v_{j}=P_{j}V_{j}^{T}v_{j}=p_{j}=J(x_{j})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}v_{j}.$
So $G_{j}$ inverts $J(x_{j})$ exactly when applied to $v_{j}$. It therefore
satisfies the _secant_ equation ([46, sec. 2.3])

$G_{j}v_{j}=p_{j}.$ (25)

This is the equivalent to the secant condition $G_{j}\Delta f_{j}=\Delta
x_{j}$ used in Broyden’s second update method. Broyden type-II methods replace
Newtons’s iteration:
$x_{j+1}=x_{j}-Df(x_{j})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}f_{j}$
with $x_{j+1}=x_{j}-G_{j}f_{j}$ where $G_{j}$ approximates the inverse of the
Jacobian $Df(x_{j})$ at $x_{j}$ by the update formula $G_{j+1}=G_{j}+(\Delta
x_{j}-G_{j}\Delta f_{j})v_{j}^{T}$ in which $v_{j}$ is defined in different
ways see [46] for details.

In addition, the update $G_{j}$ satisfies the ‘no-change’ condition:

$G_{j}q=0\quad\forall q\perp v_{j}.$ (26)

The usual no-change condition for secant methods is of the form
$(G_{j}-G_{j-m})q=0$ for $q\perp\Delta f_{j}$ which in our case would be
$(G_{j}-G_{j-m})q=0$ for $q\perp\ v_{j}$. One can therefore consider that we
are updating $G_{j-m}\equiv 0$.

It is also possible to find a link between the method proposed herein and the
Anderson acceleration, by unraveling a relation with multi-secant methods.
Note that equation (25) is satisfied for (at most) $m$ previous instances of
$j$, i.e., at step $j$ we have ($i_{0}$ defined in the algorithm)
$G_{j}v_{i}=p_{i}$ for $i=i_{0},\cdots,j.$ In other words we can also write

$G_{j}V_{j}=P_{j}.$ (27)

This is similar to the multi-secant condition $G_{j}{\cal F}_{j}={\cal X}_{j}$
of Equation (7) – see also equation (13) of [46] where ${\cal F}_{j}$ and
${\cal X}_{j}$ are defined in (2). In addition, we clearly also have a multi
secant version of the no-change condition (26) seen above, which becomes:

$G_{j}q=0\quad\forall\quad q\perp\mbox{Span}\\{V_{j}\\}.$ (28)

This is similar to the no-change condition represented by eq. (15) of [46],
which stipulates that $(G_{j}-G_{j-m})q=0$ for all $q$ orthogonal to the span
of the subspace $\mbox{Span}\\{{\cal F}_{j}\\}$ mentioned above, provided we
define $G_{j-m}=0$.

###  A.2 Complexity Analysis

Assume that the iteration number is $k$ and the model parameter size is $d$.
The full memory AA stores all previous iterations, thus the additional memory
is $2kd$. To reduce the memory overhead, the limited-memory (Truncated)
AA$(m)$ maintains the most recent $m$ iterations while discarding the older
historical information. In comparison, TGCR and NLTGCR only requires the most
recent iterate to achieve optimal performance, thus the additional memory is
$2d$. The reduced number of past iterates also saves the orthogonalization
costs from TGCR and NLTGCR compared to AA(m). In TGCR and NLTGCR, only one
orthogonalization is needed to performed which costs $O(kd)$ while AA(m)
requires $O(k^{2}d)$.

For TGCR(m), $(2d-1)$ flops are performed in Line 5, $4d$ flops are performed
in Lines 6-7 and $m(6d-1)$ flops are performed in the for loop and $2d$ flops
are performed in Line 15. If TGCR(m) performs k iterations, the computational
complexity is $((6m+8)d-1-m)k$. Thus, TGCR costs $O(mdk)$. For symmetric
problems, $m=1$ is guaranteed to generate the same iterates as $m>1$ and TGCR
costs $O(dk)$.

Then we analyze the complexity of nlTGCR(m). $m(2d-1)$ flops are performed in
Line 6, $2md$ flops are performed in Line 7, two evaluations of $F$ are
performed in Lines 8 and 11. The for loop costs $m(6d-1)$ flops and $2d$ flops
are performed in Line 15. When $k$ iterations are performed, nlTGCR costs
$O(mdk)$ plus the costs of $2k$ function evaluations of $F$. When $m=1$ is
used in nonlinear problems, nlTGCR costs $O(dk)$ plus the costs of $2k$
function evaluations of $F$.

###  A.3 The Frechet derivative

In vector analysis, derivatives provide local linear approximations. Frechet
differentiation can be used to calculate directional derivatives of gradients.
We use Frechet Differentiation to compute the directional derivative of a
gradient mapping $f$ at $x$ in direction $h$, which is $v=J(x_{j+1})p$ in
algorithm 2. We define Frechet derivative as follows,

######  Definition 1.

Let $(S,\|\cdot\|)$ and $(T,\|\cdot\|)$ be two normed spaces and let $X$ be an
open set in $(S,\|\cdot\|)$.

A function $f:X\longrightarrow T$ is Fréchet differentiable at $x_{0}$, where
$x_{0}\in X$, if there exists a linear operator
$\left(D_{X}f\right)\left(x_{0}\right):X\longrightarrow T$ such that

$\lim_{h\rightarrow
0}\frac{\left\|f\left(x_{0}+h\right)-f\left(x_{0}\right)-\left(D_{x}f\right)\left(x_{0}\right)(h)\right\|}{\|h\|}=0$

The operator $\left(D_{x}f\right)\left(x_{0}\right):X\longrightarrow T$ is
referred to the Fréchet derivative at $x_{0}$.

##  Appendix B Proofs

###  B.1 Optimality for Linear Problem

We can write the Generalized Conjugate residual formally as follows

Algorithm 3 GCR

1:Input: Matrix $A$, RHS $b$, initial $x_{0}$.

2:Set $p_{0}=r_{0}\equiv b-Ax_{0}$.

3:for $j=0,1,2,\cdots,$ Until convergence do

4: $\alpha_{j}=(r_{j},Ap_{j})/(Ap_{j},Ap_{j})$

5: $x_{j+1}=x_{j}+\alpha_{j}p_{j}$

6: $r_{j+1}=r_{j}-\alpha_{j}Ap_{j}$

7: $p_{j+1}=r_{j+1}-\sum_{i=1}^{j}\beta_{ij}p_{i}$ where
$\beta_{ij}:=(Ar_{j+1},Ap_{i})/(Ap_{i},Ap_{i})$

8:end for

######  Theorem B.1 (Lemma 6.21 in [51].).

If $\\{p_{0},\dots,p_{n-1}\\}$ is the basis of the Krylov space
$\mathcal{K}_{n}(A,r_{0})$ which are also $A^{T}A$ orthogonal . Then

$x_{n}=x_{0}+\sum_{i=0}^{n-1}\frac{\langle r_{0},Ap_{i}\rangle}{\langle
Ap_{i},Ap_{i}\rangle}p_{i}$

minimizes the residual among all the iterates with form
$x_{0}+\mathcal{K}_{n}(A,r_{0})$. Further more, we have

$\displaystyle x_{n}=x_{n-1}+\frac{\langle r_{n-1},Ap_{n-1}\rangle}{\langle
Ap_{n-1},Ap_{n-1}\rangle}p_{n-1}$

###### Proof.

We can write $x_{n}=x_{0}+\sum_{i=0}^{n-1}\beta_{i}p_{i}$ and
$r_{n}=r_{0}-\sum_{i=0}^{n-1}\beta_{i}Ap_{i}$. Since $x_{n}$ minimizes the
residual, we know the following Petrov–Galerkin condition must hold

$\displaystyle(r_{n},Ap_{j})=0,\quad j=0,\dots,n-1$

The $A^{T}A$ orthogonality gives us

$\displaystyle\beta_{i}=\frac{\langle r_{0},Ap_{i}\rangle}{\langle
Ap_{i},Ap_{i}\rangle}.$

Similarly, we can write $x_{n}=x_{n-1}+\beta_{n-1}p_{n-1}$ and
$r_{n}=r_{n-1}-\beta_{n-1}Ap_{n-1}$. Agagin, the optimality condition reads

$\displaystyle\langle r_{n},p_{n-1}\rangle=0$

which gives us

$\displaystyle\frac{\langle r_{n-1},Ap_{n-1}\rangle}{\langle
Ap_{n-1},Ap_{n-1}\rangle}$

∎

######  Theorem B.2.

When the coefficient matrix $A$ is symmetric, TGCR(m) generates exactly the
same iterates as TGCR(1) for any $m>0$.

###### Proof.

Lines 8 to 14 in Algorithm 1 computes the new direction $p_{j+1}$ – by ortho-
normalizing the vector $Ar_{j+1}$ against all previous $Ap_{i}$’s. In fact the
loop of lines 9–13, implements a modified Gram-Schmidt procedure, which in
exact arithmetic amounts simply to setting $p_{j+1}$ to

$\beta_{j+1,j}p_{j+1}:=r_{j+1}-\sum_{i=i_{0}}^{j}\beta_{ij}p_{i}\quad\mbox{where}\quad\beta_{ij}=(Ar_{j+1},Ap_{i})\
\mbox{for}\ i_{0}\leq i\leq j.$ (29)

In the above relation, $\beta_{j+1,j}$ is the scaling factor $\|v\|$ used to
normalize $p$ and $v$ in Line 14. Then, $v_{j+1}\equiv Ap_{j+1}$ is computed
accordingly as is reflected in lines 12 and 14. The update relation
$Ap_{j+1}=Ar_{j+1}-\sum_{i=i_{0},j}\beta_{ij}Ap_{i}$ (from Line 12) shows that
$Ap_{j+1}\perp Ap_{i}$ for $i=i_{0},...,j$. In addition, it can easily be
shown that in this case ($m=\infty$) the residual vectors produced by the
algorithm are $A$-conjugate in that $(r_{j+1},Ar_{i})=0$ for $i\leq j$.
Indeed, this requires a simple induction argument exploiting the equality:

$(r_{j+1},Ar_{i})=(r_{j}-\alpha_{j}Ap_{j},Ar_{i})=(r_{j},Ar_{i})-\alpha_{j}(Ap_{j},Ar_{i})$

and relation (29) which shows that $Ar_{i}=\sum\beta_{k,i-1}Ap_{k}$.

When $A$ is symmetric, exploiting the relation
$r_{i+1}=r_{i}-\alpha_{i}Ap_{i}$, we can see that the scalar $\beta_{ij}$ in
Line 11 of Algorithm 1 is

$\beta_{ij}=(Ar_{j+1},Ap_{i})=\frac{1}{\alpha_{i}}(Ar_{j+1},r_{i}-r_{i+1})=\frac{1}{\alpha_{i}}(r_{j+1},Ar_{i}-Ar_{i+1})$

which is equal to zero for $i<j$. Therefore we need to orthogonalize
$Ar_{j+1}$ against vector $Ap_{j}$ only in the loop of lines 9 to 13. This
completes the proof. ∎

######  Theorem B.3.

Let $\widehat{\mathbf{x}}_{t}$ be the approximate solution obtained at the
t-th iteration of TGCR being applied to solve
$\mathbf{A}\mathbf{x}=\mathbf{b}$, and denote the residual as
$\mathbf{r}_{t}=\mathbf{b}-\mathbf{A}\widehat{\mathbf{x}}_{t}$. Then,
$\mathbf{r}_{t}$ is of the form

$\mathbf{r}_{t}=f_{t}(\mathbf{A})\mathbf{r}_{0},$ (30)

where

$\|\mathbf{r}_{t}\|_{2}=\|f_{t}(\mathbf{A})\mathbf{r}_{0}\|_{2}=\min_{f_{t}\in\mathcal{P}_{t}}\|f_{t}(\mathbf{A})\mathbf{r}_{0}\|_{2},$
(31)

where $\mathcal{P}_{p}$ is the family of polynomials with degree p such that
$f_{p}(0)=1,\forall f_{p}\in\mathcal{P}_{p}$, which are usually called
residual polynomials.

######  Theorem B.4 (Convergence of TGCR (Indefinite Case)).

Suppose $\mathbf{A}$ is hermitian, invertible, and indefinite. Divide its
eigenvalues into positive and negative sets $\Lambda_{+}$and $\Lambda_{-}$,
and define

$\kappa_{+}=\frac{\max_{\lambda\in\Lambda_{+}}|\lambda|}{\min_{\lambda\in\Lambda_{+}}|\lambda|},\quad\kappa_{-}=\frac{\max_{\lambda\in\Lambda-}|\lambda|}{\min_{\lambda\in\Lambda_{-}}|\lambda|}$

Then $\mathrm{x}_{m}$, the $m$ th solution estimate of TGCR, satisfies

$\frac{\left\|\mathbf{r}_{m}\right\|_{2}}{\|\mathbf{b}\|_{2}}\leq
2\left(\frac{\sqrt{\kappa_{+}\kappa_{-}}-1}{\sqrt{\kappa_{+}\kappa_{-}}+1}\right)^{\lfloor
m/2\rfloor}$

where $\lfloor m/2\rfloor$ means to round $m/2$ down to the nearest integer.

###### Proof.

When A is hermitian indefinite, an estimate on the min-max approximation

$\frac{\left\|\mathbf{r}_{m}\right\|_{2}}{\left\|\mathbf{b}\right\|_{2}}\leq\min_{p\in\mathcal{P}_{m}}\max_{k}\left|p\left(\lambda_{k}\right)\right|$
(32)

that represents the worst-case TGCR convergence behavior, can be obtained by
replacing the discrete set of the eigenvalues by the union of two intervals
containing all of them and excluding the origin, say $\Lambda_{+}$and
$\Lambda_{-}$. Then the classical bound for the min-max value can be used to
obtain an estimate for the convergence of the residual [50]

$\displaystyle\min_{p\in\mathcal{P}_{m}}\max_{k}\left|p\left(\lambda_{k}\right)\right|$
$\displaystyle\leq\min_{p\in\mathcal{P}_{m}}\max_{z\in\Lambda_{+}\cup\Lambda_{-}}|p(z)|$
$\displaystyle\leq
2\left(\frac{\sqrt{\kappa_{+}\kappa_{-}}-1}{\sqrt{\kappa_{+}\kappa_{-}}+1}\right)^{\lfloor
m/2\rfloor},$

where $[m/2]$ denotes the integer part of $m/2$. ∎

The optimality of TGCR(m) is proved in Theorem B.5.

######  Theorem B.5.

Let $r$ be the residual generated by the basic TGCR (m), the following
relations hold:

  1. 1.

$\mbox{Span}\\{p_{0},\cdots,p_{m-1}\\}=\mbox{Span}\\{r_{0},\cdots,A^{m-1}r_{0}\\}\equiv
K_{m}(r_{0},A)$

  2. 2.

If $m\geq j$, the set of vectors $\\{Ap_{i}\\}_{i=1:j}$ is orthonormal.

  3. 3.

More generally: $(Ap_{i},Ap_{j})=\delta_{ij}$ for $|i-j|\leq m-1$

  4. 4.

If $m\geq j$, then $\|b-Ax_{j}\|=\min\\{\|b-Ax\|\quad|\quad\ x\ \in\
x_{0}+K_{m}(r_{0},A)\\}$

######  Proposition 3.

Assume that $J(x_{j})$ is nonsingular and that $\mu_{j}\equiv
v_{j}^{T}r_{j}\neq 0$. Then $\delta_{j}=\mu_{j}p_{j}$ is a descent direction
for the function $\phi(x)=\frac{1}{2}\|F(x)\|^{2}$ at $x_{j}$.

###### Proof.

It is known [11] that the gradient of $\phi(x)$ at $x$ is
$\nabla\phi(x)=J(x)^{T}F(x)$. In order for $p_{j}$ to be a descent direction
at $x_{j}$ it is sufficient that the inner product of $p_{j}$ and
$\nabla\phi(x_{j})$ is negative. Consider this inner product

$(\nabla\phi(x_{j}),\mu_{j}p_{j})=\mu_{j}(J(x_{j})^{T}F(x_{j}),p_{j})=\mu_{j}(F(x_{j}),J(x_{j})p_{j})=\mu_{j}(-r_{j},v_{j})=-\mu_{j}^{2}<0.$
(33)

which proves the result. ∎

###  B.2 Optimality from Quasi-Newton Viewpoint

######  Theorem B.6 (Optimality of nltgcr(m) from Quasi-Newton Viewpoint).

The matrix $G_{j}$ is the best approximation to the inverse Jacobi
$J(x_{i})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}$ of $F(x)$ at $x_{i}$
among all the matrices $G$ whose range
$\operatorname{Range}(G)=\operatorname{Span}\\{V_{j}\\}$. That is,

$\displaystyle G_{j}=\underset{\\{G\in\mathbb{R}^{d\times
d}|GV_{j}=P_{j}\\}}{\operatorname{arg\min}}\|GJ(x_{i})-I\|.$ (34)

###### Proof.

Assume $G$ is an arbitrary matrix satisfying the multisecant condition
$GV_{j}=P_{j}$. We have $(G_{j}-G)G_{j}^{T}=0$ and
$\operatorname{Range(G)}=V_{j}$. This can be derived as follows

$\displaystyle
0=P_{j}(P_{j}^{T}-P_{j}^{T})=P_{j}V_{j}^{T}(G^{T}-G_{j}^{T})=G_{j}(G^{T}-G_{j}^{T}).$

We also have $(G_{j}-G)V_{j}=0$. Then set $\Delta=G-G_{j}$, we have

$\displaystyle\|GJ(x_{i})-I\|$ $\displaystyle=\|(G_{j}+\Delta)J(x_{i})-I\|$
$\displaystyle=\|G_{j}J(x_{i})-I\|+\|\Delta
J(x_{i})\|+2\operatorname{Trace}((G_{j}J(x_{i})-I)^{T}\Delta J(x_{i}))$
$\displaystyle\geq\|G_{j}J(x_{i})-I\|+2\operatorname{Trace}((G_{j}J(x_{i})-I)^{T}\Delta
J(x_{i})).$

Then we prove $\operatorname{Trace}((G_{j}J(x_{i})-I)^{T}\Delta J(x_{i}))=0$.
In order to prove this, we compute the trace explicitly. We will denote the
natural basis in $\mathbb{R}^{d}$ by $\\{e_{l}\\}_{l=1}^{d}$ and $l-th$ column
of $J(x_{i})$ by $J_{l}$.

$\displaystyle\operatorname{Trace}((G_{j}J(x_{i})-I)^{T}\Delta J(x_{i}))$
$\displaystyle=\sum_{l=1}^{d}e_{l}^{T}((G_{j}J(x_{i})-I)^{T}\Delta
J(x_{i}))e_{l}$
$\displaystyle=\sum_{l=1}^{d}(J_{l}^{T}G_{j}^{T}-e_{l}^{T})\Delta J_{l}$

Recall we have
$\operatorname{Range}(G)=\operatorname{Range}(G_{j})=\operatorname{Span}\\{V_{j}\\}$
and $GV_{j}=G_{j}V_{j}=P_{j}$ , so

$\left\\{\begin{array}[]{cc}\Delta J_{l}=0&J_{l}\in V_{j}\\\
J_{l}^{T}G_{j}^{T}=0,~{}GJ_{l}=0&J_{l}\in V_{j}^{\perp}\end{array}\right\\}$

∎

###  B.3 Convergence Analysis

Firstly, we will show the global convergence of the Algorithm 2 from inexact
Newton perspective. Usually, some global strategies like line search or
residue check are required for inexact Newton method to converge for
$\phi(x)$.

Algorithm 4 Linesearch Algorithm

1:$\beta=\max\\{1,\epsilon^{*}\frac{|\nabla\phi(x_{n})^{T}p_{n}|}{\|p_{n}\|^{2}}\\}$.

2:If $\phi(x_{n}+\beta p_{n})\leq\phi(x_{n})+\alpha\beta\nabla\phi(x_{n})$,
then set $\beta_{n}=\beta$ and exit. Else:

3:Shrink $\beta$ to be $\beta\in[\theta_{\min}\beta,\theta_{\max}\beta]$ where
$0<\theta_{\min}\leq\theta_{\max}<1$ .Go back to Step. 2.

In [11], authors showed with the general line search algorithm 4, inexact
Newton-Krylov method can converge globally under some mild conditions.

######  Theorem B.7 (Global Convergence from Algorithm 2 with linearized
update and line search).

Assume $\phi$ is continuously differentiable and $F(x)$ is L-lipschitz.
Furthre more, the residual check is satisfied
$\|J(x_{n})P_{n}y_{n}+F(x_{n})\|\leq\eta\|\nabla f(x_{n})\|$ where
$0\leq\eta_{n}\leq\eta<1$. If $J(x_{n})$ is nonsingular and its norm is
bounded from above for all n, then $P_{n}y_{n}$ produced in line 7 of
Algorithm 2 is a descent direction and the iterates $x_{n}$ produced by
Algorithm 2 with linearized update and line search in Algorithm 4 will
converge to the minimizer:

$\displaystyle\lim_{n\rightarrow\infty}\phi(x_{n})=0$

###### Proof.

Since Algorithm 2 with linearized update is equivalent to inexact Newton
Krylov method with TGCR as the solver for the Jacobian system
$J(x_{n})p_{n}=-F(x_{n})$, the theorem is just a result of Theorem B.8. ∎

######  Theorem B.8 ([11]).

Assume $\phi$ is continuously differentiable and $F(x)$ is L-Lipschitz and let
$p_{n}$ be such that
$\|F(x_{n})+J(x_{n})p_{n}\|_{2}\leq\eta_{n}\|F(x_{n})\|_{2}$ for each
$\eta_{n}\leq\eta<1$. Further more, let the next iterate be decided by
Algorithm 4 and $J(x_{n})$ is nonsingular and bounded from above for all n.
Then

$\displaystyle\lim_{n\rightarrow\infty}\phi(x_{n})=0.$

The proof of this theorem depends on the following lemma in [11],

######  Lemma B.9 (Lemma 3.8 of [11]).

Assume $\phi$ is differentiable and $\nabla\phi$ is L-lipschitz. Let
$0<\alpha<1$ and $p_{n}$ denote a descent direction. Then the iterates
$x_{n+1}=x_{n}+\beta p_{n}$ in Algorithm 4 will generated in finite
backtracking steps and $\beta_{n}$ satisfies

$\displaystyle\beta_{n}\|p_{n}\|_{2}\geq-\frac{\nabla\phi^{T}p_{n}}{\|p_{n}\|}\min\Big{(}\epsilon^{*},\frac{1-\alpha}{L}\theta_{\min}\Big{)}.$

######  Theorem B.10 (Global convergence of nlTGCR with residual check).

Assume $\phi$ is twice differentiable and $F(x)$ is L-lipschitz. If the
residual check is satisfied
$\|J(x_{n})P_{n}y_{n}+F(x_{n})\|\leq\eta_{n}\|F(x_{n})\|$ where
$0\leq\eta_{n}\leq\eta<1$ and $J(x_{n})$ is non-singular and the norm of its
inverse is bounded from above for all n, then $P_{n}y_{n}$ produced in line 7
of Algorithm 2 is a descent direction and the iterates $x_{n}$ produced by
Algorithm 4 will converge to the minimizer $x^{*}$:

$\displaystyle\lim_{n\rightarrow\infty}\phi(x_{n})=\phi(x^{*})=0.$

###### Proof.

Denote $P_{n}y_{n}$ by $p_{n}$, $J(x_{n})p_{n}+F(x_{n})=r_{n}$. Since
$\nabla\phi(x_{n})=J(x_{n})^{T}F(x_{n})$ and
$\|r_{n}\|\leq\eta_{n}\|F(x_{n})\|$, we have
$\nabla\phi(x_{n})^{T}p_{n}=F(x_{n})^{T}r_{n}-F(x_{n})^{T}F(x_{n})\leq(\eta-1)\nabla\|F\|^{2}=-2(1-\eta)\phi$
which implies $p_{n}$ is a descent direction. To see the second part of the
theorem, we have

$\displaystyle\phi(x_{n}+\beta_{n}\alpha p_{n})$
$\displaystyle\leq\phi(x_{n})+\beta_{n}\alpha\nabla\phi(x_{n})^{T}p_{n}$ (35)
$\displaystyle\leq\phi(x_{n})-2\beta_{n}\alpha(1-\eta)\phi(x_{n})=[1-2\beta_{n}\alpha(1-\eta)]\phi(x_{n}).$
(36)

Denote $\min\Big{(}\epsilon^{*},\frac{1-\alpha}{L}\theta_{\min}\Big{)}$ by $C$
then,

$\displaystyle-\beta_{n}\|p_{n}\|\leq
C\frac{\nabla\phi(x_{n})^{T}p_{n}}{\|p_{n}\|_{2}}.$

Inserting it back to Inequality (27), we have

$\phi(x_{n+1})\leq\Big{(}1+2\alpha(1-\eta)C\frac{\nabla\phi(x_{n})^{T}p_{n}}{\|p_{n}\|_{2}^{2}}\Big{)}\phi(x_{n})$
(37)

Denote $2\alpha(1-\eta)C$ by $\lambda$ and
$\frac{\nabla\phi(x_{n})^{T}p_{n}}{\|p_{n}\|_{2}^{2}}$ by $t_{n}$, then
$\phi(x_{n+1})\leq(1+\lambda t_{n})\phi(x_{n})$. Since $\phi(x_{n})$ is
bounded from below and non-increasing by the inequality. It must converge to a
finite limit $\phi^{*}$. If $\phi^{*}=0$, we ’re done. Otherwise, dividing the
Inequality 37 by $\phi(x_{n})$ on both sides, we have

$\frac{\phi(x_{n+1})}{\phi(x_{n})}\leq(1+\lambda t_{n})\rightarrow
1,\quad\text{as}~{}n\rightarrow\infty.$ (38)

We also know $1+\lambda t_{n}\leq 1$. Therefore, $t_{n}\rightarrow 0$, as
$n\rightarrow\infty$. In the above discussion, we showed
$2(1-\eta)\phi(x_{n})\leq|t_{n}|\|p_{n}\|_{2}^{2}$ which implies
$\|p_{n}\|\rightarrow\infty$. Recall $p_{n}=J(x_{n})^{-1}(r_{n}-F(x_{n}))$, we
must have $\|p_{n}\|$ bounded. This contradicts with the fact
$\|p_{n}\|\rightarrow\infty$. Therefore, $\phi^{*}=0$ ∎

To proceed to the superlinear and quadratic convergence results, we need the
following lemma from [20]

######  Lemma B.11.

Assume F is continuously differentiable, $\\{x_{k}\\}$ is a sequence such that
$F(x_{k})\rightarrow 0$, and for each k,

$\displaystyle\|F(x_{k+1})\|\leq\|F(x_{k})\|\quad\text{and}\quad\|F(x_{k})+J(x_{k})p_{k}\|\leq\eta\|F(x_{k})\|$
(39)

where $p_{k}=x_{k+1}-x_{k}$ and $\eta>0$ is independent of k. If $x_{*}$ is a
limit point of $\\{x_{k}\\}$ such that $J(x_{*})$ is nonsingular, then
$F(x_{*})=0$ and $x_{k}\rightarrow x_{*}$. In this lemma, we don’t require
$\eta<1$.

######  Theorem B.12 (Superlinear and quadratic convergence of nlTGCR ).

With the same setting as Theorem B.10. Assume both $\nabla\phi$ and
$\nabla^{2}\phi$ are L-Lipschitz. Consider a sequence generated by Algorithm 2
such that residual check is satisfied
$\|J(x_{n})P_{n}y_{n}+F(x_{n})\|\leq\eta_{n}\|F(x_{n})\|$ where
$0\leq\eta_{n}\leq\eta<1$. Moreover, if the following conditions hold

$\displaystyle\phi(x_{n}+P_{n}y_{n})$
$\displaystyle\leq\phi(x_{n})+\alpha\nabla\phi(x_{n})^{T}P_{n}y_{n}$ (40)
$\displaystyle\phi(x_{n}+P_{n}y_{n})$
$\displaystyle\geq\phi(x_{n})+\beta\nabla\phi(x_{n})^{T}P_{n}y_{n}$ (41)

for $\alpha<\frac{1}{2}$ and $\beta>\frac{1}{2}$. If $x_{n}\rightarrow x_{*}$
with $J(x_{*})$ nonsingular , then $F(x_{*})=0$. Moreover, there exists
$N_{s}$ such that $x_{n}\rightarrow x^{*}$ superlinearly for $n\geq N_{s}$ if
$\eta_{n}\rightarrow 0$, as $n\rightarrow\infty$. Furthermore, if
$\eta_{n}=O(\|F(x_{n})\|^{2})$, the convergence is quadratic.

###### Proof.

In the proof, we denote $Pny_{n}$ by $p_{n}$ for convenience and utilize the
proof of Theorem 3.15 in [11]. According to assumptions, $x_{n}\rightarrow
x_{*}$ with $J(x_{*})$ nonsingular, then $J(x_{n})$ is nonsingular for
$n>n_{J}$ for some large enough $n_{J}$. Next, if $F(x_{n})=0$ for some $n\geq
n_{J}$, then residual check condition will imply $p_{n}=0$ which means
$x_{m}=x_{n}$ for all $m\geq n$. Then the results hold automatically because
the sequence converges in finite steps. Therefore, we can assume $J(x_{n})$ is
nonsingular and $F(x_{n})$ is nonzero for all $n$.

The residual check condition implies $p_{n}$ is a descent direction according
to Lemma B.9. That is, $\nabla\phi^{\top}p_{n}<0$. Then we can show

$\lim_{n\rightarrow\infty}\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}=0.$

To show this notice that according to 40, the following inequality holds

$\displaystyle\phi_{n}-\phi_{n+1}\geq-\alpha\nabla\phi(x_{n}+p_{n})^{\top}(x_{n+1}-x_{n})=\|p_{n}\|\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}$

Since $\phi_{n}$ is monotone decreasing, thus
$\|p_{n}\|\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\rightarrow 0$ as
$n\rightarrow\infty$. To show
$\lim_{n\rightarrow\infty}\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\rightarrow
0$. We also need to apply 41. Firstly, according to mean value theorem, there
exists a $\lambda\in(0,1)$ such that

$\displaystyle\phi_{n+1}-\phi_{n}=\nabla phi(x_{n}+\lambda
p_{n})^{\top}p_{n}.$ (42)

According to 41,

$\displaystyle\phi_{n+1}-\phi_{n}=\nabla\phi(x_{n}+\lambda
p_{n})^{\top}p_{n}\geq\beta\nabla\phi_{n}^{\top}p_{n}.$ (43)

This yields

$\displaystyle[\nabla\phi(x_{n}+\lambda
p_{n})-\nabla\phi(x_{n})]^{\top}p_{n}\geq(\beta-1)\nabla\phi_{n}^{\top}p_{n}>0.$

According to Cauchy-Schwartz inequality,

$\displaystyle(\beta-1)\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\|p_{n}\|\leq\|p_{n}\|\|\nabla\phi(x_{n}+\lambda
p_{n})-\nabla\phi(x_{n})\|\leq L\lambda\|p_{n}\|^{2}$ (44)

Therefore,

$\displaystyle\|p_{n}\|\geq\frac{(\beta-1)}{L\lambda}\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}>0.$
(45)

which means we can draw the conclusion that
$\|p_{n}\|\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\rightarrow 0$ implies
$\frac{\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\rightarrow 0$. If
$x_{n}\rightarrow x_{*}$ with $J(x_{*})$ nonsingular, then by Lemma B.11, we
know $F(x_{*})=0$. According to the definition,

$\displaystyle p_{n}=-J_{n}^{-1}F_{n}+J_{n}^{-1}(F_{n}+J_{n}p_{n}).$ (46)

This implies

$\displaystyle\|p_{n}\|\leq\|J_{n}^{-1}\|\|F_{n}\|+\|J_{n}^{-1}\|\|F_{n}+J_{n}p_{n}\|\leq(1+\eta)\|J_{n}^{-1}\|\|F_{n}\|.$
(47)

We know $\|F_{n}\|\rightarrow 0$ as $x_{n}\rightarrow x_{*}$. The above
inequality implies that $\|p_{n}\|\rightarrow 0$ as $x_{n}\rightarrow x_{*}$
since $J_{n}$ is nonsingular. Denote the residual $F_{n}+J_{n}p_{n}$ by
$r_{n}$. Then $\|r_{n}\|\leq\eta\|F_{n}\|$ and
$p_{n}=J_{n}^{-1}(r_{n}-F_{n})$. Therefore,

$\displaystyle\nabla\phi_{n}^{\top}=(J_{n}^{\top}F_{n})^{\top}J_{n}^{-1}(r_{n}-F_{n})=F^{\top}r-F^{\top}F.$
(48)

This implies

$\displaystyle\frac{|\nabla\phi_{n}^{\top}p_{n}|}{\|p_{n}\|}=\frac{|F^{\top}r-F^{\top}F|}{\|J_{n}^{-1}(r_{n}-F_{n})\|}\geq\frac{|F^{\top}F|-|F^{\top}r|}{\|J_{n}^{-1}(r_{n}-F_{n})\|}.$
(49)

Since $\|r_{n}\|\leq\|F_{n}\|$ implies $|F^{\top}r|\leq\eta\|F_{n}\|^{2}$, we
have

$\displaystyle|F^{\top}F|-|F^{\top}r|\geq(1-\eta)\|F_{n}\|^{2}.$ (50)

Moreover,

$\displaystyle\|J_{n}^{-1}(r_{n}-F_{n})\|\leq\|J_{n}^{-1}\|\|F_{n}\|+\|J_{n}^{-1}r\|_{2}\leq(1+\eta)\|J_{n}^{-1}\|\|F_{n}\|.$
(51)

Finally, we have

$\displaystyle\frac{|\nabla\phi_{n}^{\top}p_{n}|}{\|p_{n}\|}\geq\frac{(1-\eta)\|F_{n}\|^{2}}{(1+\eta)\|J_{n}^{-1}\|\|F_{n}\|}=\frac{(1-\eta)\|F_{n}\|}{(1+\eta)\|J_{n}^{-1}\|}.$
(52)

Using $\|\nabla\phi_{n}\|=\|J_{n}^{\top}F_{n}\|\leq\|J_{n}\|\|F_{n}\|$, we
have

$\displaystyle\frac{|\nabla\phi_{n}^{\top}p_{n}|}{\|\nabla\phi_{n}\|\|p_{n}\|}\geq\frac{(1-\eta)}{(1+\eta)M_{n}},$
(53)

where $M_{n}=\operatorname{cond}_{2}(J_{n})$. Therefore, we have

$\displaystyle\frac{-\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}\geq\frac{(1-\eta)}{(1+\eta)M_{n}}\|\nabla\phi_{n}\|\geq\frac{(1-\eta)}{(1+\eta)M_{n}}\|J_{n}\|^{-1}\|F_{n}\|,$
(54)

This yields,

$\displaystyle\|F_{n}\|\leq\frac{(1+\eta)M_{n}}{(1-\eta)}\|J_{n}\|\frac{-\nabla\phi_{n}^{\top}p_{n}}{\|p_{n}\|}.$
(55)

Hence,

$\displaystyle\|F_{n}\|\|p_{n}\|\leq\frac{(1+\eta)M_{n}}{(1-\eta)}\|J_{n}\|(-\nabla\phi_{n}^{\top}p_{n})=-a_{n}\nabla\phi_{n}^{\top}p_{n}),$
(56)

where $a_{n}=\frac{(1+\eta)M_{n}}{(1-\eta)}\|J_{n}\|$ . Combining 47, we have

$\displaystyle\|p_{n}\|^{2}\leq\frac{(1+\eta)^{2}M_{n}^{2}}{(1-\eta)^{2}}(-\nabla\phi_{n}^{\top}p_{n})=-b_{n}\nabla\phi_{n}^{\top}p_{n}),$
(57)

where $b_{n}=\frac{(1+\eta)^{2}M_{n}^{2}}{(1-\eta)^{2}}$. Next we show the
convergence of the algorithm with the aid of the second order Taylor
expansion. Notice

$\displaystyle\nabla\phi(x)=J^{\top}F(x)=[\nabla F_{1}(x)\dots\nabla
F_{n}(x)]\begin{bmatrix}F_{1}(x)\\\ \vdots\\\ F_{n}(x)\end{bmatrix}$ (58)

The Hessian can be computed as follows

$\displaystyle\nabla^{2}\phi(x)=J^{T}J+\sum_{i=1}^{n}\nabla^{2}F_{i}(x)F_{i}(x)=J^{T}J+G(x),$
(59)

where $\|G(x)\|=\|\sum_{i=1}^{n}\nabla^{2}F_{i}(x)F_{i}(x)\|\rightarrow 0$ as
$x_{n}\rightarrow x_{*}$ since $F(x_{*})=0$. using second order Taylor
expansion, we have

$\displaystyle\phi_{n+1}-\phi_{n}-\frac{1}{2}\nabla\phi_{n}^{\top}p_{n}=\frac{1}{2}(\nabla\phi_{n}+\nabla^{2}\phi(\bar{x})p_{n})^{\top}p_{n}.$
(60)

where $\bar{x}=\gamma x_{n}+(1-\gamma)x_{n+1}$ for some $\gamma\in(0,1)$. Then
we can have

$\displaystyle|\phi_{n+1}-\phi_{n}-\frac{1}{2}\nabla\phi_{n}^{\top}p_{n}|$
$\displaystyle=\frac{1}{2}|(\nabla\phi_{n}+\nabla^{2}\phi(\bar{x})p_{n})^{\top}p_{n}|$
(61)
$\displaystyle=\frac{1}{2}|(\nabla\phi_{n}+\nabla^{2}\phi_{n}p_{n})^{\top}p_{n}+p_{n}^{\top}(\nabla^{2}\phi(\bar{x})-\nabla^{2}\phi_{n})p_{n}|$
$\displaystyle\leq\frac{1}{2}(\|J_{n}^{\top}(F_{n}+J_{n}p_{n})\|\|p_{n}\|+\|G_{n}\|\|p_{n}\|^{2}+L\|p_{n}\|\|p_{n}\|^{2})$
$\displaystyle\leq(\eta\|J_{n}\|\|F_{n}\|\|p_{n}\|+(\|G_{n}\|+L\|p_{n}\|)\|p_{n}\|^{2})$
$\displaystyle\leq-\frac{1}{2}(a_{n}\eta_{n}\|J_{n}\|+b_{n}(\|G_{n}\|+L\|p_{n}\|))\nabla\phi_{n}^{\top}p_{n}$
$\displaystyle=-\frac{1}{2}\epsilon_{n}\phi_{n}^{\top}p_{n},$

where $\epsilon_{n}=a_{n}\eta_{n}\|J_{n}\|+b_{n}(\|G_{n}\|+L\|p_{n}\|)$.
Therefore,

$\displaystyle\frac{1}{2}(1+\epsilon_{n})\nabla\phi_{n}^{\top}p_{n}\leq\phi_{n+1}-\phi_{n}\leq\frac{1}{2}(1-\epsilon_{n})\nabla\phi_{n}^{\top}p_{n}$
(62)

Notice that $\|J_{n}\|$, $a_{n}$ and $b_{n}$ are all bounded from above and
$\eta_{n}$, $\|S_{n}\|$ and $\|p_{n}\|$ all converges to $0$ as
$x_{n}\rightarrow x_{*}$. Therefore $\epsilon_{n}\rightarrow 0$ as
$x_{n}\rightarrow x_{*}$. And choose lager enough $N$ such that for all $n\geq
N$ the following holds

$\displaystyle\epsilon_{n}\leq\min\\{1-2\alpha,2\beta-1\\}.$ (63)

Then for all $n\geq N$, the Goldsetin-Armijo condition is satisfied,

$\displaystyle\beta\phi_{n}^{\top}p_{n}\leq\phi_{n+1}-\phi_{n}\leq\alpha\phi_{n}^{\top}p_{n}.$
(64)

We then finish the proof following Theorem 3.3 in [15]. It’s easy to see

$\displaystyle
J(x_{*})(x_{k+1}-x_{*})=[I+J_{*}(J_{k}^{-1}-J_{*}^{-1})](r_{k}+[J_{k}-J_{*}](x_{k}-x_{*})-[F_{k}-F_{*}-J_{*}(x_{k}-x_{*})])$
(65)

Taking norm yields

$\displaystyle\|x_{k+1}-x_{*}\|$
$\displaystyle\leq[\|J_{*}^{-1}\|+\|J_{*}\|\|J_{k}^{-1}-J_{*}^{-1}\|][\|r_{k}\|+\|J_{k}-J_{*}\|\|x_{k}-x_{*}\|+$
(66) $\displaystyle\quad\|F_{k}-F_{*}-J_{*}(x_{k}-x_{*})\|]$
$\displaystyle=[\|J_{*}^{-1}\|+\|J_{*}\|\|J_{k}^{-1}-J_{*}^{-1}\|][\|r_{k}\|+\|J_{k}-J_{*}\|\|x_{k}-x_{*}\|$
$\displaystyle\quad+\|F_{k}-F_{*}-J_{*}(x_{k}-x_{*})\|]$
$\displaystyle=[\|J_{*}^{-1}\|+o(1)][o(F_{k})+o(1)\|x_{k}-x_{*}\|+o(\|x_{k}-x_{*}\|)]$

Therefore,

$\displaystyle\|x_{k+1}-x_{*}\|=o(F_{k})+o(1)\|x_{k}-x_{*}\|+o(\|x_{k}-x_{*}\|),\quad\text{~{}}k\rightarrow\infty.$
(67)

where we used the fact that for sufficient small $\|y-x_{*}\|$

$\displaystyle\frac{1}{\alpha}\|y-x_{*}\|\leq\|F(y)\|\leq\alpha\|y-x_{*}\|.$
(68)

which is Lemma 3.1 in [15]. Similarly, to show quadratic convergence, juts
notice

$\displaystyle\|F(y)-F(x_{*})-F(x_{*})(y-x_{*})\|\leq
L^{{}^{\prime}}\|y-x_{*}\|^{2}$ (69)

for some constant $L^{{}^{\prime}}$ and sufficient small $\|y-x_{*}\|$. For
more details, check Lemma 3.2 in [15].

∎

###  B.4 Stochastic nlTGCR

Denote the noisy gradient by $F(x;\xi_{\mathcal{G}})$ and the noisy evaluation
of Hessian along a vector $p$ by $J(x;\xi_{\mathcal{H}})p$. The subsample
exact Newton algorithm is defined in Algorithm 5. At $k$-th iteration, we
uniformly subsample $\mathcal{G}_{k},\mathcal{H}_{k}$ from full sample set to
estimate the noisy gradient and Hessian, so both of them are unbiased.

Algorithm 5 subsmaple Exact Newton

1:for $i=1,\dots,k$ do

2: Estimate $F(x_{i};\mathcal{G}_{i})$ and $J(x_{i};\mathcal{H}_{i})$

3: $x_{i+1}\leftarrow
x_{i}-s_{i}J^{-1}(x_{i};\mathcal{H}_{i})F(x_{i};\mathcal{G}_{i})$

4:end for

Before we start the theoretical analysis, we need to make some assumptions
which are usual in stochastic setting.

##### Assumptions for stochastic setting

  1. $E_{1}$

The eigenvalues of Hessian matrix for any sample $|\mathcal{H}|=\beta$ is
bounded form below and above in Loewner order

$\displaystyle\mu_{\beta}I\preceq J(x,\mathcal{H})\preceq L_{\beta}I.$ (70)

Further more, we require there is uniform lower and upper bound for all
subsmaples. That is, there exists $\hat{\mu}$ and $\hat{L}$ such that

$\displaystyle 0\leq\hat{\mu}\leq\mu_{\beta}\quad\text{and}\quad
L_{\beta}\leq\hat{L}<\infty,\quad\forall\beta\in\mathbb{N}.$ (71)

And the full Hessian is bounded below and above

$\displaystyle\mu I\preceq J(x)\preceq LI,\quad\forall x.$ (72)

  2. $E_{2}$

The variance of subsampled gradients is uniformly bounded by a constant $C$.

$\displaystyle\operatorname{tr}(Cov(F(x)))\leq C^{2},\quad\forall x$ (73)

  3. $E_{3}$

Hessian is M-Lipschitz, that is

$\displaystyle\|J(x)-J(y)\|\leq M\|x-y\|,\quad\forall x,y$ (74)

  4. $E_{4}$

The variance of subsampled Hessianis bounded by a constant $\sigma$.

$\displaystyle\|\mathbb{E}_{\mathcal{H}}[(J(x;\mathcal{H})-J(x))]\|\leq\sigma,\quad\forall
x$ (75)

  5. $E_{5}$

There exists a constant $\gamma$ such that

$\displaystyle\mathbb{E}[\|x_{n}-x^{*}\|^{2}]\leq\gamma(\mathbb{E}[\|x_{n}-x^{*}\|])^{2}.$

Firstly, we recall the few results on subsample Newton method from [5].

######  Theorem B.13 (Theorem 2.2 in [5]).

Assume $x_{n}$ is generated by Algorithm 5 with $|\mathcal{G}_{i}|=\eta^{i}$
for some $\eta>1$, $|\mathcal{H}|=\beta\geq 1$ and
$s_{i}=s=\frac{\mu_{\beta}}{L}$ and Assumptions E1-E2 hold, then

$\displaystyle\mathbb{E}_{k}[\phi(x_{k})-\phi(x^{*})]\leq\alpha\tau^{k},$ (76)

where

$\displaystyle\alpha=\max\Big{\\{}\phi(x_{0})-\phi(x^{*}),\frac{C^{2}L_{\beta}}{\mu\mu_{\beta}}\Big{\\}}\quad\text{and}\quad\tau=\max\Big{\\{}1-\frac{\mu\mu_{\beta}}{2LL_{\beta}},\frac{1}{\eta}\Big{\\}}.$

######  Theorem B.14 (Lemma 2.3 form [5]).

Assume $x_{n}$ is generated by Algorithm 5 with $s_{i}\equiv 1$ and
Assumptions E1-E3 hold. Then

$\displaystyle\mathbb{E}_{k}[\|x_{n+1}-x^{*}\|]\leq\frac{1}{\mu_{|\mathcal{H}_{n}|}}\Big{[}\frac{M}{2}\|x_{n}-x^{*}\|^{2}+\mathbb{E}_{k}[\|(J(x_{n};\xi_{\mathcal{H}_{n}})-J(x_{n}))(x_{n}-x^{*})\|]+\frac{C}{\sqrt{|\mathcal{G}_{n}|}}\Big{]}$

######  Lemma B.15 (Lemma 2.4 from [5]).

Assume the assumption E1 and E4 hold. Then

$\displaystyle\mathbb{E}_{k}[\|(J(x_{n};\xi_{\mathcal{H}_{n}})-J(x_{n}))(x_{n}-x^{*})\|]\leq\frac{\sigma}{\sqrt{\mathcal{H}_{n}}}\|x_{k}-x^{*}\|.$
(77)

######  Theorem B.16 (Convergence of stochastic version of nlTGCR ).

Assume $|\mathcal{H}_{n}|=\beta\geq\frac{16\sigma^{2}}{\mu},~{}\forall n$,
residue check is satisfied for $\eta_{n}\leq\eta\leq\frac{1}{4L}$ and
assumptions E1-E5 hold. The iterates generated by the stochastic version
Algoritrhm 2 converges to $x^{*}$ if
$\|x_{k}-x^{*}\|\leq\frac{\mu}{2M\gamma}$.

$\displaystyle\mathbb{E}\|x_{n+1}-x^{*}\|\leq\frac{3}{4}\mathbb{E}\|x_{n}-x^{*}\|$
(78)

###### Proof.

$\displaystyle\mathbb{E}_{n}[\|x_{n+1}-x^{*}\|]=\mathbb{E}_{n}[\|x_{n}-x^{*}-J(x_{n})^{-1}F(x_{n})\|]+\mathbb{E}_{n}[\|J(x_{n})^{-1}F(x_{n})+P_{n}V_{n}^{T}y_{n}\|]$

The first term can be bounded using the Theorem B.14 and Lemma B.15,

$\displaystyle\mathbb{E}_{n}[\|x_{n}-x^{*}-J(x_{n})^{-1}F(x_{n})\|]$
$\displaystyle\leq\frac{1}{\mu_{|\mathcal{H}_{k}|}}\Big{[}\frac{M}{2}\|x_{n}-x^{*}\|^{2}+$
$\displaystyle\mathbb{E}_{k}[\|(J(x_{n};\xi_{\mathcal{H}_{n}})-J(x_{n}))(x_{n}-x^{*})\|]+\frac{C}{\sqrt{|\mathcal{G}_{n}|}}\Big{]}$
$\displaystyle\leq\frac{1}{\mu}[\frac{M}{2}\|x_{n}-x^{*}\|^{2}+\frac{\sigma}{\sqrt{\mathcal{H}_{n}}}\|x_{n}-x^{*}\|]$
$\displaystyle=\frac{M}{2\mu}\|x_{n}-x^{*}\|^{2}+\frac{\sigma}{\mu\sqrt{\beta}}\|x_{n}-x^{*}\|]$

We can bound the second term through the line search, recall at each iteration
we have

$\displaystyle\mathbb{E}_{n}[\|J(x_{n})^{-1}F(x_{n})+P_{n}V_{n}^{T}y_{n}\|]$
$\displaystyle=\mathbb{E}_{n}[\|F(x_{n})+J(x_{n})P_{n}V_{n}^{T}y_{n}\|]$
$\displaystyle\leq\eta_{n}\mathbb{E}_{k}[\|F(x_{n})\|]\leq\eta_{n}L\|x_{n}-x^{*}\|.$

The last inequality comes from the assumption that eigenvalues of $J(x)$ is
uniformly upper bounded by $L$. Finally, combining the above inequalities
gives us

$\displaystyle\mathbb{E}_{n}[\|x_{n+1}-x^{*}\|]\leq\frac{M}{2\mu}\|x_{n}-x^{*}\|^{2}+\frac{\sigma}{\mu\sqrt{\beta}}\|x_{n}-x^{*}\|+\eta_{n}L\|x_{n}-x^{*}\|$

Taking the total expectation on both sides leads to

$\displaystyle\mathbb{E}\mathbb{E}_{n}[\|x_{n+1}-x^{*}\|]$
$\displaystyle=\mathbb{E}\|x_{n+1}-x^{*}\|\leq\frac{M}{2\mu}\mathbb{E}[\|x_{n}-x^{*}\|^{2}+(\frac{\sigma}{\mu\sqrt{\beta}}+\eta_{n}L)\mathbb{E}[\|x_{n}-x^{*}\|]$
$\displaystyle\leq\frac{M\gamma}{2\mu}\mathbb{E}\|x_{n}-x^{*}\|\mathbb{E}\|x_{n}-x^{*}\|+(\frac{\sigma}{\mu\sqrt{\beta}}+\eta_{n}L)\mathbb{E}[\|x_{n}-x^{*}\|]$

We prove the convergence by induction, notice that

$\displaystyle\mathbb{E}\|x_{1}-x^{*}\|$
$\displaystyle\leq\frac{M\gamma}{2\mu}\mathbb{E}\|x_{0}-x^{*}\|\mathbb{E}\|x_{0}-x^{*}\|+(\frac{\sigma}{\mu\sqrt{\beta}}+\eta_{n}L)\mathbb{E}[\|x_{0}-x^{*}\|]$
$\displaystyle\leq\Big{(}\frac{M\gamma}{2\mu}\mathbb{E}\|x_{0}-x^{*}\|+\frac{\sigma}{\mu\sqrt{\beta}}+\eta
L\Big{)}\mathbb{E}[\|x_{0}-x^{*}\|]$
$\displaystyle\leq\Big{(}\frac{M\gamma}{2\mu}*\frac{\mu}{2M\gamma}+\frac{\sigma}{\mu\sqrt{\frac{16\sigma^{2}}{\mu}}}+\frac{1}{4L}*L\Big{)}\mathbb{E}\|x_{0}-x^{*}\|=\frac{3}{4}\mathbb{E}\|x_{0}-x^{*}\|$

Now assume inequality 78 holds for $n-th$ iteration, we prove it for $n+1$-th
iteration

$\displaystyle\mathbb{E}\mathbb{E}_{n}[\|x_{n+1}-x^{*}\|]$
$\displaystyle\leq\Big{(}\frac{M\gamma}{2\mu}\mathbb{E}\|x_{n}-x^{*}\|+\frac{\sigma}{\mu\sqrt{\beta}}+\eta
L\Big{)}\mathbb{E}[\|x_{n}-x^{*}\|]\leq\frac{3}{4}\mathbb{E}[\|x_{n}-x^{*}\|]$

∎

##  Appendix C Experimental Details and More Experiments

In this section, we first include more experimental details that could not be
placed in the main paper due to the space limitation. We then present more
experimental results of NLTGCR for different settings and difficult problems.

###  C.1 Experimental Details

We provide codes implemented in both Matlab and Python. All experiments were
run on a Dual Socket Intel E5-2683v3 2.00GHz CPU with 64 GB memory and NVIDIA
GeForce RTX 3090.

For linear problems considered in Section 4.1,
$\mathbf{A},\mathbf{b},\mathbf{c}$, and initial points are generated using
normally distributed random number. We use
$\mathbf{A}^{T}\mathbf{A}+\alpha\mathbf{I}$ to generate symmetric matrices.
The step size is set as 1 after rescaling $\mathbf{A}$ to have the unit
2-norm. For solving linear equations, we depict convergence by use of the norm
of residual, which is defined as
$\lVert\mathbf{b}-\mathbf{A}\mathbf{x}\rVert$. For solving bilinear games, we
depict convergence by use of the norm of distance to optima, which is defined
as $\lVert\mathbf{w}^{\ast}-\mathbf{w}_{t}\rVert$. For most baselines, we use
the Matlab official implementation.

The softmax regression problem considered in Section 4.2 is defined as
follows,

$f=-\frac{1}{s}\sum_{i=1}^{s}\log\left(\frac{e^{w_{y_{j}}^{T}x^{(i)}}}{\sum_{j=1}^{k}e^{w_{j}^{T}x^{(i)}}}\right),$
(79)

where $s$ is the total number of sample, $k$ is the total number of classes,
$x^{(i)}$ is vector of all features of sample $i$, $w_{j}$ is the weights for
the $j^{th}$ class, and $y_{j}$ is the correct class for the $i^{th}$ sample.

###  C.2 TGCR(1) for linear system

We first test the robustness of TGCR(1) for solving linear systems by running
with 50 different initials. Figure 3 indicates TGCR converge well regardless
of initialization. We then compare the performance on bilinear games with
Anderson Acceleration as [26] shows AA outperforms existing methods on such
problems.

(a) $\mathbf{A}$ is SPD

(b) $\mathbf{A}$ is symmetric indefinite

Figure 3: Linear Systems
$\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{A}\in\mathbb{R}^{100\times 100}$:
Comparison in terms of iteration over 50 random runs. We can observe that
TGCR(1) match full memory GMRES in the first stage and consistently converge
faster than AA(10).

##### Minimax Optimization.

Next, we test TGCR on the following zero-sum bilinear games:

$\min_{\mathbf{x}\in\mathbb{R}^{n}}\max_{\mathbf{y}\in\mathbb{R}^{n}}f(\mathbf{x},\mathbf{y})=\mathbf{x}^{T}\mathbf{A}\mathbf{y}+\mathbf{b}^{T}\mathbf{x}+\mathbf{c}^{T}\mathbf{y},\quad\mathbf{A}\textrm{
is full rank.}$ (80)

Bilinear games are often regarded as an important but simple class of problems
for theoretically analyzing and understanding algorithms for solving general
minimax problems [66, 2]. Here we consider simultaneous GDA mapping for
minimax bilinear games
$\bigl{(}\begin{smallmatrix}\mathbf{I}&-\eta\mathbf{A}\\\
\eta\mathbf{A}^{T}&\mathbf{I}\end{smallmatrix}\bigr{)}$ [26]. Although this
mapping is skew-symmetric, TGCR can still exploit the short-term recurrence.
It can be observed from Figure 4 that Krylov subspace methods such as TGCR and
AA converge fast for bilinear problem when $A$ is either SPD or random
generated. More importantly, Figure 4 demonstrates that TGCR $(1)$ exhibits a
_superlinear_ convergence rate and converges to optimal significantly faster
than AA$(m)$ in terms of both iteration number and computation time.

(a) $\mathbf{A}$ is SPD

(b) Time Comparison (4(a))

(c) $\mathbf{A}$ is random generated

Figure 4: Bilinear Problems: 4(a): Distance to optimal vs. Iterations; 4(b):
Distance to optimal vs. Time; 4(c): Distance to optimal vs. Iterations. It can
be observed that the short-term property (Theorem 11) holds as long as the
mapping is symmetric or skew-symmetric.

###  C.3 TGCR(1) for nonsymmetric quadratic minimization and linear system

A quadratic form is simply a scalar, quadratic function of a vector with the
form

$f(x)=\frac{1}{2}\mathbf{x}^{T}\mathbf{A}\mathbf{x}-\mathbf{b}^{T}\mathbf{x}+c,$
(81)

where $\mathbf{A}$ is a matrix, $\mathbf{x}$ and $\mathbf{b}$ are vectors, and
$c$ is a scalar constant. When $\mathbf{A}$ is symmetric and positive-
definite, $f(x)$ is minimized by the solution to
$\mathbf{A}\mathbf{x}=\mathbf{b}$.

(a) Quadratic form (81)

(b) Linear System $\mathbf{A}\mathbf{x}=\mathbf{b}$

Figure 5: $\mathbf{A}$ is nonsymmetric. In Figure 5(a), we can find TGCR (1)
still converges fast because the Hessian of Equation 81 is
$\frac{1}{2}(A^{T}+A)$, which is still symmetric. In Figure 5(b), we can find
TGCR(100) has the same convergence rate with GMRES. This shows that TGCR(m)
can still converge for solving nonsymmetric linear systems with $m>1$ and is
mathematically equivalent to non-restart GMRES when $m$ is equal to the matrix
size.

###  C.4 Investigation of Stochastic NLTGCR

In the main paper, we report the effectiveness of NLTGCR for softmax
regression on MNIST. We give the result about using deterministic and
stochastic gradients in Figure 2 and find that NLTGCR only requires a small
batch size and table size ($1$) . In this section, we further investigate the
effectiveness of NLTGCR in a stochastic setting, provide additional
experimental results in Figure 6 and Figure 7. From Figure 6, we can observe
that using the same batch size and table size $m=1$, stochastic NLTGCR
consistently outperforms stochastic AA. Also, it can be observed from Figure 7
that stochastic NLTGCR outperforms stochastic AA for different table size
using a fixed batch size of 1500. Although the short-term property does not
strictly hold in a stochastic setting, we can still observe that NLTGCR
outperforms AA with a smaller variance. In addition, it is worth noting that a
smaller table size works better for both NLTGCR and AA. We suspect this is due
to the accumulation of inaccurate gradient estimates.

(a) Loss

(b) Test Accuracy

Figure 6: Softmax Regression: Effects of batch size, $m=1$

(a)  Loss

(b) Test Accuracy

Figure 7: Softmax Regression: Effects of table size in a stochastic setting,
$B=1000$

##### Compatible with Momentum

Another important technique in optimization is momentum, which speeds up
convergence significantly both in theory and in practice. We experimentally
show that it is possible to further accelerate the convergence of NLTGCR by
using Momentum. We run stochastic NLTGCR with different momentum term and
present results in Figure 8. It suggests that by incorporating momentum into
NLTGCR momentum further accelerates the convergence, although the variance
gets larger for a large momentum term $v$. We leave the theoretical analysis
of of NLTGCR with momentum for future work.

(a)  Loss

(b) Test Accuracy

Figure 8: Softmax Regression: Compatibile with momentum. We further test the
acceleration effect of momentum on stochastic NLTGCR. It shows that stochastic
NLTGCR with momentum converges faster than the one without momentum.

###  C.5 Results for Deep learning applications

####  C.5.1 Image classification using CNN

In a more realistic setting, we test our algorithm for neural networks on an
image classification task. Particularly, we use the standard MNIST dataset
222http://yann.lecun.com/exdb/mnist/. The architecture of the network is based
on the official PyTorch implementation 333implementation
https://github.com/pytorch/examples/blob/master/mnist.. We tried our best to
ensure that the baselines had the best performance in the tests.
Hyperparamters are selected after grid search. We use a batch size of $64$.
For SGD, Nestrov ($v=0.9$), Adam (default $\beta_{1}$ and $\beta_{2}$), and
NLTGCR $(m=1)$, we use a learning rate of $1\times 10^{-2}$, $1\times
10^{-1}$, $1\times 10^{-3}$, and $1\times 10^{-3}$, respectively. Figure 9(a)
shows the curves of loss for training the neural network on MNIST. Figure 9(b)
shows the curves of test accuracy on MNIST. It can be found that NLTGCR
outperforms SGD and Nestrov and is comparable to Adam. In addition, we conduct
experiments on the effects of table size for NLTGCR and present results in
Figure 10. Although Figure 10(b), shows $m=10$ does slightly better, we found
that $m=1$ generally yields robust results. In addition, $m=1$ significantly
reduces the memory and computation overhead. As a result, we would like to
suggest $m=1$ for general experiments. These preliminary results provide
insights of the effectiveness of our algorithm for training neural networks.
Our algorithm is comparable with the widely used optimizer, Adam. In addition,
our algorithm is more memory and computation efficient than other nonlinear
acceleration methods including AA and RNA. As a result, it is worth
investigating the performance of NLTGCR for different tasks and more complex
networks. We leave it for our future work.

(a) Loss

(b) Test Accuracy

Figure 9: Training on MNIST. Averaged on 5 runs, $m=1$ for our algorithm.

(a) Loss

(b) Test Accuracy

Figure 10: Effects of table size $m$ We run experiments with a fixed seed. It
shows S-NLTGCR(10) does slightly better than Adam and S-NLTGCR($m<10)$. We
would suggest to use $m=1$ for saving memory and computation.

####  C.5.2 Image classification using ResNet

We now perform our tests on ResNet32 [28]
444https://github.com/akamaster/pytorch_resnet_cifar10 using CIFAR10 [37]. We
randomly split the training set of all the datasets into two subsets, train
and validation. The former is used to train the neural network, whereas the
latter is used for measuring the performance of the learned model.
Hyperparameters are selected after a grid search. We use a batch size of
$128$. For Nesterov ($v=0.9$), Adam (default $\beta_{1}$ and $\beta_{2}$), and
NLTGCR $(m=1)$, we use a learning rate of $3\times 10^{-4}$, $1\times
10^{-3}$, and $1\times 10^{-1}$, respectively. For better visualization,
figure 11 shows the curves of training loss and validation accuracy using 50
epochs. It can be observed that nlTGCR(1) consistently outperforms baselines
and has a smaller variance.

(a) Training MSE

(b) Validation Accuracy

Figure 11: Validation MSE of training ResNet32 on CIFAR10. Averaged on 5 runs
(manual random seed 0 to 5 for all methods), $m=1$ for our algorithm.
nlTGCR(1) consistently outperforms baselines and has a smaller variance.

####  C.5.3 Time series forecasting using LSTM

Next, we test our algorithm for Long Short-Term Memory [32] on time series
forecasting task using Airplane Passengers and Shampoo Sales Dataset
555https://github.com/spdin/time-series-prediction-lstm-pytorch. We use a
learning rate of 0.04 and the mean squared error (MSE) as our loss function
and evaluation metric. Figure 12 depicts the MSE on validation set during
training. It shows TGCR converges better than baselines. It also suggests TGCR
is capable of optimizing complex deep learning architectures.

(a) Airline-passengers

(b) shampoo

Figure 12: Validation MSE of LSTM on two datasets. Averaged on 5 runs, $m=1$
for our algorithm. It can be observed that nlTGCR(1) outperforms baselines.

####  C.5.4 Semi-supervised classification of graph data using GCN

We also test the effectiveness of nlTGCR on GCN [36] using Cora Dataset. It
consists of 2708 scientific publications classified into one of seven
different classes. The citation network consists of 5429 links. Each
publication in the dataset is described by a 0/1-valued word vector indicating
the absence/presence of the corresponding word from the dictionary. The
objective is to accurately predict the subject of a paper given its words and
citation network, as known as node classification. Hyperparameters are
selected after a grid search. For Nesterov ($v=0.9$), Adam (default
$\beta_{1}$ and $\beta_{2}$), and NLTGCR $(m=1)$, we use a learning rate of
$1\times 10^{-1}$, $1\times 10^{-2}$, and $1\times 10^{-1}$, respectively.
Figure 13(a) and 13(b) depict the training loss and validation accuracy
averaged on 5 runs. It shows TGCR converges faster and achieves higher
accuracy than baselines, which demonstrates the effectiveness of TGCR
optimizing complex deep learning architectures.

(a) Train Loss

(b) Validation Accuracy

Figure 13: Validation accuracy on CORA using GCN. Averaged on 5 runs, $m=1$
for our algorithm. It can be observed that nlTGCR(1) significantly outperforms
baselines.