# Analysis of an Adaptive Safeguarded Newton-Anderson Algorithm with
Applications to Fluid Problems

Matt Dallas
Department of Mathematics, University of Florida
Sara Pollock
Department of Mathematics, University of Florida
Leo Rebholz
Department of Mathematics, Clemson University

###### Abstract

The purpose of this paper is to develop a practical strategy to accelerate
Newton’s method in the vicinity of singular points. We do this by developing
an adaptive safeguarding scheme, which we call $\gamma$-safeguarding, that one
can apply to Anderson accelerated Newton’s method when solving problems near
singular points. The key features of adaptive $\gamma$-safeguarding are that
it converges locally for singular problems, and it can detect nonsingular
problems, in which case the Newton-Anderson iterates are scaled towards a
standard Newton step. This leads to faster local convergence compared to both
Newton’s method and Newton-Anderson without safeguarding, at no additional
computational cost. We demonstrate three strategies one can use when
implementing Newton-Anderson and $\gamma$-safeguarded Newton-Anderson to solve
parameter-dependent problems near singular points. For our benchmark problems,
we take two parameter-dependent incompressible flow systems: flow in a channel
and Rayleigh-Bénard convection.

##  1 Introduction

Nonlinear systems of equations of the form $f(x)=0$ arise frequently in the
sciences, e.g., when computing numerical solutions to nonlinear PDEs. Many of
these systems are parameter-dependent, and exhibit bifurcations for particular
parameter values. A consequence of the Implicit Function Theorem is that the
Jacobian is necessarily singular at bifurcation points [23, p. 8]. Any point
$x^{*}$ at which $f^{\prime}(x^{*})$ is singular is referred to as a singular
point. Hence bifurcation points are necessarily singular points. Solving
problems efficiently near bifurcation points can be challenging, with popular
methods such as Newton’s method failing to converge or converging slowly. Here
we introduce a technique that can both accelerate convergence of Newton’s
method, and recover convergence when Newton’s method fails, at no additional
computational cost.

Standard methods for solving problems at or near singular points include
continuation methods [35] and deflation [4, 13]. Continuation methods are
computationally expensive as they require solving a sequence of nonlinear
problems, varying a parameter between solves. With deflation, one modifies the
nonlinear system in question by penalizing against previously identified
solutions, thereby identifying new ones. Newton’s method can be used to solve
this modified system [4], but can take many iterations to converge. The slow
convergence of Newton’s method, if it converges at all, at and near singular
points motivates the development of techniques to accelerate Newton solves in
the context of bifurcation problems.

Newton’s method applied to singular problems has been studied in great detail
[8, 9, 16, 31, 32], and various acceleration schemes have been proposed and
analyzed to improve convergence in this setting [2, 3, 8, 10, 20, 21, 22, 33].
The method of interest in this paper is Anderson acceleration. Also known as
Anderson extrapolation or Anderson mixing [36], Anderson acceleration was
first introduced in 1965 by D.G. Anderson in [1] to improve the convergence of
fixed-point iterations applied to integral equations. The algorithm combines
the previous $m$ iterates and update steps into a new iterate at each step of
the solve. The number $m$ is commonly known as the algorithmic depth. The
combination of the $m$ iterates often involves solving a least-squares
problem, but since $m$ is typically small, the computational cost of this step
is in general orders of magnitude less than that of a single linear solve.
There are problems for which taking $m$ much larger can be beneficial [26,
36], and later in this paper we’ll see that increasing $m$ can improve or
recover convergence near bifurcation points. With greater depths, the least-
squares problem may suffer from ill-conditioning if proper care is not taken
in the implementation [27].

In [7], two of the authors developed a convergence and acceleration theory for
Anderson accelerated Newton’s method, referred to as Newton-Anderson, applied
to singular problems. The authors introduced a novel safeguarding scheme
called $\gamma$-safeguarding. Convergence was proven for $\gamma$-safeguarded
Newton-Anderson, and it was observed numerically to perform better or no worse
compared to standard Newton-Anderson, particularly when applied to nonsingular
problems.

The purpose of this paper is to develop an adaptive version of
$\gamma$-safeguarding, and to demonstrate the effectiveness of Newton-Anderson
and adaptive $\gamma$-safeguarded Newton-Anderson at solving parameter-
dependent PDEs near bifurcation points in fluid problems. This new adaptive
scheme is proven to be locally convergent under the same conditions as the
non-adaptive scheme, and is able to detect when the problem is nonsingular,
which enables local quadratic convergence enjoyed by Newton for nonsingular
problems. Such a property is desirable when solving nonlinear problems near
bifurcation points, because even if the problem itself is not singular,
convergence can still be affected if it is close to a singular problem, i.e.,
close to a bifurcation point [11]. One would like to enjoy the benefits of
Newton-Anderson in the preasymptotic regime such as a larger domain of
convergence, but not lose quadratic convergence in the asymptotic regime if
the problem is nonsingular. Indeed, Anderson acceleration can reduce the order
of convergence when applied to superlinearly converging iterations such as
Newton’s method applied to a nonsingular problem [30] . This problem is solved
with adaptive $\gamma$-safeguarded Newton-Anderson at no additional
computational cost. We also show numerically that increasing the algorithmic
depth of the Newton-Anderson algorithm can recover convergence from Newton
near bifurcation points, but only for specific choices of $m$.

The algorithms of interest in this paper, Newton, Newton-Anderson with
algorithmic depth $m$ (NA($m$)), $\gamma$-safeguarded Newton-Anderson
($\gamma$NA($r$)), and adaptive $\gamma$-safeguarded Newton-Anderson
($\operatorname{\gamma\text{NAA}(\hat{r})}$) are defined in algorithms 1, 2,
4, and 5 in appendix A. There should be no confusion between NA($m$) and
$\gamma$NA($r$) or $\operatorname{\gamma\text{NAA}(\hat{r})}$ since
$\gamma$-safeguarding is currently only developed for depth $m=1$.

##  2 The $\gamma$-safeguarding algorithm

For the theory discussed in sections 2 and 3, we take
$f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ to be a $C^{3}$ function, $f(x^{*})=0$,
$N=\text{null}\left(f^{\prime}(x^{*})\right)$,
$R=\text{range}\left(f^{\prime}(x^{*})\right)$, $R\perp N$, and $\dim N=1$.
The assumption that $R\perp N$ is common in the singular Newton literature [8,
15], and no generality is lost in finite dimensions. If
$f^{\prime}(x^{*})=U\Sigma V^{T}$ is a singular value decomposition of
$f^{\prime}(x^{*})$, the Jacobian of $F(x):=U^{T}f(x)V$ at $x^{*}$ satisfies
$R\perp N$. Since Newton’s method is essentially invariant under nonsingular
affine transformations of the domain and nonsingular linear transformations of
the range [12], it suffices to study $F(x)$. Thus we can assume that $R\perp
N$ without loss of generality. We let $P_{N}$ and $P_{R}$ denote the
orthogonal projections onto $N$ and $R$ respectively. Lastly, let
$e_{k}=x_{k}-x^{*}$, $w_{k+1}=-f^{\prime}(x_{k})^{-1}f(x_{k})$, and

$\displaystyle\theta_{k+1}=\frac{\|w_{k+1}-\gamma_{k+1}(w_{k+1}-w_{k})\|}{\|w_{k+1}\|},$
(1)

where $\gamma_{k+1}$ is computed via algorithm 2. The term $\theta_{k+1}$ in
(1) is referred to in the Anderson literature as the optimization gain, and is
key to determining when Anderson acceleration is successful both in the
singular and nonsingular cases [7, 26]. In the $\gamma$-safeguarding
algorithm, algorithm 4, the NA steps are scaled towards a Newton step based on
a user-chosen parameter, $r$. This parameter determines how strongly the NA
steps are scaled towards a Newton step. When $r\approx 0$, the safeguarded NA
steps will be heavily scaled towards the latest Newton step, and when
$r\approx 1$, the iterates will behave more like standard NA. This
safeguarding strategy was introduced in [7] to deal with some of the technical
challenges involved with proving convergence of Newton-like methods near
singular points. Specifically, $\gamma$-safeguarding ensures that an Anderson-
accelerated Newton step is not accelerated so much that the iterate leaves the
region of invertibility. As discussed in [7] and references therein, Newton
iterates near the null space $N$ will remain near $N$. This is due to linear
convergence of $P_{N}e_{k}$ and the quadratic convergence of $P_{R}e_{k}$,
from which one can show, under appropriate conditions [8], that there exists a
constant $\hat{\sigma}>0$ for which

$\displaystyle\frac{\|P_{R}e_{k}\|}{\|P_{N}e_{k}\|}:=\sigma_{k}\leq\hat{\sigma}$
(2)

for all $k$ provided $\|e_{0}\|$ and $\|\sigma_{0}\|$ are sufficiently small.
Thus accelerating Newton means accelerating the component of the iterates
along the null space. The risk is that as $P_{N}e_{k}$ is accelerated,
$\sigma_{k}$ increases, potentially greater than $\hat{\sigma}$. In the case
of $\dim N=1$, there is potential for Anderson to accelerate Newton
significantly, especially in the asymptotic regime. This is demonstrated by
the following proposition.

######  Proposition 2.1.

Suppose $x_{k+1}=x_{k}+w_{k+1}-\gamma_{k+1}(x_{k}-x_{k-1}+w_{k+1}-w_{k})$ is
well-defined, $w_{k+1}\neq w_{k}$, and $\dim N=1$. Then there is a number
$C_{k}$ depending on $\sigma_{k}$, $\|e_{k}\|$, $\gamma_{k+1}$, $\hat{\gamma}$
(defined in proof) and $f$ such that

$\displaystyle\|P_{N}e_{k+1}\|\leq
C_{k}\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}.$ (3)

The proof relies on the following two results from [7]. The second, lemma 2.2,
says that if $x_{k}$ and $x_{k-1}$ are both near $N$, then $x_{k+1}$ is
compatible. This is a notion developed in [7] to categorize Newton-Anderson
iterates whose acceleration is determined by $\theta_{k+1}$. These are the
iterates in the singular case that behave in a similar way to the iterates in
the nonsingular case [26].

######  Lemma 2.1.

Let $\hat{D}(x)(\cdot):=P_{N}f^{\prime\prime}(x^{*})(x-x^{*},P_{N}(\cdot))$ be
invertible as a linear map on $N$, and let $\|x-x^{*}\|$ be sufficiently small
so that $P_{R}f^{\prime}(x)P_{R}(\cdot)$ is invertible. We will let
$T_{k}(\cdot):=(1/2)\hat{D}^{-1}(x_{k})f^{\prime\prime}(x_{k})(e_{k},\cdot)$,
and $I$ denote the identity map on $\mathbb{R}^{n}$. Then

$\displaystyle
w_{k+1}=-\frac{1}{2}P_{N}e_{k}+\left(T_{k}-I\right)P_{R}e_{k}+\mathcal{O}(\|e_{k}\|^{2}).$
(4)

Note that the range of $T_{k}(\cdot)$ lies in $N$. Applying (4) to
$w_{k+1}^{\alpha}:=w_{k+1}-\gamma_{k+1}(w_{k+1}-w_{k})$, and letting
$q_{k-1}^{k}$ denote a term where $\|q_{k-1}^{k}\|\leq
c_{1}|1-\gamma_{k+1}|\,\|e_{k}\|^{2}+c_{2}|\gamma_{k+1}|\,\|e_{k-1}\|^{2}$, we
get

$\displaystyle
w_{k+1}^{\alpha}=-\frac{1}{2}P_{N}e_{k}^{\alpha}+\left(\left(T_{k}-I\right)P_{R}e_{k}\right)^{\alpha}+q_{k-1}^{k}.$
(5)

######  Lemma 2.2.

Let the assumptions of lemma 2.1 hold, and let
$\sigma_{k-1,k}=\max\\{\sigma_{k-1},\sigma_{k}\\}$ and
$\|e_{k-1,k}\|=\max\\{\|e_{k-1}\|,\|e_{k}\|\\}$. If both $x_{k}$ and $x_{k-1}$
lie sufficiently close to $N$, then

$\displaystyle\|P_{N}e_{k+1}\|\leq\left(\frac{1+C_{1}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}{1-C_{2}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}\right)\theta_{k+1}\|w_{k+1}\|.$
(6)

We remark that $C_{1}$ and $C_{2}$ tend to one as $\|\sigma_{k-1,k}\|$ and
$\|e_{k-1,k}\|$ tend to zero [7]. Now we prove proposition 2.1

###### Proof.

(Proposition 2.1)

Using the optimality of $\gamma_{k+1}$ and orthogonality of $R$ and $N$, we
have that for any $\gamma\in\mathbb{R}$,

$\displaystyle\|w_{k+1}^{\alpha}\|_{2}^{2}\leq\|P_{N}w_{k+1}-\gamma(P_{N}w_{k+1}-P_{N}w_{k})\|_{2}^{2}+\|P_{R}w_{k+1}-\gamma(P_{R}w_{k+1}-P_{R}w_{k})\|_{2}^{2}.$
(7)

Observe that $\|P_{N}w_{k+1}-\gamma(P_{N}w_{k+1}-P_{N}w_{k})\|_{2}^{2}$ is a
one-dimensional expression. So we can pick $\gamma=\hat{\gamma}$ such that
$\|P_{N}w_{k+1}-\hat{\gamma}(P_{N}w_{k+1}-P_{N}w_{k})\|_{2}^{2}=0$. It then
follows from lemmas 2.1 and 2.2 that

$\displaystyle\|P_{N}e_{k+1}\|$
$\displaystyle\leq\left(\frac{1+C_{1}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}{1-C_{2}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}\right)\theta_{k+1}\|w_{k+1}\|$
(8)
$\displaystyle=\left(\frac{1+C_{1}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}{1-C_{2}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}\right)\|w_{k+1}^{\alpha}\|$
(9)
$\displaystyle\leq\left(\frac{1+C_{1}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}{1-C_{2}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}\right)q_{k-1}^{k}.$
(10)

Note that the $q_{k-1}^{k}$ in the last line comes from 5, and that
$\|q_{k-1}^{k}\|\leq
c_{1}|1-\hat{\gamma}|\,\|e_{k}\|^{2}+c_{2}|\hat{\gamma}|\,\|e_{k-1}\|^{2}$
here. Taking

$\displaystyle
C_{k}:=\left(\frac{1+C_{1}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}{1-C_{2}\left(\sigma_{k-1,k},\|e_{k-1,k}\|\right)}\right)\max\\{c_{1}|1-\hat{\gamma}|,c_{2}|\hat{\gamma}|\\}$
(11)

completes the proof. ∎

Remark: Note that $\hat{\gamma}$ can be very large when $P_{N}w_{k+1}\approx
P_{N}w_{k}$, which may correspond to the terminal phase of the solve. So this
bound is meaningful in the asymptotic regime when there is still a significant
decrease in the residual at each step.

Such a bound is good for a single step, but as discussed above, such dramatic
acceleration of $P_{N}e_{k+1}$ could place $x_{k+1}$ outside the domain of
invertibility. This is where $\gamma$-safeguarding is useful. In [7], it is
shown that

######  Lemma 2.3.

For $\lambda\in(0,1]$, let

$\displaystyle\nu_{k+1}=\frac{\min\\{|1-\lambda\gamma_{k+1}|\,\|P_{N}(e_{k}+w_{k+1})\|,|\lambda\gamma_{k+1}|\,\|P_{N}(e_{k-1}+w_{k})\|\\}}{\max\\{|1-\lambda\gamma_{k+1}|\,\|P_{N}(e_{k}+w_{k+1})\|,|\lambda\gamma_{k+1}|\,\|P_{N}(e_{k-1}+w_{k})\|\\}}.$
(12)

Fix $0<r<1$. Assume $|\gamma_{k+1}|<1$ and nonzero. Then given
$\sigma_{i}<\hat{\sigma}$ and $\rho_{i}<\hat{\rho}$, $i=k,k-1$, there exists a
number $\lambda\in(0,1]$ such that $\nu_{k+1}\leq r<1$ and
$\min\\{|1-\lambda\gamma_{k+1}|\,\|P_{N}(e_{k}+w_{k+1}),|\lambda\gamma_{k+1}|\,\|P_{N}(e_{k-1}+w_{k})\|\\}=|\lambda\gamma_{k+1}|\,\|P_{N}(e_{k-1}+w_{k})\|$.

We are also able to bound

$\displaystyle\overset{\sim}{\nu}_{k+1}=\frac{\min\\{|1-\lambda\gamma_{k+1}|\,\|P_{N}e_{k}\|,|\lambda\gamma_{k+1}|\,\|P_{N}e_{k-1}\|\\}}{\max\\{|1-\lambda\gamma_{k+1}|\,\|P_{N}e_{k}\|,|\lambda\gamma_{k+1}|\,\|P_{N}e_{k-1}\|\\}}$
(13)

by $\gamma$-safeguarding. This leads to the local convergence of
$\gamma$-Newton-Anderson presented in [7]. It was observed in [7] that
$\gamma$-safeguarding performs well compared to Newton-Anderson for singular
problems, and it could outperform Newton-Anderson when applied to nonsingular
problems with the right choice of parameter $r$. This motivated the
development of an adaptive form of $\gamma$-Newton-Anderson. This is the focus
of the next section.

##  3 Adaptive $\gamma$-safeguarding

Algorithm 4 is very useful in proving local convergence of $\gamma$NA($r$),
and it performed competitively with standard NA(1) in [7]. For example,
$\gamma$NA(0.5) could outperform NA(1) when applied to certain nonsingular
problems. This is not surprising. It had previously been observed in [28] that
NA does not necessarily improve convergence when applied to nonsingular
problems, and recent results in [30] show that under certain conditions NA
reduces the order of a convergence for superlinearly converging sequences.
Thus, by setting the $\gamma$-safeguarding parameter $r=0.5$, thereby scaling
the iterates closer to Newton steps, we see less of the effect of a full NA
step on the order of convergence.

One may not always know a priori if a given problem is singular or
nonsingular, and even if this is known, setting a fixed $r$ means that the
effect of NA is never completely eliminated which could lead to a smaller
order of convergence relative to Newton for nonsingular problems. This
motivates the development of an adaptive form of $\gamma$NA($r$) that can
detect nonsingular problems and scale an NA step accordingly. We will denote
this adaptive choice of $r$ from algorithm 4 by $r_{k}$, with $k$ denoting the
iteration count. For reference, we record three criteria that a satisfactory
choice of $r_{k}$ should satisfy.

######  Criteria 3.1.

An adaptive $\gamma$-safeguarding tolerance $r_{k}$ should satisfy the
following.

  1. 1.

$r_{k}{\ll}1$ if $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|{\ll}1$;

  2. 2.

$r_{k}\approx 1$ if $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|\approx 1$; and

  3. 3.

$\lim_{k\to\infty}r_{k}=0$ if $f^{\prime}(x^{*})$ is nonsingular.

Criterion (1) says that if $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|$ is very small,
then we want to scale the NA step generated from $x_{k}$ and $x_{k-1}$ heavily
towards $x_{k}+w_{k+1}$. In the singular case, this will (locally) keep
$x_{k+1}$ within the domain of invertibility. Alternatively, if the problem is
nonsingular, but close to a singular problem, then scaling NA towards a Newton
step is also preferred when $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|$ is small since
then we do not slow Newton’s fast local quadratic convergence. Criterion (2)
says that if $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|$ is close to one, then the error
is not decreasing significantly, and we want to allow NA to act on $x_{k}$ and
$x_{k-1}$ without significant scaling of $\gamma_{k+1}$ from safeguarding.
Lastly, criterion (3) is important because if $f^{\prime}(x^{*})$ is
nonsingular, Newton’s method will converge quadratically in a neighborhood of
$x^{*}$. We therefore want to “turn off” NA near $x^{*}$, and insisting that
$r_{k}\to 0$ asymptotically achieves this.

Observing criteria (1)-(3), one may note that we essentially want $r_{k}$ to
behave like $\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|$ within the domain of
convergence. Of course, we can not actually compute
$\|P_{N}e_{k}\|/\|P_{N}e_{k-1}\|$, but lemma 2.1 says that in the singular
domain of convergence $w_{k+1}\approx P_{N}e_{k}$. So if we take
$r_{k}=\|w_{k+1}\|/\|w_{k}\|$, we can expect criteria (1)-(2) to be enforced
locally. For (3), if $f^{\prime}(x^{*})$ is nonsingular, then locally we will
have $\|w_{k+1}\|/\|w_{k}\|\to 0$. Thus $\|w_{k+1}\|/\|w_{k}\|$ satisfies the
three criteria within the domain of convergence.

Henceforth we will let $\eta_{k+1}=\|w_{k+1}\|/\|w_{k}\|$. Adaptive
$\gamma$-safeguarding differs by algorithm 4 only in the first line. The user
still sets an $\hat{r}\in(0,1)$, but rather than use this $r$ in subsequent
steps, adaptive $\gamma$-safeguarding uses
$r_{k}=\min\\{\eta_{k+1},\hat{r}\\}$. The choice of $\hat{r}$ here sets the
weakest safeguarding the user wants to impose. Thus we are always safeguarding
at least as strictly as standard $\gamma$-safeguarding with $r=\hat{r}$.
Stated concisely, we have $r_{k}\leq\hat{r}$. Local convergence of
$\operatorname{\gamma\text{NAA}(\hat{r})}$ then follows from theorem $6.1$ in
[7]. To state this precisely, let $\lambda_{k+1}$ be the value computed by
algorithm 5 at step $k+1$,
$\theta_{k+1}^{\lambda}=\|w_{k+1}-\lambda_{k+1}\gamma_{k+1}(w_{k+1}-w_{k})\|/\|w_{k+1}\|$,
and
$L_{k+1}^{\lambda}(x,y)=(1-\lambda_{k+1}\gamma_{k+1})x+\lambda_{k+1}\gamma_{k+1}y$.
Then we have the following theorem.

######  Theorem 3.1.

Let $\dim N=1$, and let $\hat{D}_{N}(x)$ be invertible as a map on $N$ for all
$P_{N}(x-x^{*})\neq 0$. Let $W_{k}=W(\|e_{k}\|,\sigma_{k},x^{*})$. If $x_{0}$
is chosen so that $\sigma_{0}<\hat{\sigma}$ and $\|e_{0}\|<\hat{\rho}$, for
sufficiently small $\hat{\sigma}$ and $\hat{\rho}$, $x_{1}=x_{0}+w_{1}$, and
$x_{k+1}=L_{k+1}^{\lambda}(x_{k}+w_{k+1},x_{k-1}+w_{k})$ for $k\geq 1$, then
$W_{k+1}\subset W_{0}$ for all $k\geq 0$ and $x_{k}\to x^{*}$. That is,
$\\{x_{k}\\}$ remains well-defined and converges to $x^{*}$. Furthermore,
there exist constants $c_{4}>0$ and $\kappa\in(1/2,1)$ such that

$\displaystyle\|P_{R}e_{k+1}\|$ $\displaystyle\leq
c_{4}\max\\{|1-\lambda_{k+1}\gamma_{k+1}|\,\|e_{k}\|^{2},|\lambda_{k+1}\gamma_{k+1}|\,\|e_{k-1}\|^{2}\\}$
(14) $\displaystyle\|P_{N}e_{k+1}\|$
$\displaystyle<\kappa\theta_{k+1}^{\lambda}\|P_{N}e_{k}\|$ (15)

for all $k\geq 1$.

###### Proof.

As noted in the paragraph before theorem 3.1, the local convergence of
$\operatorname{\gamma\text{NAA}(\hat{r})}$ when $f^{\prime}(x^{*})$ is
singular with $\dim N=1$ and $\hat{D}_{N}(x)$ be invertible as a map on $N$
for all $P_{N}(x-x^{*})\neq 0$ follows from theorem $6.1$ in [7]. The key role
$\gamma$-safeguarding played in the proof was to ensure that that $\nu_{k+1}$
and $\overset{\sim}{\nu}_{k+1}$ are bounded above by some $r<1$. This is the
parameter set by the user in algorithm 4. With adaptive $\gamma$-safeguarding,
we have that the adaptive parameter $r_{k}\leq\hat{r}<1$, and therefore
convergence follows from the same argument as in [7]. ∎

Our choice of $r_{k}$ in algorithm 5 is partially motivated by criterion (3)
in 3.1. That is, in the case of a nonsingular problem, we prefer to use Newton
asymptotically rather than Newton-Anderson. Hence we want $r_{k}$ to tend to
zero as our solver converges, thereby scaling the iterates heavily towards
pure Newton steps in the asymptotic regime and enjoying quadratic convergence
locally. The remainder of this section is dedicated to quantifying how close
an iterate generated by algorithm 5 is to a standard Newton iterate in the
asymptotic regime when $f^{\prime}(x^{*})$ is nonsingular. For the remainder
of this section, we take the asymptotic regime to be all $x_{k}$ such that
$\eta_{k+1}<1$.

The main result is theorem 3.2 below which bounds
$\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|$ locally, where
$x_{k+1}^{\operatorname{Newt}}:=x_{k}+w_{k+1}$, and
$x_{k+1}^{\operatorname{NA}}:=x_{k+1}^{\operatorname{Newt}}-\lambda_{k+1}\gamma_{k+1}(x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}})$.
This notation is introduced to emphasize the Newton and Newton-Anderson
iterates in the comparison. We will also define
$e_{k}^{\operatorname{Newt}}:=x_{k}^{\operatorname{Newt}}-x^{*}=x_{k-1}+w_{k}-x^{*}$.
The following lemma will be used in the proof of theorem 3.2. Lemma 3.1 bounds
$|\lambda_{k+1}\gamma_{k+1}|$, the scaled $\gamma_{k+1}$ returned by
$\operatorname{\gamma\text{NAA}(\hat{r})}$, in terms of $\eta_{k+1}$ and
$r_{k}=\min\\{\eta_{k+1},\hat{r}\\}$.

######  Lemma 3.1.

Let $\eta_{k+1}=\|w_{k+1}\|/\|w_{k}\|$, $\lambda_{k+1}$ be the value computed
by $\operatorname{\gamma\text{NAA}(\hat{r})}$. Set $\hat{r}\in(0,1)$, define
$r_{k}:=\min\\{\eta_{k+1},\hat{r}\\}$, and $\beta_{k+1}:=r_{k}\eta_{k+1}$. If
$\eta_{k+1}<1$, then the value $\lambda_{k+1}\gamma_{k+1}$ returned by
$\operatorname{\gamma\text{NAA}(\hat{r})}$ (algorithm 5) satisfies

$\displaystyle|\lambda_{k+1}\gamma_{k+1}|\leq\frac{\beta_{k+1}}{1-\beta_{k+1}}.$
(16)

when $\lambda_{k+1}=1$, and

$\displaystyle|\lambda_{k+1}\gamma_{k+1}|=\frac{\beta_{k+1}}{1+\operatorname{sign}(\gamma_{k+1})\beta_{k+1}}.$
(17)

when $\lambda_{k+1}<1$.

###### Proof.

The proof is essentially just walking through the cases in
$\operatorname{\gamma\text{NAA}(\hat{r})}$.
$\operatorname{\gamma\text{NAA}(\hat{r})}$, algorithm 5, returns
$\lambda_{k+1}\gamma_{k+1}$, where $\lambda_{k+1}$ is computed based on
whether $|\gamma_{k+1}|/|1-\gamma_{k+1}|\leq\beta_{k+1}$ or if
$|\gamma_{k+1}|/|1-\gamma_{k+1}|>\beta_{k+1}$. In the former,
$\lambda_{k+1}=1$ and the value of $\gamma_{k+1}$ is unchanged. Thus

$\displaystyle|\lambda_{k+1}\gamma_{k+1}|=|\gamma_{k+1}|\leq\beta_{k+1}|1-\gamma_{k+1}|\leq\beta_{k+1}+\beta_{k+1}|\gamma_{k+1}|.$
(18)

Now (16) follows from (18). Now suppose
$|\gamma_{k+1}|/|1-\gamma_{k+1}|>\beta_{k+1}$. Then one can show that either
$\beta_{k+1}/(\gamma_{k+1}(1+\beta_{k+1}))<1$ or
$\beta_{k+1}/(\gamma_{k+1}(\beta_{k+1}-1))<1$ depending on whether
$\gamma_{k+1}$ is positive or negative respectively. Equivalently, either the
condition in line 7 or in line 10 of algorithm 5 is met when
$|\gamma_{k+1}|/|1-\gamma_{k+1}|>\beta_{k+1}$. These conditions imply
$\lambda_{k+1}\gamma_{k+1}=\beta_{k+1}/(\beta_{k+1}+1)$ when $\gamma_{k+1}>0$
and $\lambda_{k+1}\gamma_{k+1}=\beta_{k+1}/(\beta_{k+1}-1)$ when
$\gamma_{k+1}<0$. Taking absolute values, we can concisely express both cases
at once by

$\displaystyle|\lambda_{k+1}\gamma_{k+1}|=\frac{\beta_{k+1}}{1+\operatorname{sign}(\gamma_{k+1})\beta_{k+1}},$
(19)

which proves (17). ∎

Recall that we set
$x_{k+1}^{\operatorname{NA}}=x_{k+1}^{{\operatorname{Newt}}}-\lambda_{k+1}\gamma_{k+1}(x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}})$.
That is, $x_{k+1}^{\operatorname{NA}}$ is the iterate returned by
$\operatorname{\gamma\text{NAA}(\hat{r})}$. Using lemma 3.1, we can bound
$\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|$ in terms of
$\|e_{k}\|$, $\|e_{k-1}\|$, and $\eta_{k+1}$.

######  Theorem 3.2.

If $f^{\prime}(x^{*})$ is nonsingular, then there exists an $R>0$ and a
constant $C$ depending only on $f$ such that for $x_{k}$ and $x_{k-1}$ in
$B_{R}(x^{*})$ and $\eta_{k+1}<1$,

$\displaystyle\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|\leq
C\left(\frac{\beta_{k+1}}{1-\beta_{k+1}}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$
(20)

when $\lambda_{k+1}=1$, and

$\displaystyle\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|\leq
C\left(\frac{\beta_{k+1}}{1+\operatorname{sign}(\gamma_{k+1})\beta_{k+1}}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$
(21)

when $\lambda_{k+1}<1$.

###### Proof.

Using our notation from the discussion preceding lemma 3.1, an iterate
generated by $\operatorname{\gamma\text{NAA}(\hat{r})}$ takes the form
$x_{k+1}^{\operatorname{NA}}=x_{k+1}^{\operatorname{Newt}}-\lambda_{k+1}\gamma_{k+1}\left(x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}}\right)$.
Therefore
$\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|=|\lambda_{k+1}\gamma_{k+1}|\,\|x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}}\|$.
Since $f^{\prime}(x^{*})$ is nonsingular, we can take $R$ sufficiently small
to ensure that $\|e_{j+1}^{\operatorname{Newt}}\|\leq C\|e_{j}\|^{2}$, where
$C$ is a constant determined by $f$, when $\|e_{j}\|<R$ for $j=k,k-1$. Hence,
upon adding and subtracing $x^{*}$ we obtain
$\|x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}}\|=\|e_{k+1}^{\operatorname{Newt}}-e_{k}^{\operatorname{Newt}}\|\leq
2C\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$. To complete the proof, we write
$2C=C$ and apply lemma 3.1 to bound $|\lambda_{k+1}\gamma_{k+1}|$ for the
cases $\lambda_{k+1}=1$ and $\lambda_{k+1}<1$. ∎

We conclude this section by discussing two noteworthy corollaries following
from theorem 3.2. The first corollary we consider bounds the error
$\|e_{k+1}\|$ of an iterate $x_{k+1}$ generated by
$\operatorname{\gamma\text{NAA}(\hat{r})}$ in terms of $\|e_{k}\|$,
$\|e_{k-1}\|$, and $\beta_{k+1}$.

######  Corollary 3.1.

Suppose $f^{\prime}(x^{*})$ is nonsingular and
$x_{k+1}=x_{k+1}^{\operatorname{Newt}}-\lambda_{k+1}\gamma_{k+1}(x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}})$.
If $x_{k}$ and $x_{k-1}$ are sufficiently close to $x^{*}$ with
$\eta_{k+1}<1$, then there are constants $C_{1}$ and $C_{2}$ depending only on
$f$ such that

$\displaystyle\|e_{k+1}\|\leq\left(C_{1}\left(\frac{\beta_{k+1}}{1-\beta_{k+1}}\right)+C_{2}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$
(22)

when $\lambda_{k+1}=1$ and

$\displaystyle\|e_{k+1}\|\leq\left(C_{1}\left(\frac{\beta_{k+1}}{1+\operatorname{sign}(\gamma_{k+1})\beta_{k+1}}\right)+C_{2}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$
(23)

when $\lambda_{k+1}<1$.

###### Proof.

The proof is exactly the same for two cases $\lambda_{k+1}=1$ and
$\lambda_{k+1}<1$. We’ll therefore assume $\lambda_{k+1}<1$. Note that

$\displaystyle\|e_{k+1}\|\leq\|e_{k+1}-e_{k+1}^{\operatorname{Newt}}\|+\|e_{k+1}^{\operatorname{Newt}}\|=\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|+\|e_{k+1}^{\operatorname{Newt}}\|.$
(24)

Since the problem is nonsingular, $\|e_{k+1}^{\operatorname{Newt}}\|\leq
C_{2}\|e_{k}\|^{2}\leq C_{2}\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$ for
sufficiently small $\|e_{k}\|$. To complete the proof, apply theorem 3.2 to
$\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|$ to obtain

$\displaystyle\|e_{k+1}\|\leq
C_{1}\left(\frac{\beta_{k+1}}{1+\operatorname{sign}(\gamma_{k+1})\beta_{k+1}}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}+C_{2}\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}.$
(25)

∎

Considering the previous results, one may note that, locally, the bounds in
theorem 3.2 and corollary 3.1 are particularly strong when $\gamma_{k+1}>0$,
which occurs when $w_{k+1}^{T}(w_{k+1}-w_{k})>0$, or, equivalently, when the
direction cosine between $w_{k+1}$ and $w_{k+1}-w_{k}$ is positive. This is
the case highlighted in the following corollary. Its proof is immediate from
theorem 3.2 and corollary 3.1.

######  Corollary 3.2.

Let the assumptions of theorem 3.2 and corollary 3.1 hold. If
$\gamma_{k+1}>0$, then the following bounds hold

$\displaystyle\|x_{k+1}^{\operatorname{NA}}-x_{k+1}^{\operatorname{Newt}}\|$
$\displaystyle\leq C\beta_{k+1}\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$ (26)
$\displaystyle\|e_{k+1}\|$
$\displaystyle\leq\left(C_{2}\beta_{k+1}+C_{1}\right)\max\\{\|e_{k}\|^{2},\|e_{k-1}\|^{2}\\}$
(27)

##  4 Numerics

In this section, we demonstrate the effectiveness of NA and
$\operatorname{\gamma\text{NAA}(\hat{r})}$ near bifurcation points by applying
these algorithms to the following parameter-dependent PDEs. All computations
are performed on an M1 MacBook with GNU Octave 8.2.0.

###  4.1 Test Problems

  1. 1.

Navier-Stokes Flow in a Channel

${\small\left\\{\begin{split}-\mu\Delta\operatorname{\mathbf{u}}+\operatorname{\mathbf{u}}\cdot\nabla\operatorname{\mathbf{u}}+\nabla
p&=\mathbf{0}\\\
\nabla\cdot\operatorname{\mathbf{u}}&=\mathbf{0}\end{split}\qquad\qquad\begin{split}\mathbf{u}&=\mathbf{u}_{\text{in}},\hskip
4.5pt\Gamma_{\text{in}}\\\ \mathbf{u}&=0,\hskip 4.5pt\Gamma_{\text{wall}}\\\
-p\mathbf{n}+(\mu\nabla\mathbf{u})\mathbf{n}&=0,\hskip
4.5pt\Gamma_{\text{out}}\end{split}\right.}$ (28)

  2. 2.

Rayleigh-Bénard Convection

${\small\left\\{\begin{split}-\mu\Delta\operatorname{\mathbf{u}}+\operatorname{\mathbf{u}}\cdot\nabla\operatorname{\mathbf{u}}+\nabla
p-\text{Ri}\,T\mathbf{e}_{y}&=\mathbf{0}\\\
\nabla\cdot\operatorname{\mathbf{u}}&=0\\\ -\kappa\Delta
T+\operatorname{\mathbf{u}}\cdot\nabla T&=0\\\
\end{split}\qquad\qquad\begin{split}T&=1,\hskip
2.25pt\operatorname{\mathbf{x}}\in\Gamma_{1}:=\\{1\\}\times(0,1),\\\
T&=0,\hskip
2.25pt\operatorname{\mathbf{x}}\in\Gamma_{2}:=\\{0\\}\times(0,1),\\\ \nabla
T\cdot\operatorname{\mathbf{n}}&=0,\hskip
2.25pt\operatorname{\mathbf{x}}\in\Gamma_{3}:=(0,1)\times\\{0,1\\},\\\
\operatorname{\mathbf{u}}&=\mathbf{0},\hskip
2.25pt\operatorname{\mathbf{x}}\in\partial\Omega=\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{3}.\end{split}\right.}$
(29)

In both models, $\mathbf{u}$ denotes the fluid velocity, $p$ the pressure, and
$\mathbf{n}$ the outward normal. In (28), $\mu$ denotes the viscosity
parameter. In (29), $T$ denotes the temperature of the fluid, and
$\operatorname{Ri}$ denotes the Richardson number, which is the parameter of
interest for (29). We set $\mu=\kappa=10^{-2}$. These parameters values are
chosen so as to replicate the results seen in [14] as $\operatorname{Ri}$
ranges from $3.0$ to $3.5$. For the flow in a channel, model (28), we use
$\mathcal{P}_{2}-\mathcal{P}_{1}$ Taylor-Hood elements [5, p. 164]. The
channel, shown in figure 1, is arranged such that the left most boundary lies
at $x=0$, the right most boundary lies at $x=50$, and the boundary components
are given by $\Gamma_{\text{in}}=\\{0\\}\times[2.5,5]$,
$\Gamma_{\text{out}}=\\{50\\}\times[0,7.5]$, and
$\Gamma_{\text{wall}}=[0,10]\times(\\{2.5\\}\cup\\{5\\})\cup\\{10\\}\times([0,2.5]\cup[5,7.5])\cup[10,40]\times(\\{0\\}\cup\\{7.5\\})$.
For the Rayleigh-Béndard model, we use
$\mathcal{P}_{2}-\mathcal{P}_{1}^{\text{disc}}$ Scott-Vogelius elements [17],
where $\mathcal{P}_{1}^{\text{disc}}$ denotes piece-wise linear discontinuous
elements. The $\mathcal{P}_{2}-\mathcal{P}_{1}^{\text{disc}}$ elements are
stable on Alfeld-split, also known as barycenter-split, triangulations [29, p.
77]. The meshes used for each model are shown below.

Figure 1:  Top: Mesh used for Rayleigh-Bénard model. Bottom: mesh used for
flow in a channel.

For (28), it is known [25] that there exists a critical viscosity
$\mu^{*}\in(0.9,1)$ at which a bifurcation occurs. For $\mu>\mu^{*}$, the
stable velocity solution is symmetric about the center horizontal ($y=3.75$)
as seen in the top plot of figure 2. For $\mu<\mu^{*}$, there is still a
symmetric solution, but it is unstable. Stability is inherited by two
asymmetric solutions seen in the bottom two plots of figure 2.

Figure 2:  Top: Representative symmetric solution for $\mu=1.0$. Middle:
Representative asymmetric solution with positive vertical velocity upon
exiting the narrow channel for $\mu=0.9$. Bottom: Representative asymmetric
solution with negative vertical velocity upon exiting the narrow channel for
$\mu=0.9$.

The parameter region of interest for (29) is $\operatorname{Ri}\in[3,3.5]$. In
this range, the flow appears to be in transition from a single eddy in the
center of the domain to two eddies as seen in figure 3 below.

Figure 3:  Velocity streamlines for velocity $u$ from model (29) showing
transition from one eddy to two eddies. Top Left: Ri = 3.0. Top Right: Ri =
3.1. Bottom Left: Ri = 3.2. Bottom Right: Ri = 3.4.

For (28), we take the zero vector as our initial iterate. For (29), we
initialize our iterate by applying a Picard step and then applying
$\operatorname{\gamma\text{NAA}(\hat{r})}$. The $H^{1}$ seminorm is used in
our implementations [5].

###  4.2 General Discussion of Results

The following experiments demonstrate three strategies for solving nonlinear
problems near bifurcations using NA and
$\operatorname{\gamma\text{NAA}(\hat{r})}$. The first two are asymptotic
safeguarding and preasymptotic safeguarding. With asymptotic safeguarding,
$\operatorname{\gamma\text{NAA}(\hat{r})}$, algorithm 5, is only activated
when the residual is smaller than some user-chosen threshold. Hence the solve
will behave like NA until the last few iterations when
$\operatorname{\gamma\text{NAA}(\hat{r})}$ is activated. For the experiments
in section 4.3, we activated algorithm 5 when $\|w_{k+1}\|<10^{-1}$. We chose
$10^{-1}$ since we want to activate $\gamma$-safeguarding as early as
possible. This allows for earlier detection of a nonsingular problem, and thus
faster convergence. We found that activating $\gamma$-safeguarding for
$\|w_{k+1}\|<\tau<10^{-1}$ will not necessarily break convergence, but if the
problem is nonsingular, this will not be detected as early. Activating
$\operatorname{\gamma\text{NAA}(\hat{r})}$ when $\|w_{k+1}\|<1$ can break
convergence, though this seems to be problem-dependent. Convergence for (28)
was virtually unaffected with threshold $\|w_{k+1}\|<1$, but for (29) with
$\operatorname{Ri}=3.5$, setting the activation threshold to $1$ caused
$\gamma\text{NAA}(0.9)$ to diverge when it had converged with threshold $0.1$.
With activation threshold $0.1$, asymptotic safeguarding is shown to be
effective close to the bifurcation point, where Newton’s method can fail to
converge. Preasymptotic safeguarding, on the other hand, applies
$\operatorname{\gamma\text{NAA}(\hat{r})}$ at each step of the solve.
Interestingly, we observe that $\operatorname{\gamma\text{NAA}(\hat{r})}$
applied preasymptotically can outperform NA when applied to (29), and even
recover convergence when both Newton and $\operatorname{NA}$ diverge (see
figure 11). A theoretical explanation for this requires a better understanding
of these methods in the preasymptotic regime. In particular, a better
understanding of the descent properties of Anderson acceleration would be of
great value. It is known [34] that for singular problems in the preasymptotic
regime, where $\|f(x)\|$ is not small, the Newton update step
$s_{k}=-f^{\prime}(x_{k})^{-1}f(x_{k})$, or $w_{k+1}$ in our notation, can be
large and nearly orthogonal to the gradient of $\|f\|_{2}^{2}$. With
$\operatorname{NA}$, our update step takes the form
$s_{k}^{\operatorname{NA}}=w_{k+1}-\gamma_{k+1}(x_{k+1}^{\operatorname{Newt}}-x_{k}^{\operatorname{Newt}})$.
It is clear that $s_{k}^{\operatorname{NA}}$ is a descent direction for
sufficiently small $\gamma_{k+1}$, since in this case it is nearly $w_{k+1}$.
It is possible that for certain values of $\gamma_{k+1}$,
$s_{k}^{\operatorname{NA}}$ is a stronger descent direction than $s_{k}$.

Another possible explanation as to why
$\operatorname{\gamma\text{NAA}(\hat{r})}$ can outperform $\operatorname{NA}$
in some cases is its resemblance to restarted Anderson acceleration methods.
Restarted versions of Anderson acceleration are often applied in various forms
for depth $m>1$, and have been shown to effective in practice [6, 18, 19, 24].
In the special case of Newton-Anderson with depth $m=1$, every odd iterate is
simply a Newton step, rather than a combination of the previous two Newton
steps. Hence the algorithm is “restarted” every other step. Explicitly, we
have for $k\geq 1$,

$\displaystyle x_{2k-1}$ $\displaystyle=x_{2k-1}^{\operatorname{Newt}}$
$\displaystyle x_{2k}$
$\displaystyle=x_{2k}^{\operatorname{Newt}}-\gamma_{2k}(x_{2k}^{\operatorname{Newt}}-x_{2k-1}^{\operatorname{Newt}})$

In other words, $\gamma_{2k-1}=0$ for all $k\geq 1$. With
$\operatorname{\gamma\text{NAA}(\hat{r})}$, $\gamma_{k+1}$ is not necessarily
set to zero, but it is scaled towards zero, significantly so depending on
$\beta_{k+1}$. In this way, one may think of
$\operatorname{\gamma\text{NAA}(\hat{r})}$ as a quasi-restarted Anderson
scheme when the depth $m=1$. At the moment, $\gamma$-safeguarding is not
developed for $m>1$, but this interpretation of
$\operatorname{\gamma\text{NAA}(\hat{r})}$ as a quasi-restarted method could
lead to such a development. Presently, these are only heuristics, but they
provide interesting questions for future projects.

The third technique we demonstrate to solve these problems near bifurcation
points is increasing the depth $m$. Evidently, the right choice of $m$ can
significantly improve convergence by reducing the number of iterations to
convergence by half, and increase the domain of convergence with respect to
the parameter. We found, however, that such performance was very sensitive to
the choice of $m$. So while the results suggest this could be developed into a
viable strategy, more work is required to achieve this.

###  4.3 Asymptotic Safeguarding

Under the assumptions of 3.1 , $\operatorname{\gamma\text{NAA}(\hat{r})}$ is
guaranteed to converge locally. This motivates the strategy of this
subsection. As discussed above, asymptotic safeguarding is when we only apply
$\operatorname{\gamma\text{NAA}(\hat{r})}$ once the residual is smaller than a
set threshold. This allows one to take full advantage of NA in the
preasymptotic regime, and ensures fast quadratic convergence for nonsingular
problems. In practice, this means fast local convergence provided that NA
reaches the domain of convergence. We activate
$\operatorname{\gamma\text{NAA}(\hat{r})}$ when $\|w_{k+1}\|<10^{-1}$ in the
experiments below. When $\|w_{k+1}\|>10^{-1}$, we run NA without safeguarding.

####  4.3.1 Results for Channel Flow Model

The results when applied to model (28) are shown below in figures 4, 5, and 6.
The takeaway is that when applied asymptotically, convergence is
$\operatorname{\gamma\text{NAA}(\hat{r})}$ is not as sensitive to the choice
of $\hat{r}$ as it is when applied in the preasymptotic regime (see section
4.4), and adaptive $\gamma$-safeguarding is working as intended by detecting
that the problem is nonsingular. This is seen in plot on the right of figures
4, 5, and 6. Recall that $r_{k}$ is the adaptive parameter in algorithm 5 that
determines how significantly an NA step is scaled to a Newton step. When
$r_{k}\approx 0$, and the criteria in algorithm 5, the NA step $x_{k+1}$ will
be heavily scaled towards a Newton step. When $r_{k}\approx 1$, this scaling
will be much less severe. From our discussion in section 3, we want $r_{k}\to
0$ when the problem is nonsingular in order to enjoy local quadratic
convergence. This is not guaranteed with NA [30]. Observing the $r_{k}$ plots
in figures 4, 5, and 6, one notes that $r_{k}\to 0$ as the solver converges.
Since $r_{k}=\min\\{\eta_{k+1},\hat{r}\\}$, this is equivalent to
$\eta_{k+1}\to 0$ as the solve converges. Thus by theorem 3.2, the
$\operatorname{\gamma\text{NAA}(\hat{r})}$ iterates converge to
$\operatorname{NA}$ iterates asymptotically. This is precisely what
$\operatorname{\gamma\text{NAA}(\hat{r})}$ was designed to do: detect
nonsingular problems, and respond by scaling the iterates towards
$\operatorname{NA}$ asymptotically. Since algorithm 5 is only activated when
$\|w_{k+1}\|<10^{-1}$ in these examples of asymptotic safeguarding, $r_{k}$ is
only computed, and plotted, for the last two iterations. In the next section
on preasymptotic safeguarding, a more interesting $r_{k}$ history is seen. Our
methods, including $\operatorname{NA}$, failed to converge for $\mu<0.92$. If
one wishes to solve a problem at a particular parameter, and a direct solve
fails like we see here for $\mu<0.92$, one could still employ
$\operatorname{NA}$ or $\operatorname{\gamma\text{NAA}(\hat{r})}$ to solve the
problem directly for a parameter value close to the desired one to obtain an
initial guess for continuation. The benefit here is that the continuation
would be required in a smaller parameter range, thus reducing the total number
of solves required. We will see in section 4.5 that increasing $m$ can lead to
convergence for a wider range of parameters.

Figure 4:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
asymptotically with Newton and NA applied to (28) with $\mu=0.96$.

Figure 5:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
asymptotically with Newton and NA applied to (28) with $\mu=0.94$.

Figure 6:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
asymptotically with Newton and NA applied to (28) with $\mu=0.92$.

####  4.3.2 Results for Rayleigh-Bénard Model

The results of asymptotic safeguarding with activation threshold $0.1$ applied
to (29) are similar to those of (28) seen in the previous section. In this
case, however, NA diverged for $\operatorname{Ri}=3.0$ and
$\operatorname{Ri}=3.2$. For $\operatorname{Ri}=3.0$, Newton’s method
converges, but all methods diverged for $\operatorname{Ri}=3.2$. In section
4.4, we are able to recover convergence for $\operatorname{Ri}=3.2$ with
preasymptotic safeguarding. For $\operatorname{Ri}=3.1$,
$\operatorname{Ri}=3.3$, $\operatorname{Ri}=3.4$, and $\operatorname{Ri}=3.5$,
Newton diverged while NA and $\operatorname{\gamma\text{NAA}(\hat{r})}$
converged. Like with (28), we see $r_{k}\to 0$ as the solve converges. The
results for $\operatorname{Ri}=3.4$, shown below, are representative of the
others for which NA converged.

Figure 7:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
asymptotically with Newton and NA applied to (29) with
$\operatorname{Ri}=3.4$.

###  4.4 Preasymptotic Safeguarding

In this section, we demonstrate the preasymptotic safeguarding strategy, where
$\operatorname{\gamma\text{NAA}(\hat{r})}$ is activated starting with $k=2$,
the first iteration where $\operatorname{NA}$ can be applied. Compared to
asymptotic safeguarding, preasymptotic safeguarding is more sensitive to the
choice of $\hat{r}$, but, with the right choice of $\hat{r}$, it can recover
convergence when both Newton and $\operatorname{NA}$ fail.

####  4.4.1 Results for Channel Flow Model

Here we present the results from applying the preasymptotic safeguarding
strategy to the channel flow model (28). Evidently, Newton’s method can still
converge quickly with the right initial guess when $\mu=0.96$. Figure 8
demonstrates that in this case, $\operatorname{NA}$ and
$\operatorname{\gamma\text{NAA}(\hat{r})}$ also perform well. The right-most
plot in figure 8 demonstrates that the
$\operatorname{\gamma\text{NAA}(\hat{r})}$ is working as intended. That is,
$r_{k}\to 0$, and therefore $\operatorname{\gamma\text{NAA}(\hat{r})}$ is
detecting that Newton is converging quickly, and responds by scaling NA
towards a pure Newton iteration. With $\mu=0.94$, we are closer to the
bifurcation point, and observing the left plot in figure 9, we see that Newton
takes many more iterations to converge. However, after a long preasymptotic
phase, Newton does eventually converge quickly, which suggests that the
problem is not truly singular for $\mu=0.94$. This is again detected by
$\operatorname{\gamma\text{NAA}(\hat{r})}$. Even though the
$\operatorname{\gamma\text{NAA}(\hat{r})}$ algorithms take a few more
iterations to converge than NA with $\mu=0.94$, the terminal order of
convergence of $\operatorname{\gamma\text{NAA}(\hat{r})}$ is greater than that
of NA. Approximating the rate by $q_{k+1}=\log(\|w_{k+1}\|)/\log(\|w_{k}\|)$
at each step $k$, and letting $q_{\text{term}}$ denote the terminal order, we
found that $q_{\text{term}}=1.537$ for NA, $q_{\text{term}}=3.386$ for
$\gamma\text{NAA}(0.1)$, and $q_{\text{term}}=2.091$ for
$\gamma\text{NAA}(0.9)$. In the right most plot in figure 9, we take as our
initial iterate the zero vector, but with the fifth entry set to 50. From this
perturbed initial guess, Newton’s method is seen to perform better than NA.
Further, $\gamma\text{NAA}(0.1)$ and $\gamma\text{NAA}(0.5)$ outperform both
Newton and NA, with $\gamma\text{NAA}(0.5)$ converging in about half as many
iterations as NA. This demonstrates the flexibility offered by
$\operatorname{\gamma\text{NAA}(\hat{r})}$. That is, whether Newton or NA is
the best choice for a particular problem and initial guess,
$\operatorname{\gamma\text{NAA}(\hat{r})}$ is more agnostic to these choices,
and can perform well in either case. There is still, however, sensitivity to
$\hat{r}$. We also tested $\operatorname{\gamma\text{NAA}(\hat{r})}$ with
other choices of $\hat{r}$. The method failed to converge for $\hat{r}=0.3$,
$0.5$, $0.6$, $0.7$, and $0.8$. $\operatorname{\gamma\text{NAA}(\hat{r})}$
converged with $\hat{r}=0.4$, but only after 74 iterations. We observed that
$\gamma$NA(0.5), the non-adaptive version of $\gamma$-safeguarding, managed to
converge. The reason for this variation is likely due to the complex behavior
of Newton and NA in the preasymptotic regime, leading to sensitivity to the
choice of $\hat{r}$. With preasymptotic safeguarding, our methods failed for
$\mu<0.94$, hence continuation may still be required in some cases, but NA and
$\operatorname{\gamma\text{NAA}(\hat{r})}$ can be used to efficiently solve
the problem closer to the bifurcation point compared to Newton, thereby
reducing the overall computational cost. The point is that applying
$\gamma$-safeguarding and adaptive $\gamma$-safeguarding in the preasymptotic
regime can still converge faster than Newton, but this convergence is again
sensistive to the choice of $\hat{r}$, and further work is required to
understand this sensitivity.

Figure 8:  Left: Convergence history for Newton, Newton-Anderson, and adaptive
$\gamma$-Newton-Anderson for $\mu=0.96$. Right: The value of $r_{k}$ at each
iteration for the adaptive methods.

Figure 9:  Left: Convergence history for Newton, Newton-Anderson, and adaptive
$\gamma$-Newton-Anderson for $\mu=0.94$ with zero vector initial guess. Right:
Convergence history for Newton, Newton-Anderson, and adaptive $\gamma$-Newton-
Anderson for $\mu=0.94$ with perturbed initial guess.

####  4.4.2 Results for Rayleigh-Bénard Model

With preasymptotic safeguarding applied to model (29), we again see more
varied behavior since we have $\gamma$-safeguarding activated from the
beginning of the solve. The results are shown in figures 10, 11, and 12. The
$r_{k}$ plots are only shown for those values of $\hat{r}$ for which
$\operatorname{\gamma\text{NAA}(\hat{r})}$ converged. We found that for
$\operatorname{Ri}=3.0$ and $\operatorname{Ri}=3.1$,
$\operatorname{\gamma\text{NAA}(\hat{r})}$ failed to converge with
$\hat{r}=0.1$, $0.5$, or $0.9$. Convergence is recovered with $\hat{r}=0.4$
and $\hat{r}=0.6$ respectively. We also ran
$\operatorname{\gamma\text{NAA}(\hat{r})}$ with these $\hat{r}$ values for
$\operatorname{Ri}=3.2$, $3.3$, $3.4$, and $3.5$. The results were similar for
$\operatorname{Ri}=3.3$, $3.4$, and $3.5$. Hence we only show the results for
$\operatorname{Ri}=3.4$ in figure 12. The theme demonstrated in figures 10,
11, and 12 is that when preasymptotic safeguarding is employed, it is possible
to converge faster than standard NA. Moreover, as seen in figure 11,
$\operatorname{\gamma\text{NAA}(\hat{r})}$ can converge when both Newton and
NA diverge. One may note that in the right-most plot in figure 10, prior to
the asymptotic regime where $r_{k}\to 0$, we observe $r_{k}<\hat{r}=0.6$ only
twice. Using the restarted Anderson interpretation discussed in section 4.2,
we could say that there are two quasi-restarts prior to the asymptotic regime.
Evidently, these two quasi-restarts are essential, since we found that non-
adaptive $\gamma$-safeguarding with $\hat{r}=0.6$, algorithm 4, diverges.
Similar behavior of $r_{k}$ is seen in figures 11 and 12. It remains unclear
precisely how the choice of $\hat{r}$ affects convergence, e.g., in figure 10,
why does $\gamma\text{NAA}(0.6)$ converge, but
$\operatorname{\gamma\text{NAA}(\hat{r})}$ diverges for $\hat{r}=0.1$, $0.4$,
$0.5$, and $0.9$? As previously discussed, a better understanding of these
methods in the preasymptotic regime would help answer questions like these,
and this will be the focus of future studies.

Figure 10:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
preasymptotically with Newton and NA applied to (29) with
$\operatorname{Ri}=3.1$.

Figure 11:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
preasymptotically with Newton and NA applied to (29) with
$\operatorname{Ri}=3.2$.

Figure 12:  Comparison of $\operatorname{\gamma\text{NAA}(\hat{r})}$ applied
preasymptotically with Newton and NA applied to (29) with
$\operatorname{Ri}=3.4$.

###  4.5 Increasing Anderson Depth

The strategy employed in this section is to increase the depth of the NA from
$m=1$ used in previous sections.

####  4.5.1 Results for Channel Flow Model

The plots in figure 13 demonstrate the effectiveness of increasing the
algorithmic depth $m$ to solve (28) near the bifurcation point. We also
experimented with applying $\operatorname{\gamma\text{NAA}(\hat{r})}$
asymptotically. When $\|w_{k+1}\|<1$, we set $m=1$ and activated
$\operatorname{\gamma\text{NAA}(\hat{r})}$ with $\hat{r}=0.9$. These results
are seen as the dashed lines in the left-most plot in figure 13. The
philosophy is similar to that of preasymptotic safeguarding. We use NA(3) to
reach the asymptotic regime, and then allow adaptive $\gamma$-safeguarding to
detect if the problem is nonsingular. We set $m=1$ because, for the present,
$\gamma$-safeguarding is only designed for $m=1$. The left-most plot in figure
13 shows how, from the same initial iterate, we are able to solve (28) for a
wider range of $\mu$ values, including in the regime where Newton,
$\operatorname{NA}(1)$, and $\operatorname{\gamma\text{NAA}(\hat{r})}$ failed
to converge.

The right-most plot in figure 13 focuses on the results of applying
$\operatorname{NA}(3)$ to (28) with $\mu=0.92$. The point here is that even in
the regime where $\operatorname{NA}(1)$ converges, $\operatorname{NA}(3)$
converges in about half as many iterations. Thus increasing the algorithmic
depth can lead to faster convergence. The improved convergence seen with
increasing the depth $m$ suggests that a generalization of
$\gamma$-safeguarding for greater depths could be useful as a generalization
of the strategies presented in previous sections. However, for this strategy
to be effective in general, further study is needed on the proper choice of
depth $m$. In the chosen parameter regime, $m=3$ was the only value of
$m\in\\{1,2,...,10\\}$ observed to improve convergence when NA(1) failed.

Figure 13:  Left: Convergence history of NA(3) for various choices of $\mu$
near $\mu^{*}$. Solid lines are NA$(3)$, and dashed lines are NA$(3)$ with
asymptotic $\operatorname{\gamma\text{NAA}(\hat{r})}$. Right: Convergence
history for NA(1) and NA(3) with $\mu=0.92$.

####  4.5.2 Results for Rayleigh-Beńard Model

The results of increasing $m$ to solve (29) are shown below in figure 14. As
with (28), we experimented with reducing $m=1$ and activating
$\gamma$-safeguarding asymptotically. The dashed lines in 14 are the results
of these experiments. For this problem, activation occurred when
$\|w_{k+1}\|<10^{-1}$ since we found that activating
$\operatorname{\gamma\text{NAA}(\hat{r})}$ with $\hat{r}=0.9$ when
$\|w_{k+1}\|<1$ broke convergence like it did in section 4.3.2. We once again
observe that the right choice of $m>1$ can improve convergence significantly.
This can be seen in figure 14 for $\operatorname{Ri}=3.3$. Moreover, figure 14
demonstrates that increasing $m$ can recover convergence when
$\operatorname{NA}(1)$ fails to converge for $\operatorname{Ri}=3.0$ and
$\operatorname{Ri}=3.2$. However, for $\operatorname{Ri}=$ 3.1, 3.4, and 3.5,
there was no significant improvement gained from increasing $m>1$. The results
for $\operatorname{Ri}=3.4$ shown below in figure 14 are representative of the
results for $\operatorname{Ri}=3.1$ and $\operatorname{Ri}=3.5$.

The most significant difference between $\operatorname{NA}(m)$ without
asymptotic $\operatorname{\gamma\text{NAA}(\hat{r})}$, and
$\operatorname{NA}(m)$ with asymptotic
$\operatorname{\gamma\text{NAA}(\hat{r})}$, is seen with $m=10$. It also
appears that $m=10$ is more sensitive to the activation threshold than smaller
choices of $m$. This is seen in figure 14 for $\operatorname{Ri}=3.3$ and
$\operatorname{Ri}=3.4$. These results are promising, and motivate further
investigation to fully understand how the choice of $m$ affects convergence
near singularities, and in particular, near bifurcation points.

Figure 14:  Convergence history of $\operatorname{NA}(m)$ for $m=$ 1,2,3,4,5,
and 10 for $\operatorname{Ri}=$3.0, 3.2, 3.3, and 3.4. Dashed lines denote
$\operatorname{NA}(m)$ with asymptotic
$\operatorname{\gamma\text{NAA}(\hat{r})}$ with $\hat{r}=0.9$ activated when
$\|w_{k+1}\|<10^{-1}$.

##  5 Conclusion

We have presented a modification of Anderson accelerated Newton’s method for
solving nonlinear equations near bifurcation points. We proved that, locally,
this modified scheme can detect nonsingular problems and scale the iterates
towards a pure Newton step, which leads to faster local convergence compared
to standard Newton-Anderson. We numerically demonstrated two strategies one
can employ when using our modified Newton-Anderson scheme to solve nonlinear
problems near bifurcation points, with our test problems being two Navier-
Stokes type parameter-dependent PDEs. Asymptotic safeguarding was shown to
recover local quadratic convergence when the problem is nonsingular, and it
shows virtually no sensitivity to the choice of parameter $\hat{r}$. It can,
however, be sensitive to the choice of activation threshold. Preasymptotic
safeguarding is shown to significantly improve convergence, and can recover
convergence when both Newton and $\operatorname{NA}$ fail. There is strong
sensitivity to the choice of $\hat{r}$ though, and future work will clarify
this dependence. We also demonstrated that increasing the Anderson depth $m$
can improve convergence, and increase the domain of convergence with respect
to the problem parameter. Future projects will further study how the choice of
$m$ impacts convergence, and work towards developing $\gamma$-safeguarding for
greater algorithmic depths.

##  6 Acknowledgements

MD and SP are supported in part by the National Science Foundation under Grant
No. DMS-2011519. LR is supported in part by the National Science Foundation
under Grant No. DMS-2011490. This material is based upon work supported by the
National Science Foundation under Grant No. DMS-1929284 while the authors were
in residence at the Institute for Computational and Experimental Research in
Mathematics in Providence, RI, during the Acceleration and Extrapolation
Methods (MD, SP and LR), and the Numerical PDEs: Analysis, Algorithms and Data
Challenges (SP and LR), programs.

## References

  * [1] D.G. Anderson.  Iterative procedures for nonlinear integral equations.  J. Assoc. Comput. Mach., 12(4):547–560, 1965.
  * [2] R. Behling and A. Fishcer.  A unified local convergence analysis of inexact constrained Levenberg-Marquardt methods.  Optimization Letters, 6:927–940, 2012.
  * [3] S. Bellavia and B. Morini.  Strong local convergence properties of adaptive regularized methods for nonlinear least squares.  IMA Journal of Numerical Analysis, 35, 2015.
  * [4] N. Boullé, V. Dallas, and P.E. Farrell.  Bifurcation analysis of two-dimensional Rayleigh-Bénard convection using deflation.  Phys. Rev. E, 105, 2022.
  * [5] D. Braess.  Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics.  Cambridge University Press, 2007.
  * [6] M. Chupin, M. Dupuy, G. Legendre, and E. Séré.  Convergence analysis of adaptive DIIS algorithms with application to electronic ground state calculations.  ESAIM: Mathematical Modelling and Numerical Analysis, 55(6):2785–2825, 2021.
  * [7] M. Dallas and S. Pollock.  Newton-Anderson at singular points.  Int. J. Numer. Anal. and Mod., 20(5):667–692, 2023.
  * [8] D.W. Decker, H.B. Keller, and C. T. Kelley.  Convergence rates for Newton’s method at singular points.  SIAM J. Numer. Anal., 20(2):296–314, 1983.
  * [9] D.W. Decker and C.T. Kelley.  Newton’s method at singular points I.  SIAM J. Numer. Anal., 17(1):66–70, 1980.
  * [10] D.W. Decker and C.T. Kelley.  Convergence acceleration for Newton’s method at singular points.  SIAM J. Numer. Anal., 19(1):219–229, 1982.
  * [11] D.W. Decker and C.T. Kelley.  Expanded convergence domains for Newton’s method at nearly singular roots.  SIAM J. Sci. Stat. Comput., 6(4), 1985.
  * [12] P. Deuflhard.  Newton Methods for Nonlinear Problems.  Springer Series in Computational Mathematics. Springer, 2005.
  * [13] P.E. Farrell, A. Birkisson, and S.W. Funke.  Deflation techniques for finding distinct solutions of nonlinear partial differential equations.  SIAM J. Sci. Comput., 37, 2015.
  * [14] J.K. Galvin, A. Linke, L. Rebholz, and N. Wilson.  Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection.  Comput. Methods Appl. Mech. Engrg., pages 166–176, 2012.
  * [15] A.O. Griewank.  Analysis and Modification of Newton’s Method at Singularities.  PhD thesis, 1980.
  * [16] A.O. Griewank.  Starlike domains of convergence for Newton’s method at singularities.  Numer. Math., 35:95–111, 1980.
  * [17] J. Guzḿan and L.R. Scott.  The Scott-Vogelius finite elements revisited.  Mathematics of Computation, 88(316):515–529, 2019.
  * [18] H. He, S. Zhao, Y. Xi, Joyce C. Ho, and Yousef Saad.  Solve minimax optimization by Anderson acceleration.  International Conference on Learning Representations, 2022.
  * [19] N.C. Henderson and R. Varadhan.  Damped Anderson acceleration with restarts and monotonicity control for accelerating EM and EM-like algorithms.  Journal of computational and graphical statistics, 28(4):834–846, 2019.
  * [20] J.L. Hueso, E.M., and J.R. Torregrosa.  Modified Newton’s method for systems of nonlinear equations with singular Jacobian.  J. Comput. Appl. Math., 224(1):77–83, 2008.
  * [21] C. Kanzow, N. Yamashita, and M. Fukushima.  Levenberg–Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints.  J. Comput. Appl. Math, 172:375–397, 2004.
  * [22] C.T. Kelley and R. Suresh.  A new acceleration method for Newton’s method at singular points.  SIAM J. Numer. Anal, 20(5):1001–1009, 1983.
  * [23] H. Kielhöfer.  Bifurcation Theory: An Introduction with Applications to Partial Differential Equations.  Springer Science+Business Media, 2012.
  * [24] M.L. Pasini.  Convergence analysis of Anderson-type acceleration of Richardson’s iteration.  Numer. Linear Algebra Appl., 26, 2019.
  * [25] F. Pichi.  Reduced order models for parametric bifurcation problems in nonlinear PDEs.  PhD thesis, 2020.
  * [26] S. Pollock and L.G. Rebholz.  Anderson acceleration for contractive and noncontractive operators.  IMA J. Numer. Anal., 41(4):2841–2872, 2021.
  * [27] S. Pollock and L.G. Rebholz.  Filtering for Anderson acceleration.  SIAM J. Sci. Comput., 45(4):A1571–A1590, 2023.
  * [28] S. Pollock and H. Schwartz.  Benchmarking results for the Newton–Anderson method.  Results Appl. Math., 8:100095, 2020.
  * [29] J. Qin.  On the convergence of some low order mixed finite elements for incompressible fluids.  PhD thesis, 1994.
  * [30] L. Rebholz and M. Xiao.  The effect of Anderson acceleration on superlinear and sublinear convergence.  Journal of Scientific Computing, 96(2), 2023.
  * [31] G.W. Reddien.  On Newton’s method for singular problems.  SIAM J. Numer. Anal, 15(5):993–996, 1978.
  * [32] G.W. Reddien.  Newton’s method and high order singularities.  Comput. Math. Appl., 5(2):79–86, 1979.
  * [33] R.B. Schnabel and P.D. Frank.  Tensor methods for nonlinear equations.  SIAM J. Numer. Anal., 21(5):815–843, 1984.
  * [34] R.S. Tuminaro, H.F. Walker, and J. N. Shadid.  On backtracking failure in Newton-GMRES methods with demonstration for the Navier-Stokes equations.  Journal of Computational Physics, 180:549–558, 2002.
  * [35] H. Uecker.  Continuation and bifurcation in nonlinear pdes – algorithms, applications, and experiments.  Jahresbericht der Deutschen Mathematiker-Vereinigung, 124:43–80, 2021.
  * [36] H.F. Walker and P. Ni.  Anderson acceleration for fixed-point iterations.  SIAM J. Numer. Anal., 49(4):1715–1735, 2011.

##  Appendix A Algorithms

1: Choose $x_{0}\in\mathbb{R}^{n}$.

2: for k=1,2,… do

3: $w_{k+1}\leftarrow-f^{\prime}(x_{k})^{-1}f(x_{k})$

4: $x_{k+1}\leftarrow x_{k}+w_{k+1}$

5: end for

Algorithm 1 Newton

1: Choose $x_{0}\in\mathbb{R}^{n}$ and $m\geq 0$. Set
$w_{1}=-f^{\prime}(x_{0})^{-1}f(x_{0})$, and $x_{1}=x_{0}+w_{1}$.

2: for k=1,2,… do

3: $m_{k}\leftarrow\min\\{k,m\\}$

4: $w_{k+1}\leftarrow-f^{\prime}(x_{k})^{-1}f(x_{k})$

5: $F_{k}=\big{(}(w_{k+1}-w_{k})\cdots(w_{k-m+2}-w_{k-m+1})\big{)}$

6: $E_{k}=\big{(}(x_{k}-x_{k-1})\cdots(x_{k-m+1}-x_{k-m})\big{)}$

7:
$\gamma_{k+1}\leftarrow\text{argmin}_{\gamma\in\mathbb{R}^{m}}\|w_{k+1}-F_{k}\gamma\|_{2}^{2}$

8: $x_{k+1}\leftarrow x_{k}+w_{k+1}-(E_{k}+F_{k})\gamma_{k+1}$

9: end for

Algorithm 2 Newton-Anderson(m)

1: Choose $x_{0}\in\mathbb{R}^{n}$. Set
$w_{1}=-f^{\prime}(x_{0})^{-1}f(x_{0})$, and $x_{1}=x_{0}+w_{1}$.

2: for k=1,2,… do

3: $w_{k+1}\leftarrow-f^{\prime}(x_{k})^{-1}f(x_{k})$

4:
$\gamma_{k+1}\leftarrow(w_{k+1}-w_{k})^{T}w_{k+1}/\|w_{k+1}-w_{k}\|_{2}^{2}$

5: $x_{k+1}\leftarrow x_{k}+w_{k+1}-\gamma_{k+1}(x_{k}-x_{k-1}+w_{k+1}-w_{k})$

6: end for

Algorithm 3 Newton-Anderson(1) Algorithm 4 $\gamma$-safeguarding

1: Given $x_{k}$, $x_{k-1}$, $w_{k+1}$, $w_{k}$, and $\gamma_{k+1}$. Set
$r\in(0,1)$ and $\lambda=1$.

2: $\beta\leftarrow r\|w_{k+1}\|/\|w_{k}\|$

3: if $\gamma_{k+1}=0$ or  $\gamma_{k+1}\geq 1$ then

4: $x_{k+1}\leftarrow x_{k}+w_{k+1}$

5: else

6: if $|\gamma_{k+1}|/|1-\gamma_{k+1}|>\beta$ then

7: if $\gamma_{k+1}>0\textbf{ and }\beta/(\gamma_{k+1}(1+\beta))<1$ then

8: $\lambda\leftarrow\beta/(\gamma_{k+1}(1+\beta))$

9: end if

10: if $\gamma_{k+1}<0\textbf{ and }0\leq\beta/(\gamma_{k+1}(\beta-1))<1$ then

11: $\lambda\leftarrow\beta/(\gamma_{k+1}(\beta-1))$

12: end if

13: end if

14: $\gamma_{k+1}\leftarrow\lambda\gamma_{k+1}$

15: $x_{k+1}\leftarrow
x_{k}+w_{k+1}-\gamma_{k+1}(x_{k}-x_{k-1}+w_{k+1}-w_{k})$

16: end if

For a discussion of notation and theory for algorithm 5, see section 3.

Algorithm 5 Adaptive $\gamma$-safeguarding

1: Given $x_{k}$, $x_{k-1}$, $w_{k+1}$, $w_{k}$, and $\gamma_{k+1}$, set
$\hat{r}\in(0,1)$, $\eta_{k+1}=\|w_{k+1}\|/\|w_{k}\|$, and
$r_{k}=\min\\{\eta_{k+1},\hat{r}\\}$, and $\lambda^{a}=1$.

2: $\beta\leftarrow r_{k}\|w_{k+1}\|/\|w_{k}\|$

3: if $\gamma_{k+1}=0$ or  $\gamma_{k+1}\geq 1$ then

4: $x_{k+1}\leftarrow x_{k}+w_{k+1}$

5: else

6: if $|\gamma_{k+1}|/|1-\gamma_{k+1}|>\beta$ then

7: if $\gamma_{k+1}>0\textbf{ and }\beta/(\gamma_{k+1}(1+\beta))<1$ then

8: $\lambda^{a}\leftarrow\beta/(\gamma_{k+1}(1+\beta))$

9: end if

10: if $\gamma_{k+1}<0\textbf{ and }0\leq\beta/(\gamma_{k+1}(\beta-1))<1$ then

11: $\lambda^{a}\leftarrow\beta/(\gamma_{k+1}(\beta-1))$

12: end if

13: end if

14: $\gamma_{k+1}\leftarrow\lambda^{a}\gamma_{k+1}$

15: $x_{k+1}\leftarrow
x_{k}+w_{k+1}-\gamma_{k+1}(x_{k}-x_{k-1}+w_{k+1}-w_{k})$

16: end if