
\section{Introduction}
Efficient fixed-point schemes are needed in many complex physical applications such as (i) iterative algorithms to solve sparse linear systems, (ii) Picard iterations \cite{Kewang2022} to solve systems of non-linear partial differential equations as in computational fluid dynamics \cite{rebholz, 10.1093/imanum/draa095}, semiconductor modelling \cite{doi:10.1137/20M1365922, LAIU2020109567}, and astrophysics \cite{Paul_Laiu_2020, Paul_Laiu_2021}, and (iii) solving extended-space full-waveform inversion (FWI) \cite{dx.doi.org/10.1190/geo2021-0409.1}.

For scientific applications, the standard fixed-point scheme often converges slowly either because the fixed-point operator is weakly globally contractive or is not globally contractive. 
The Anderson acceleration (AA) \cite{Anderson} is a multi-secant method \cite{brezinski_shanks_2018, 10.1093/imanum/drab061} that has been widely used either to improve the convergence rate of convergent fixed-point schemes or to restore convergence when the original fixed-point scheme is not convergent \cite{Fang, walker_anderson_2011, potra_characterization_2013, zhang2020, saad2021}. In particular, the convergence of AA has been studied with respect to specific properties of different scientific applications \cite{toth_convergence_2015, toth_local_2017, Mai2019, doi:10.1137/20M132938X}.

AA requires solving a least-squares at each fixed-point iteration, and this can be computationally expensive when scientific applications involve large-scale calculations, especially when the solver is distributed on high-performance computing (HPC) platforms.
To reduce the computational cost of AA, a recently proposed variant of AA called Alternating Anderson Richardson (AAR) \cite{banerjee_periodic_2016, AAJ, AAR} performs multiple fixed-point iterations between two consecutive AA steps. Recent theoretical results have shown that AAR is more robust than the standard AA against stagnations \cite{lupo_pasini_convergence_2019}. 
Even if performed at periodic intervals, computing the Anderson mixing via least-squares can still be computationally expensive when the number of rows in the tall and skinny matrix that defined the left-hand side (LHS) of the least-squares problem is large. 
One possibility to reduce the computational cost consists in projecting the least-squares onto a subspace. However, the projection of the least-squares onto a subspace must be performed judiciously in order not to compromise convergence. 
Although for specific problems it has been shown that projecting the least-squares to perform AA onto an appropriate projection subspace results in computational savings without affecting accuracy \cite{doi:10.1137/19M1290097, doi.org/10.1111/cgf.14081}, physics driven guidelines to identify a projection subspace are difficult (sometimes impossible) to determine because of lack of structure in the fixed-point operator that could guide an appropriate definition of the projection subspace. In this situation, general guidelines (not necessarily physics driven) to identify a projection subspace are needed. 

We provide rigorous theoretical bounds for AA on linear fixed-point iterations that enable to project the least-squares on a subspace and perform the Anderson mixing onto that subspace while still ensuring the final residual to drop below a user defined convergence threshold. 
By interpreting the projection of the least-squares problem onto a subspace as perturbation to the original least-squares, we assess how much accuracy can be sacrificed to limit the communication and computational burden without compromising convergence. On the same lines as
previously published theoretical results in \cite{simoncini} for well-established Krylov methods \cite{Saad2003}, our analysis concludes that also for AAR one can progressively reduce the accuracy with which calculations (e.g., evaluations of the fixed-point operator or solving least-squares for AA) are performed. 

Besides being effective to solve sparse linear systems, the error bounds open a path to efficiently use AAR also for non-linear fixed-point iterations arising also from quantum mechanics, plasma astrophysics, and computational fluid dynamics where function evaluations are also expensive to perform \cite{toth_local_2017}. In this context, the use of AAR allows to leverage less accurate but also less expensive function evaluations (which affect the accuracy of the calculations of the residual) without affecting the final attainable accuracy of the non-linear physics solver. 

Since the theoretical bounds are defined in terms of quantities that are not computationally convenient to retrieve, we also propose heuristics that are inexpensive to compute while still maintaining a close relation to the rigorous theoretical estimate. The heuristics allow to dynamically tune the size of the projection subspace at each iteration without introducing additional computational overhead. 
We combine the heuristics with a check on the monotonic reduction of the residual across two consecutive Anderson mixing steps to ensure that the projection of the least squares onto a subspace does not compromise the convergence of the numerical scheme. 
The monotonicity check allows one to replace the tight theoretical bound with an empirical one, and dynamically adjust (increase) the value of the heuristics if non-monotonicity is observed across two consecutive AA steps. 

Numerical results confirm that approximate calculations to solve the least-squares for AA saves computational time without sacrificing convergence within a desired accuracy. 
The proposed method has appealing properties for high-performance computing since it reduces computational requirements, inter-process communications, and storage requirements to solve least-squares required for AA when the fixed-point operator is distributed across multi-node processes. 

The remainder of the paper is organized as follows. In Section \ref{section:background} we briefly recall the stationary Richardson's method that solves linear systems by recasting them as linear fixed-point iterations. 
In Section \ref{section:error_bounds} we provide error bounds that allow approximate calculations of AA for linear fixed-point iterations while still maintaining convergence.
In Section \ref{section:reduced_aar} we introduce the new formulation of AAR that enables the projection of the least-squares onto a projection subspace, called Reduced AAR, followed in Section \ref{section:numerical_results} by numerical implementations to illustrate the analyses conducted in the previous sections for two different scientific applications that we use as representatives of three different types of fixed-point iterations. The applications we consider are: (i) linear deterministic fixed-point iterations arising from the Richardson's scheme to solve linear systems and (ii) non-linear deterministic fixed-point iterations arising from non-linear time-dependent Boltzmann equations. We conclude the paper with remarks on the state-of-the-art and comments about future developments in Section \ref{section:conclusions}. 

\section{Fixed-point iteration}
The standard iterative method for solving a fixed-point problem
$
\mathbf{x} = G(\mathbf{x})$ with $\mathbf{x}\in\mathbb{R}^n$ and $G:\mathbb{R}^n \rightarrow \mathbb{R}^n$ is the fixed-point iteration:
\begin{equation}
\mathbf{x}^{k+1} = G(\mathbf{x}^{k}), \quad k=0,1,\ldots,
\label{fixed_point}
\end{equation}
and the residual $\mathbf{r}^k \in \mathbb{R}^n$ of the fixed-point iteration is defined as 
\begin{equation}
\mathbf{r}^k = G(\mathbf{x}^k)-\mathbf{x}^{k} , \quad k=0,1,\ldots.
\end{equation}
The convergence of the fixed-point iteration relies on the global contraction property of the nonlinear fixed-point operator $G$. When the contraction constant of $G$ is close to one, the fixed-point iteration may converge at an unacceptably slow rate.

\subsection{Anderson acceleration and the Alternating Anderson acceleration}
Anderson acceleration (AA) was proposed in \cite{Anderson} to accelerate the fixed-point iteration. There exist several equivalent formulations of AA \cite{Fang, toth_convergence_2015}. The formulation we adopt for convenience in the analysis is described in Algorithm \ref{AAR_algorithm}.
Define
\begin{equation}
    X_k = [(\mathbf{x}^{k-\ell+1} - \mathbf{x}^{k-\ell}), \ldots, (\mathbf{x}^k-\mathbf{x}^{k-1}) ]\in \mathbb{R}^{n\times \ell}
\end{equation}
and 
\begin{equation}
    R_k = [(\mathbf{r}^{k-\ell+1} - \mathbf{r}^{k-\ell}), \ldots, (\mathbf{r}^k-\mathbf{r}^{k-1}) ]\in \mathbb{R}^{n\times \ell}
\end{equation}
where $\ell=\min\{k,m\}$ and $m$ is an integer that describes the maximum number of previous terms of the sequence $\{\mathbf{x}^k\}_{k=0}^{\infty}$ used to compute the AA update. Denoting with $\lVert \cdot \rVert_2$ the $\ell^2$-norm of a vector, the vector $\mathbf{g}^k = [g^{(k)}_1, \ldots, g^{(k)}_{\ell}]^T \in \mathbb{R}^\ell$ defined as
\begin{equation}
    \mathbf{g}^k=\underset{\mathbf{g}\in \mathbb{R}^{\ell}}{\operatorname{argmin}}\lVert R_k\mathbf{g} - \mathbf{r}^k\rVert_2^2
   \label{least_squares_matrix_form}
\end{equation}
is used to update the sequence $\{\mathbf{x}^k\}_{k=0}^{\infty}$ through an Anderson mixing as follows
\begin{equation}
    \mathbf{x}^{k+1} = \mathbf{x}^k - (X_k+R_k) \mathbf{g}^k.
\end{equation}
A relaxation parameter $0\le \omega \le 1$ can be used to smooth the AA update by performing a weighted average between a standard Picard iteration and AA 
\begin{equation}
    \mathbf{x}^{k+1} = (1-\omega)\mathbf{x}^k + \omega [ \mathbf{r}^k-(X_k+ R_k)\mathbf{g}^k].
\end{equation}

\begin{algorithm}
\caption{(Alternating) Anderson Acceleration}
\label{AAR_algorithm}
\begin{algorithmic}
\STATE{Data: $\mathbf{x}^0\in\mathbb{R}^n$, $\omega\in[0,1]$, $p\in\mathbb{N}$; \hfill \COMMENT{$// \texttt{AA: } p=1,\,\,\texttt{Alternating AA: } p>1$}}
\STATE{Compute $\mathbf{r}^0 = G(\mathbf{x}^0)-\mathbf{x}^0$ and $\mathbf{x}^1=\mathbf{x}^0+\mathbf{r}^0$;}
\STATE{Set $k=1$}
\WHILE{$\displaystyle \lVert \mathbf{r}^{k-1} \rVert_2>tol$}
 \STATE{Compute $\ell=\min\{k,m\}$}
 \STATE{Compute $\mathbf{r}^k=G(\mathbf{x}^k)-\mathbf{x}^k$}
 \STATE{Set $X_k = [(\mathbf{x}^{k-\ell+1} - \mathbf{x}^{k-\ell}), \ldots, (\mathbf{x}^k-\mathbf{x}^{k-1}) ]\in \mathbb{R}^{n\times \ell}$}
 \STATE{Set $R_k = [(\mathbf{r}^{k-\ell+1} - \mathbf{r}^{k-\ell}), \ldots, (\mathbf{r}^k-\mathbf{r}^{k-1}) ]\in \mathbb{R}^{n\times \ell}$}
 \IF{$k\pmod p\ne 0$}
 \STATE{Update $\mathbf{x}^{k+1}=\mathbf{x}^k + \omega\mathbf{r}^k$}
 \ELSE
 \STATE{Compute  
 $\displaystyle
 \mathbf{g}^k=\underset{\mathbf{g}\in \mathbb{R}^{\ell}}{\operatorname{argmin}}\lVert R_k\mathbf{g} - \mathbf{r}^k\rVert_2^2
 $}
 \STATE{Update $\mathbf{x}^{k+1} = (1-\omega)\mathbf{x}^k + \omega [ \mathbf{r}^k-(X_k+ R_k)\mathbf{g}^k]$}
 \ENDIF
 \STATE{$k=k+1$}
\ENDWHILE
 \RETURN{$\mathbf{x}^{k+1}$ }
\end{algorithmic}
\end{algorithm}

When compared to the fixed-point iteration, AA often requires fewer iterations to converge thus resulting in less computation time. On the other hand, AA also introduces the overhead of solving a least-squares problem.
This computation overhead is outweighed by the benefit from fewer iterations when solving problems in which evaluating the operator $G$ incurs the dominant cost. However, there are many problems where the cost of solving the least-squares problems incurs the main cost.

Solving a least-squares problem at each iteration is computationally expensive, and requires global communications which introduce severe bottlenecks for the parallelization in HPC environments. 
P.~Suryanarayana and collaborators recently proposed to compute
an Anderson mixing after multiple Picard iterations, so that the cost of solving successive least-squares problems is reduced. This new AA variant, called \textit{Alternating AA}, has been shown to effectively accelerate both linear and non-linear fixed-point iterations \cite{AAJ, AAR}. Algorithm \ref{AAR_algorithm} describes Alternating AA, where the parameter $p$ represents the number of Picard iterations separating two consecutive AA steps. Letting $p\rightarrow \infty$ reduces Alternating AA to a standard fixed-point scheme with $\omega$ the relaxation parameter, while $p=1$ makes it coincide with the standard AA formulation. 


In following sections \ref{section:background} and \ref{AA_linear} we focus on the case when $G$ is linear and arises from iteratively solving a sparse linear system. 

\subsection{Stationary Richardson}
\label{section:background}
Consider a nonsingular sparse linear system 
\begin{equation}
A\mathbf{x}=\mathbf{b},
\label{linsys}
\end{equation}
where $A\in\mathbb{R}^{n\times n}$ and $\mathbf{x}$, $\mathbf{b} \in \mathbb{R}^n$. 
We assume that left preconditioning has already been incorporated.
Let $H := I-A$, adding $H\mathbf{x}$ on both sides of \eqref{linsys} leads to a linear fixed-point problem with the fixed-point iteration
\begin{equation}
\mathbf{x}^{k+1} = G(\mathbf{x}^k) := H\mathbf{x}^k + \mathbf{b},\quad k=0,1,\ldots.
\label{linear_fixed_iteration}
\end{equation}
A commonly used equivalent representation of \eqref{linear_fixed_iteration} is the correction form
\begin{equation}
\mathbf{x}^{k+1} =\mathbf{x}^k + \mathbf{r}^k,\quad k=0,1,\ldots,
\label{fixedpoint}
\end{equation}
where $\mathbf{r}^k = \mathbf{b}-A\mathbf{x}^k$ is the residual at the $k$th iteration. 
The scheme \eqref{fixedpoint} converges to the solution if and only if the spectral radius $\rho(H)<1$. 
To ensure that this condition on the spectral radius is respected, the scheme \eqref{linear_fixed_iteration} can be generalized by introducing a positive weighing parameter $\omega$:  
\begin{equation}
\mathbf{x}^{k+1}=(1-\omega)\mathbf{x}^k + \omega( H\mathbf{\mathbf{x}}^k+\mathbf{b}),
\label{weighted}
\end{equation}
which in correction form is equivalently represented as
\begin{equation}
\mathbf{x}^{k+1} = \mathbf{x}^k + \omega \mathbf{r}^{k}.
\label{weighted2}
\end{equation}
Equations \eqref{weighted} and \eqref{weighted2} are known as the \textit{stationary Richardson scheme}. 


\subsection{Anderson acceleration for linear systems}
\label{AA_linear}
AA is well suited to potentially accelerate Richardson schemes. 
In the case of a linear fixed-point iteration as in Equation \eqref{linear_fixed_iteration}, AA can be split in two steps. The first step consists in calculating the Anderson mixing
\begin{equation}
\overline{\mathbf{x}}^k = \mathbf{x}^k - X_k\mathbf{g}^k,
\label{anderson_correction}
\end{equation}
where the vector $\mathbf{g}^k$ is computed by solving the over-determined least-squares in Equation \eqref{least_squares_matrix_form}.
Under the assumption that $R_k$ is full rank%
\footnote{When $R_{k}$ is not full rank, using $R_\ell$ to compute an 
Anderson mixing as in Equation \eqref{anderson_correction} is 
still legitimate by solving least-squares problems. Indeed, albeit some numerical examples presented in Section \ref{section:numerical_results} generated 
a matrix $R_\ell$ with linearly dependent columns at some iterations, this did not hinder convergence.}, 
we have
$ \mathbf{g}^k  = (R^T_kR_k)^{-1}R^T_k\mathbf{r}^{k}$.
Since $R_k = -AX_k$ for a linear fixed-point iteration, \eqref{anderson_correction} can be recast as 
\begin{equation}
\overline{\mathbf{x}}^k 
= \mathbf{x}^k - X_k(R^T_kR_k)^{-1}R^T_k\mathbf{r}^{k} 
= \mathbf{x}^k + X_k(R^T_k A X_k)^{-1}R^T_k\mathbf{r}^{k}.
\end{equation}
Therefore, the Anderson mixing can be interpreted as a step of an \textit{oblique projection method} (see Chapter 5 of \cite{Saad}) 
with the subspace of corrections chosen as
$
\mathcal{V}_k = \mathcal{R}(X_k)
$
and the subspace of constraints as
$
\mathcal{W}_k = \mathcal{R}(R_k) = A \mathcal{R}(X_k) = A \mathcal{V}_k
$.
After the Anderson mixing is applied, $\overline{\mathbf{x}}^k$ is used to compute a new update via a standard Richardson's step
\begin{equation}
 \mathbf{x}^{k+1} = \overline{\mathbf{x}}^k + \omega \overline{\mathbf{r}}^k.
 \label{relaxation_step}
\end{equation}
If AA is applied at each fixed-point iteration ($p=1$), the scheme resulting from step \eqref{anderson_correction} and step \eqref{relaxation_step} is called \textit{Anderson-Richardson} (AR) \cite{Anderson, Potra, Walker}.
Studies about AR have highlighted similarities between this method and GMRES \cite{Fang, fang_two_2009, Potra, walker_anderson_2011}. 
When $p>1$, we refer to this approach as \textit{Alternating Anderson-Richardson} (AAR for short). Recent studies have shown also that AAR is more robust than AR against stagnations \cite{lupo_pasini_convergence_2019}. 
It should be noted that the least-squares problems solved to compute the Anderson mixing are not structured. This is a disadvantage in comparison with other standard linear solvers like GMRES \cite{FullGMRES}, where the least-squares problems are less expensive to solve because the matrices are structured as upper Hessenberg matrices. 

\subsubsection{Computational cost of AAR}
\label{cost_aar}
Denoting the number of non-zero entries per row of a sparse $n\times n$ matrix as $a$ and assuming for simplicity that it is constant across the rows, the computational complexity of a sparse matrix-vector multiplication is $\mathcal{O}(an)$. 
The least-squares problem in Equation \eqref{least_squares_matrix_form} used to perform the Anderson acceleration is computed using the QR factorization \cite{golub13} of the $n\times m$ matrix $R_k$. 
The main components used in Algorithm \ref{AAR_algorithm} to perform AA with their computational costs expressed with the big-$\mathcal{O}$ are:
\begin{itemize}
    \item QR factorization to perform AA: $\mathcal{O}(nm^2)$
    \item matrix-vector multiplication: $\mathcal{O}(an)$
\end{itemize}
AAR mitigates the computational cost to perform an AA step by interleaving successive AA step with $p$ relaxation sweeps. 
Relaxing the frequency at which the least-squares for Anderson is performed reduces the computational cost per iteration, but it is recommend not to excessively relax the frequency to avoid severely deteriorating the convergence rate of the overall numerical scheme. 
The computational complexity of AAR is
\begin{equation}
    \mathcal{O}\bigg( an + \frac{1}{p} nm^2\bigg).
\end{equation}
We can see the $\frac{1}{p}$ factor in front of the computational cost of the QR factorization allows to mitigate the impact that the least-squares solve has over the total computational cost of the numerical scheme.

In large-scale distributed computational environments, solving the global least-squares at periodic intervals may still not be sufficient to avoid severe bottlenecks in the computation. 
Reducing the computational cost to solve the least-squares is possible by replacing expensive and accurate calculations with inexpensive approximate ones while ensuring that convergence is maintained. 
To this end, in the following section we conduct an error analysis that provides theoretical bounds to estimate the error affordable at each iteration of AAR, which will allow us to develop numerical strategies to reduce the cost of the least-squares solves. 



\section{Backward stability analysis of approximate least-squares solves for AAR}
\label{section:error_bounds}
The solution obtained by solving approximately the least-squares problem \eqref{least_squares_matrix_form} can be treated as the exact solution of a different least-squares problem, which can be looked at as a yet unknown perturbation of the original problem \eqref{least_squares_matrix_form}. 
The object of this section is to estimate the difference between the original least-squares problem \eqref{least_squares_matrix_form} and its perturbation. This is important in order to ensure backward stability \cite{saunders, doi:10.1137/17M1157106}, which is a property of the numerical scheme to solve a problem as opposed to the forward stability, which is a property of the problem itself. Approximate calculations may arise from approximate evaluations of the fixed-point operator or from approximate calculations to solve the least-squares to perform an AA step. 

Analogous to the approach adopted in \cite{lupo_pasini_convergence_2019}, we provide the backward stability analysis under the assumption that the AAR is not truncated, i.e., $m=k$, which from herein we refer to as ``Full AAR". 
We focus on the backward stability analysis of Full AAR in the case that the mixing coefficients $\mathbf{g}^k$ solves an inexact or perturbed version of \eqref{least_squares_matrix_form}.
Specifically, we consider the following two cases:
(i) when the backward error affects only the matrix $R_k$ on the left-hand side (LHS) of the least squares system, which covers the scenario when the numerical methods performs an approximate factorization of $R_k$ (e.g., approximate QR or low-rank singular value decomposition), and
(ii) when the backward error affects both the LHS matrix $R_k$ and the residual vector $\mathbf{r}^k$ on the right-hand side (RHS), when occurs when the least-squares \eqref{least_squares_matrix_form} is approximated by projection onto a subspace. 
The backward stability analysis for cases (i) and (ii) are given in Sections~\ref{inaccuracy_lhs} and \ref{inaccuracy_both}, respectively.
In the following analysis, we assume that, in Full AAR, there is no stagnation that exceeds $p$ consecutive iterations, which is the number of Richardson sweeps between two AA steps.%
\footnote{
From the convergence results in \cite{lupo_pasini_convergence_2019}, we know that Full AAR converges in exact arithmetic only if the iterative scheme has not stagnated for more than $p$ consecutive steps.
If the iterative scheme has stagnated for more than $p$ consecutive steps, theoretical results in \cite{lupo_pasini_convergence_2019} indicate that Full AAR does not converge. 
In scientific applications where $m$ is fixed at $m<k$ and the history window of updates is cyclically erased, the theoretical results provided in this work can still be applied by restarting the iteration count from scratch every time the memory is erased.     
}

Let us denote with $\mathcal{E}_k \in \mathbb{R}^{n\times k}$ any perturbation to the matrix $R_k$ and $\delta \mathbf{r}_k \in\mathbb{R}^n$ any perturbation to the residual $\mathbf{r}_k$.
For convenience, we set $\omega=1$ in the analysis, but the obtained theoretical results are applicable to any $\omega\in[0,1]$.
The $i$th column of $\mathcal{E}_k$ can be modeled as an error matrix $E_i\in\mathbb{R}^{n\times n}$ applied to the $i$th columns of $R_k$, i.e.,
$    \mathcal{E}_k = [E_0{ (\mathbf{r}^1-\mathbf{r}^0)}, \ldots, E_k{  (\mathbf{r}^{k}-\mathbf{r}^{k-1})}]\:.
$
We then denote the perturbed LHS as
\begin{equation}
    \hat{R}_k = R_k + \mathcal{E}_k = R_k + [E_0{ (\mathbf{r}^1-\mathbf{r}^0)}, \ldots, E_k{  (\mathbf{r}^{k}-\mathbf{r}^{k-1})}].
\end{equation}
Using $\mathbf{r}^i - \mathbf{r}^{i-1} = -A[\mathbf{x}^i - \mathbf{x}^{i-1}]$ for $i=1,\ldots,k$, 
we obtain
\begin{equation}
    \hat{R}_k = R_k {  -} [E_0{ A}(\mathbf{x}^1-\mathbf{x}^0), \ldots, E_k{  A }(\mathbf{x}^{k}-\mathbf{x}^{k-1})].
    \label{eq:perturbed_LHS}
\end{equation}
As for the perturbed RHS, we denote it as 
\begin{equation}
    \hat{\mathbf{r}}_k = \mathbf{r}_k + \delta \mathbf{r}_k\:.
    \label{eq:perturbed_RHS}
\end{equation}

\subsection{Backward stability analysis of AAR with perturbed LHS}
\label{inaccuracy_lhs}
In this subsection, we consider the case that the mixing coefficients $\mathbf{g}^k$ in \eqref{anderson_correction} is given by 
\begin{equation}
  \mathbf{g}^k =\underset{\mathbf{g}\in \mathbb{R}^{m}}{\operatorname{argmin}}\lVert  \hat{R}_k  \mathbf{g} - \mathbf{r}_k  \rVert_2^2,
  \label{least_squares_perturbed_matrix}
\end{equation}
where the perturbed LHS matrix $\hat{R}_k$ is defined in \eqref{eq:perturbed_LHS}.
The optimal backward error of the least-squares in \eqref{eq:perturbed_LHS} \cite{saunders} is defined as 
\begin{equation}
    \min_{\mathcal{E}_k}\{ \lVert \mathcal{E}_k \rVert_F \; \text{:} \; \lVert  (R_k + \mathcal{E}_k)   \mathbf{g}^k - \mathbf{r}_k  \rVert_2^2 = \min\},
\end{equation}
where $\lVert \cdot \rVert_F$ denotes the Frobenius norm. 
When the Anderson mixing is performed at iteration $k=\ell p$, 
we define the perturbation from the backward error as 
\begin{equation}
    \delta_k = \| \mathcal{E}_k \mathbf{g}^k \|_2 = 
    \lVert [E_0{ A}(\mathbf{x}^1-\mathbf{x}^0), \ldots, E_k{ A}(\mathbf{x}^{k}-\mathbf{x}^{k-1})]{ \mathbf{g}^k} \rVert_2.
    \label{distance}
\end{equation}
Assume that $k=\ell p$ iterations of the AAR have been carried out with approximate residual evaluations and the iteration has not stagnated for more than $p$ consecutive steps. 
The perturbation $\delta_k$ is bounded by 
\begin{equation}
    \delta_k 
    \le \sum_{i=1}^k \lvert g_i^{(k)}\rvert \lVert E_i \rVert_2 { \lVert A \rVert_2} \lVert \mathbf{x}^i-\mathbf{x}^{i-1}\rVert_2,
\end{equation}
where the $\lVert \cdot \rVert_2$ of a matrix is defined as the matrix norm induced by the $\ell^2$-norm of a vector. 
Following the backward stability analysis, the final attainable residual norm should be monitored in terms of the relative quantity $$\displaystyle \frac{\lVert [E_0 { A} (\mathbf{x}^1-\mathbf{x}^0), \ldots, E_k{ A}(\mathbf{x}^{k}-\mathbf{x}^{k-1})]\mathbf{g}^k \rVert_2}{\lVert \mathbf{g}^k \rVert_2}.$$ It is thus natural to impose relative conditions on the norm of the error matrices. Given $\epsilon \in (0,1)$, we can consider the perturbation matrices $E_1, \ldots, E_k$, satisfying
\begin{equation}
\lVert [E_0{ A}(\mathbf{x}^1-\mathbf{x}^0), \ldots, E_k{ A}(\mathbf{x}^{k}-\mathbf{x}^{k-1})] \rVert_2 \le \epsilon { \lVert A \rVert_2} \lVert X^k \rVert_2 \le \sum_{i=1}^k\epsilon { \lVert A \rVert_2} \lVert \mathbf{x}^i-\mathbf{x}^{i-1} \rVert_2.    
\end{equation}
{  To simplify the presentation, we shall incorporate the norm of $A$ in the relaxation parameter, thus defining $\varepsilon = \epsilon \lVert A \rVert_2$, $\varepsilon \in (0,\lVert A \rVert_2)$}. 
The following theoretical results provide error bounds on the affordable perturbation introduced into the least-squares problem by approximate calculations to ensure that the final numerical solution to the linear system \eqref{linsys} drops the final residual below a user-defined threshold $\epsilon$.
\begin{lemma}
\label{lemma1}
Assume that $k=p\ell$ iterations of Full AAR have been carried out and the iteration has not stagnated for more than $p$ consecutive steps. Let ${ \mathbf{g}}^{k}$ be the Anderson mixing computed by solving the perturbed least-squares problem defined in Equation \eqref{least_squares_perturbed_matrix}. Then, the following inequality holds
\begin{equation}
    \lvert g_i^{(k)} \rvert \le  \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} \lVert \mathbf{r}^k \rVert_2 \quad \text{for any}\quad i= 1,\ldots,k. 
    \label{estimate_lemma1}
\end{equation}
\begin{proof}
We denote the QR factorization \cite{Golub} of the matrix $\hat{R}_k$ computed with a backward stable numerical method such as modified Gram-Schmidt or Householder transformations with
\begin{equation}
    \hat{R}_k = Q_k \hat{T}_k,
\end{equation}
where $Q_k \in \mathbb{R}^{n\times k}$ is a rectangular matrix with orthogonal columns and $\hat{T}_k \in \mathbb{R}^{k\times k}$ is an upper triangular matrix. 
By applying $(\hat{T}_k)^{-1} Q_k^T$ on both sides of the least-squares equation
$    \hat{R}_k { \mathbf{g}^k} = \mathbf{r}_k$,
we obtain 
\begin{equation}
    { \mathbf{g}^k} = (\hat{T}_k)^{-1} Q_k^T \mathbf{r}^k.
\end{equation}
It then follows that
\begin{equation}
\begin{aligned}
\lVert { \mathbf{g}^k}\rVert_2 & = \lVert (\hat{T}_k)^{-1} Q_k^T \hat{\mathbf{r}}^k \rVert_2 \le  
\lVert (\hat{T}_k)^{-1} \rVert_2  \lVert Q_k^T \hat{\mathbf{r}}^k \rVert_2 = \lVert (\hat{T}_k)^{-1} \rVert_2  \lVert \mathbf{r}^k \rVert_2
\end{aligned}
\end{equation}
where the last equality is obtained using the property of orthogonal matrices acting as isometric transformations. 
Since $\lvert { g_i^{(k)}} \rvert \le \lVert { \mathbf{g}^k} \rVert_2$ for $i=1,\ldots,k,$ and
\begin{equation}
    \lVert (\hat{T}_k)^{-1}  \rVert_2 = \left({\sigma_{\text{min}}(\hat{T}_k)}\right)^{-1}
\end{equation}
with $\sigma_{\text{min}}(\hat{T}_k)$ the smallest non-zero singular value of $\hat{T}_k$, 
the lemma follows.
\end{proof}
\end{lemma}

Lemma \ref{lemma1} supports the demonstration of the following theorem, which provides a dynamic criterion for adjusting the accuracy of linear algebra operations that involve the matrix $\hat{R}_k$ in solving the least-squares, thus allowing for more computationally inexpensive computations to solve \eqref{least_squares_perturbed_matrix} while maintaining accuracy. 

\begin{theorem}
\label{thm1}
Let $\epsilon>0$ and $\varepsilon = \epsilon \lVert A \rVert_2$. Let $\mathbf{r}^{k}$ be the residual of AAR after $k$ iterations. Assume that $k=p\ell$ iterations of Full AAR have been carried out and the iteration has not stagnated for more than $p$ consecutive steps. If for every $i\le k$,
\begin{equation}
    \lVert E_i \rVert_2 \le \frac{\sigma_{\text{min}}(\hat{T}_k)}{k} 
    \frac{1}{\lVert \mathbf{r}^{k}\rVert_2}
    \frac{1}{\lVert \mathbf{x}^{i} - \mathbf{x}^{i-1} \rVert_2}\epsilon,
    \label{error_bound1}
\end{equation}
then $\delta_k < \varepsilon$. 
\end{theorem}
\begin{proof}
The definition of $\delta_k$ in Equation \eqref{distance} and the upper bound for $\lvert g_i^{(k)}\rvert$ from Lemma \ref{lemma1} lead to the following chain of inequalities
\begin{equation}
     \delta_k  \le \sum_{i=1}^k \lvert { g_i^{(k)}}\rvert \lVert E_i \rVert_2  \lVert \mathbf{x}^i-\mathbf{x}^{i-1}\rVert_2 \le 
    \sum_{i=1}^k \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} \lVert \mathbf{r}^k \rVert_2 { \lVert A \rVert_2} \lVert E_i \rVert_2  \lVert \mathbf{x}^i-\mathbf{x}^{i-1}\rVert_2.
\end{equation}
Therefore, using the assumption \eqref{error_bound1},
we obtain 
\begin{equation}
     \delta_k  \le 
    \sum_{i=1}^k \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} \lVert \mathbf{r}^k \rVert_2 { \lVert A \rVert_2} \lVert E_i \rVert_2  \lVert \mathbf{x}^i-\mathbf{x}^{i-1}\rVert_2  \le  \sum_{i=1}^k \frac{1}{k} \epsilon \lVert A \rVert_2 = \varepsilon.
\end{equation}
\end{proof}
\begin{remark}
Lemma \ref{lemma1} and Theorem \ref{thm1} reproduce for Full AAR similar theoretical results presented in \cite{simoncini} for GMRES and full orthogonalization method (FOM). 
Theorem \ref{thm1} states that one can afford performing more and more inaccurate evaluations of the residual throughout successive iterations and still ensure that Full AAR converges with the final residual below a prescribed threshold $\varepsilon$.
In particular, when the calculation of the residual is approximate due to the restriction of the residual onto a subspace, the theorem states that the restriction can be more and more aggressive as convergence is reached. 
\end{remark}

\begin{remark}
Theorem \ref{thm1} is of important theoretical value, but the estimate is inconvenient for the following reasons: 
\begin{itemize}
\item the theorem does not provide an easily computable criterion to bound $\lVert E_i\rVert_2$ because the matrix $\hat{T}_k$ and the computed residual $\hat{\mathbf{r}}^k$ are not available at $i<k$ iterations.
\item $\hat{T}_k$ differs depending on the magnitude of the perturbations occurring in the approximate evaluations of the residual during $k$ steps.  
\item computing the minimum singular value $\sigma_{min}(\hat{T}_k)$ is as expensive as the computation we want to avoid and therefore defeats the purpose.
\end{itemize}
\end{remark}
To address the first impracticality of the error bound in Equation \eqref{error_bound1}, we provide the following corollary that replaces $\mathbf{r}^k$ with $\mathbf{r}^i$, which is computable at iteration $i<k$. 
\begin{corollary}
\label{corollary1}
Let $\epsilon>0$ and $\varepsilon = \epsilon \lVert A \rVert$. Assume that $k=p\ell$ iterations of Full AAR have been carried out and the iteration has not stagnated for more than $p$ consecutive steps. If for every $i\le k$,
\begin{equation}
    \lVert E_i \rVert_2 \le \frac{\sigma_{\text{min}}(\hat{T}_k)}{k} 
    \frac{1}{\lVert \mathbf{r}^{i}\rVert_2}
    \frac{1}{\lVert \mathbf{x}^{i} - \mathbf{x}^{i-1} \rVert_2}\epsilon,
        \label{error_bound3}
\end{equation}
then $\delta_k < \varepsilon$. 
\end{corollary}
\begin{remark}
We point out that the error bound \eqref{error_bound3} for any iteration $i$ relies on the iteration index $k$. Convenient values of $k=k^*$ that can be chosen are the maximum number of iterations allowed or the size of the matrix. 
\end{remark}
\begin{remark}
The singular value $\sigma_{min}(\hat{T}_k)$ could go to zero when converging since the columns may become linearly dependent. Projecting the residual matrix onto a subspace, as proposed later in Section~\ref{section:reduced_aar}, may help recover the linear independence of the columns and thus increase the minimum singular value of $\hat{T}_k$.
\end{remark}
The error bound in Equation \eqref{error_bound3} is still impractical because $\sigma_{\text{min}}(\hat{T}_k)$ is not available at iteration $i$, for $i<k$.  
We propose to replace $\sigma_{\text{min}}(\hat{T}_k)$ in \eqref{error_bound3} with a practical heuristic bound. Suppose we choose 1 as the heuristic bound.
This choice is analogous to what the authors propose in \cite{simoncini}, where it is mentioned that this heuristic bound may be occasionally too far from the theoretical bound, causing the true residual to reach its final attainable accuracy above the requested tolerance $\varepsilon$. To address this limitation, the authors in \cite{simoncini} propose replacing 1 with an inexpensive estimate that resorts to an approximation of the minimum non-zero singular value of the matrix $A$. Instead, we propose to combine the original heuristic bound that uses 1 at the numerator with a check on the monotonicity of the residual norm across consecutive AA steps. In fact, theoretical results in \cite{lupo_pasini_convergence_2019} ensure that Full AAR produces a sequence of AA updates that monotonically decreases the residual if the iteration has not stagnated for more than $p$ consecutive steps. If the residual norm stagnates between two consecutive AA updates, we use this as an indication that the last approximate AA step may have injected too large errors in the calculations. To attempt recovering the monotonic decrease, we thus propose to discard the last AA step, revert to the previous AA update, and increase the accuracy of the approximate calculations performed to solve the least-squares and obtain the new AA update. 
This consideration allows us to state the following (both rigorous and practical) theoretical result.
\begin{corollary}
\label{heuristic}
Assume that $k=p\ell$ iterations of Full AAR have been carried out and the iteration has not stagnated for more than $p$ consecutive steps. 
Let us denote the appropriately adjusted heuristic error estimate by 
\begin{equation}
B_i = \displaystyle \frac{\gamma_i}{k} 
    \frac{1}{\lVert \mathbf{r}^{i}\rVert_2}
    \frac{1}{\lVert \mathbf{x}^{i} - \mathbf{x}^{i-1}\rVert_2}, 
\end{equation}
where $0<\gamma_i<1$ for any iteration $i$ with $i\le k$.
If the following conditions hold: 
\begin{itemize}
    \item $\lVert E_i \rVert_2 \le B_i \epsilon$  for any $i$ with $i\le k$
    \item $\lVert \mathbf{r}^{k}\rVert_2 < \lVert \mathbf{r}^{\ell(p-1)}\rVert_2$ 
\end{itemize}
then $\delta_k<\varepsilon$.
\end{corollary}

\subsection{Backward stability analysis of AAR with perturbed LHS and RHS}
\label{inaccuracy_both}
In this subsection, we consider the case that the mixing coefficients $\mathbf{g}^k$ in \eqref{anderson_correction} is given by 
\begin{equation}
  \mathbf{g}^k =\underset{\mathbf{g}\in \mathbb{R}^{m}}{\operatorname{argmin}}\lVert  \hat{R}_k  \mathbf{g} - \hat{\mathbf{r}}_k  \rVert_2^2,
  \label{least_squares_perturbed_matrix_and_rhs}
\end{equation}
where $\hat{R}_k$ and $\hat{\mathbf{r}}_k$ are defined in \eqref{eq:perturbed_LHS} and \eqref{eq:perturbed_RHS}, respectively.
When both LHS and RHS of the least-squares are perturbed \cite{saunders}, the optimal backward error to solve the least-squares in Equation \eqref{least_squares_perturbed_matrix_and_rhs} at a generic iteration $k$ is
\begin{equation}
    \min_{\mathcal{E}_k, \delta\mathbf{r}^k} \{ \lVert \mathcal{E}_k, \delta\mathbf{r}^k \rVert_F \; \text{:} \; \lVert  (R_k + \mathcal{E}_k)   \mathbf{g}^k - (\mathbf{r}_k+\delta\mathbf{r}^k)  \rVert_2^2 = \min \}.
\end{equation}
The analysis conducted in this section assumes that the error committed on the RHS of the least-squares can be tuned in a practical way. For instance, when the least-squares is projected onto a subspace, the distance between the original residual vector and the projected residual vector, namely $\lVert \delta \mathbf{r}^k \rVert_2$, can be measured in a straightforward way and the projection techniques can be adjusted accordingly. As for the error committed on the LHS, namely $\delta_k$, this is less straightforward to control as it generally results from a sequence of several manipulations performed on the matrix $R_k$. The backward error analysis performed in this section gives the requirements on the inaccuracy of the operations performed on the matrix $R_k$, given a prescribed inaccuracy in the evaluation of the RHS. 
In this context, we define for convenience
\begin{equation}
    \tilde{\varepsilon} = \epsilon \max( \lVert A \rVert_2, \lVert \mathbf{r}_k \rVert_2)
\end{equation}
and we study the conditions required to ensure that $\delta_k<\tilde{\varepsilon}$. 

To obtain an upper bound on $|g_i^{(k)}|$, we apply the analysis in Lemma~\ref{lemma1} to \eqref{least_squares_perturbed_matrix_and_rhs} and obtain, for $i= 1,\ldots,k$,
\begin{equation}
    \lvert g_i^{(k)} \rvert \le  \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} \lVert \hat{\mathbf{r}}^k \rVert_2 \leq \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} \bigg ( \lVert \mathbf{r}^k \rVert_2 + \lVert \delta \mathbf{r}^k \rVert_2 \bigg),
\end{equation}
where the second inequality follows from \eqref{eq:perturbed_RHS} and the triangle inequality.
Assuming that $\|\delta \mathbf{r}^k\|_2\leq \epsilon\|\mathbf{r}_k \|_2$, we obtain 
\begin{equation}
    \lvert g_i^{(k)} \rvert 
    \le \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} (1+\epsilon) \lVert \mathbf{r}_k \rVert_2,  \quad \text{for any} \quad i= 1,\ldots,k.
    \label{eq:g_bound}
\end{equation}
By following the same linear algebra steps applied in Theorem \ref{thm1} where we replace $\lVert \mathbf{r}_k \rVert_2$ with $\lVert \mathbf{r}_k \rVert_2 + \lVert \delta \mathbf{r}_k \rVert_2$, we obtain the following theoretical result.


\begin{theorem}
\label{thm2}
Let $\epsilon>0$ and let 
$\varepsilon = \epsilon \max (\lVert \mathbf{r}_k \rVert_2 , \lVert A \rVert_2)$. 
Let $\hat{\mathbf{r}}^{k}$ be the residual of AAR after $k$ iterations. Under the same hypotheses and notation of Lemma \ref{lemma1} and assuming that $\|\delta \mathbf{r}^k\|_2\leq \epsilon\|\mathbf{r}_k \|_2$, if for every $i\le k$,
\begin{equation}
    \lVert E_i \rVert_2 \le \frac{\sigma_{\text{min}}(\hat{T}_k)}{k} 
    \frac{1}{\lVert \mathbf{r}^{k} \rVert_2}
    \frac{1}{\lVert \mathbf{x}^{i} - \mathbf{x}^{i-1} \rVert_2}\frac{\epsilon}{1+\epsilon},
    \label{error_bound2}
\end{equation}
then $\delta_k < \tilde\varepsilon$. 
\end{theorem}
\begin{proof}
The definition of $\delta_k$ in Equation \eqref{distance} and the upper bound for $\lvert g_i^{(k)}\rVert_2$ in \eqref{eq:g_bound} lead to the following inequality 
\begin{equation}
\begin{aligned}
     \delta_k  & 
         \le  \sum_{i=1}^k \frac{1}{\sigma_{\text{min}}(\hat{T}_k)} (1+\epsilon) \lVert \mathbf{r}^k \rVert_2 { \lVert A \rVert_2} \lVert E_i \rVert_2  \lVert \mathbf{x}^i-\mathbf{x}^{i-1}\rVert_2\:.
\end{aligned}
\end{equation}
Therefore, applying the assumption \eqref{error_bound2} to the right-hand size of above then gives
\begin{equation}
 \delta_k 
 \le \sum_{i=1}^k \frac{1}{k} (1+\epsilon) \frac{\epsilon}{(1+\epsilon)} \lVert A \rVert_2 
 \le \tilde\varepsilon.
\end{equation}
\end{proof}


Theorem \ref{thm2} states that one can afford performing more and more inaccurate evaluations of the residual throughout successive iterations and still ensure that Full AAR drops the residual norm below a prescribed threshold $\tilde{\varepsilon}$. 
In particular, when the calculation of the residual is approximate due to the restriction of the residual onto the subspace, the theorem states that the restriction can be more and more aggressive as convergence is reached. 

\begin{remark}
Theorem \ref{thm2} is valid under a slightly more stringent requirement with respect to the theoretical bound provided in Theorem \ref{thm1}. This is reasonably expected because the hypotheses of Theorem \ref{thm2} assume that the perturbation both affects LHS and RHS, differently from the hypotheses of Theorem \ref{thm1} where only the LHS was assumed to be perturbed. Since the hypothesis of Theorem \ref{thm2} assumes that more sources of error are injected in the original definition of the least-squares, a slightly stronger requirement needs to be imposed on the backward error to preserve backward stability in performing the Anderson acceleration.
\end{remark}

The considerations provided in section~\ref{inaccuracy_lhs} that motivate the need for a heuristic quantity to dynamically adjust the inaccuracy of the least-squares solve can be extended to this section as well. 

\section{Reduced Alternating AA}
\label{section:reduced_aar}
Now we are ready to use the guidelines of the theoretical results to propose a new variant of AA which reduces the computational effort to compute the Anderson mixing by judiciously projecting the least-squares onto a subspace. We call this method \textit{Reduced Alternating AA}.
At iteration $k=\ell p$, Reduced Alternating AA computes approximate mixing coefficients $\hat{\mathbf{g}}^k$ by solving
\begin{equation}
\hat{\mathbf{g}}^k = \underset{\mathbf{g}\in \mathbb{R}^{m}}{\operatorname{argmin}}\lVert 
P_k (R_k\mathbf{g} - \mathbf{r}^k)\rVert^2_{2} \:,   
\label{eq:projected_LS}
\end{equation}
where $P_k\in\mathbb{R}^{n\times n}$ of rank $s<n$ is a projection operator onto an $s$-dimensional Euclidean subspace $\mathcal{S}_k\subset \mathbb{R}^n$.
To keep the computation cost of the projection low, we focus on the simplified case in which the projection is restricted to a row selection procedure, i.e., 
\begin{equation}
    P_k = S_k S_k^T\quad\text{with}\quad S_k =[\mathbf{e}_{j_1}, \ldots, \mathbf{e}_{j_s}]\in \mathbb{R}^{n\times s}.
\end{equation}
Here $\mathbf{e}_{j}$ denotes the $j$-th canonical basis, and $S_k$ is the restriction operator that selects only the $j_i$-th rows, $i=1,\dots,s$.
In this paper, we propose the following two strategies to select the index set $\{j_i\}_{i=1}^s$.
\begin{itemize}[leftmargin=20pt]
    \item {\it Subselected Alternating AA} -- the index set $\{j_i\}_{i=1}^s$ is selected to be the $s$ indices corresponding to the $s$ entries in the residual vector $\mathbf{r}^k$ with largest magnitudes.
    \item {\it Randomized Alternating AA} -- the index set $\{j_i\}_{i=1}^s$ is drawn from the full index set $\{1,\dots,n\}$ under a uniform distribution without replacement. 
\end{itemize}
With the same size $s$, the two strategies lead to similar reduction in terms of the computation time for solving \eqref{eq:projected_LS}, with minor difference in the computation overhead in the index selection process. 
The dimensionality $s$ of the projection subspace $\mathcal{S}_k$ is adaptively tuned at each according to values of the heuritstic that estimates the error bound in Equation \eqref{error_bound2}.

\subsection{Computational cost of Reduced Alternating AA for linear fixed-point iterations}
Calling $s$ the total number of entries of the residual vector retained for the projection of the least-squares onto a subspace, since $s\le n$, the computational cost to solve the least-squares reduces from $\mathcal{O}(nm^2)$ to $\mathcal{O}( s m^2)$.
Therefore, the computational complexity of Reduced AAR is less than the computational cost of AAR
\begin{equation}
    \mathcal{O}\bigg( an + \frac{1}{p} s m^2\bigg) \le \mathcal{O}\bigg( an + \frac{1}{p} n m^2\bigg).
\end{equation}
While the dimension $s$ of the subspace $\mathcal{S}_k$ is tuned according to the heuristic proposed in Equation \eqref{error_bound2}, the value of $p$ still needs to be empirically and judiciously tuned. While increasing $p$ reduces the total computational cost of the scheme, it may also expose the scheme to a higher risk of numerical instabilities and deteriorated convergence rate. 

\section{Numerical results}
\label{section:numerical_results}

In this section we present numerical results where fixed-point iterations are performed inaccurately, but the Anderson acceleration still succeeds in converging by bringing the residual of the fixed-point iteration below a prescribed tolerance. We focus on two types of fixed-point iterations: linear deterministic and non-linear deterministic. 
The numerical results for the linear deterministic fixed-point iteration describe the use of AA to accelerate the Richardson scheme to solve linear systems. The numerical results for the non-linear deterministic fixed-point iteration describe the use of AA to accelerate the solution of an implicit system from a simplified Boltzmann equation at a given time-step.


\subsection{Linear deterministic fixed-point iteration: iterative linear solver}
The numerical results used to illustrate AA on linear deterministic fixed-point iterations to solve linear systems are split in two parts. First, we study how the convergence of AR is affected by random perturbation with tunable intensity applied to the LHS of the least-squares solved to compute the Anderson mixing. Then, we show how restricting the least-squares onto a projection subspace reduces the computational time to solve a linear system while still converging.    

\subsubsection{Injection of random noise on LHS of least-squares for each AA step}
We consider the diagonal matrix $A = \text{diag}(10^{-4}, 2, 3, \ldots, 100)$. 
We add random noise the the LHS of the least squares to compute the Anderson mixing:
\begin{equation}
    \bigg (R_k + \epsilon_k \lVert R_k\rVert_2 \lVert E_k \rVert_2  \bigg)\mathbf{g} = \mathbf{r}^k
\end{equation}
The perturbation matrices $E_k$'s, with $\lVert E_k\rVert_2=1$ for $k>0$ are random $100 \times 100$ matrices generated with normally distributed values using the \texttt{Matlab} function \textbf{randn}. The quantity $\epsilon_k$ is defined as 
\begin{equation}
    \epsilon_k =  \frac{\epsilon}{k^*} \frac{\sigma_\text{min} (T_k)}{\lVert \mathbf{r}^k \rVert_2 \lVert \mathbf{x}^k - \mathbf{x}^{k-1} \rVert_2 },
    \label{test_error}
\end{equation}
where $T_k$ is the triangular factor of the QR factorization of $R_k$. 
The quantity $\epsilon_k$ is used to tune the magnitude of the noise injected in the LHS of Equation \eqref{test_error}, and we set $k^* = 100$ according to the size of the linear system we are solving. We use four different values of $\epsilon$ to tune the amount of noise injected in the LHS of the least-squares, $\epsilon = 1e-8$, $\epsilon = 1e-6$, $\epsilon = 1e-4$. $\epsilon = 1.0$. 
Numerical results in Figure \ref{random_noise} show the magnitude $\epsilon_k$ of the random noise injected in the LHS as a function of the iteration count. The residual norm history of AR for both the unperturbed least-squares and the perturbed least-squares is shown as well to illustrate how the random perturbation affects the convergence of AR for each value of $\epsilon$ tested. As expected, increasing values of $\epsilon$ lead to a progressive deterioration of the convergence rate, but final convergence within the prescribed tolerance of the residual is still preserved. 

\begin{figure}
  \centering
  \subfloat[]{\label{profile_1e8}\includegraphics[width=0.5\textwidth]{./figures/eps_1e8}}
  \subfloat[]{\label{profile_1e6}\includegraphics[width=0.5\textwidth]{./figures/eps_1e6}} 
  \newline
  \subfloat[]{\label{profile_1e4}\includegraphics[width=0.5\textwidth]{./figures/eps_1e4}}
  \subfloat[]{\label{profile_1e0}\includegraphics[width=0.5\textwidth]{./figures/eps_1e0}}  
  \caption{\label{random_noise}Approximate AA convergence history.}
\end{figure}

\subsubsection{Reduced Alternating AA to solve sparse linear systems}

We now use Reduced AAR to solve a set of linear systems where the matrices are open source and pulled from the SuiteSparse Matrix Collection (formerly known as the University of Florida Sparse Matrix Collection) \cite{UFL}, the Matrix Market Collection \cite{MM}. 
In Table \ref{matrices}, we report the matrices and their most significant properties. The sources used to retrieve the matrices are specified in Table \ref{matrices} as well. The notation MM is used to refer to the Matrix Market Collection, and SS for the SuiteSparse Matrix Collection. 

\begin{table}
\center
\begin{tabular}{|c|c|c|c|c|c|}
    \hline
        \textbf{Matrix} & \textbf{Type} & \textbf{Size} & 
         \textbf{Structure} & \textbf{Pos. def.} & \textbf{Source}\\
    \hline  
    fidap029 & real & 2,870 & nonsymmetric & yes & MM\\
    \hline
    raefsky5 & real & 6,316 & nonsymmetric & yes & SS\\
    \hline
    bcsstk29* & real & 13,992 & symmetric & no & SS\\    
    \hline
    sherman3* & real & 5,005 & nonsymmetric & no & MM \\
    \hline
    sherman5* & real & 3,312 & nonsymmetric & no & MM\\
    \hline
    fidap008* & real & 3,096 & nonsymmetric & no & MM \\
    \hline
    chipcool0 & real & 20,082 & nonsymmetric & no & SS\\
    \hline
    e20r0000* & real & 4,241 & nonsymmetric & no & MM \\
    \hline
    spmsrtls & real & 29,995 & nonsymmetric &no & SS \\
    \hline
    garon1* & real & 3,175 & nonsymmetric & no & SS \\
    \hline
    garon2* & real & 13,535& nonsymmetric & no & SS\\
    \hline
    memplus & real & 17,758 & nonsymmetric & no & SS\\
    \hline
    saylr4 & real & 3,564 & nonsymmetric & no & MM \\     
    \hline
    xenon1* & real & 48,600 & nonsymmetric & no & SS \\
    \hline
    xenon2* & real & 157,464 & nonsymmetric & no & SS\\
    \hline
    venkat01 & real & 62,424 & nonsymmetric & no & SS \\
    \hline
    QC2534 & complex & 2,534 & non-Hermitian & no & SS \\
    \hline
    mplate* & complex & 5,962 & non-Hermitian & no & SS \\
    \hline
    light\_in\_tissue & complex & 29,282 & non-Hermitian & no & SS \\
    \hline
    kim1 & complex & 38,415 & non-Hermitian & no & SS \\
    \hline
    chevron2 & complex & 90,249 & non-Hermitian & no & SS \\   
    \hline    
\end{tabular}
\caption{List of matrices used for numerical experiments performed in \sc{Matlab}.}
\label{matrices}
\end{table}

The experiments compare the performance of GMRES, Alternating AA, Subselected Alternating AA and Randomized Alternating AA. For each problem, two preconditioning techniques are used: Zero Fill-In Incomplete LU Factorization (ILU(0) for short) and Incomplete LU Factorization with Threshold (ILUT($\tau$) for short) \cite{Chan1997}.
When ILUT($\tau$) is used, the tolerance $\tau$ controls the level of sparsity in the incomplete LU factors. We set $\tau=10^{-4}$ and any zeros on the diagonal of the upper triangular factor are replaced by the local drop tolerance. The ILUT($\tau$) preconditioner is computed with a column pivoting and the pivot is chosen as the maximum magnitude entry in the column.
For more details about ILU(0) and ILUT($\tau$) we refer to \cite{Benzi}, \cite[pp.~287--307]{Saad}. 
A symmetric reverse Cuthill-McKee reordering \cite{rcm} has been applied to some matrices to allow a stable construction of the ILU factors.
The matrices that needed a reordering for the ILU preconditioner are marked with an asterisk close to their name in Table \ref{matrices}. 
Since the matrices \texttt{xenon1} and \texttt{xenon2} were poorly scaled, the ILU factors have been computed for these matrices
only after a diagonal scaling. The diagonal scaling has been applied via a diagonal preconditioner on these matrices. Those situations where the matrices had zeros on the main diagonal were treated by setting to one the associated entries in the diagonal preconditioner.
The Richardson relaxation parameter is set to $\omega = 0.2$.  The number of Richardson iterations between two Anderson updates for the AAR is $p=3$. 
The least-squares problems related to Anderson mixing steps
are solved via QR factorization with column pivoting of the rectangular matrix to the LHS. 
The RHS of the linear systems is obtained by multiplying the solution vector by the coefficient matrix $A$. The variation across single runs never exceeded $10\%$ in terms of computational time. 
A threshold of $10^{-8}$ is used as a stopping criterion on the relative residual $\ell^2$-norm.
The heuristic proposed in Corollary \ref{heuristic} is used to automatically tune the dimensionality of the projection subspace at each AA step. To avoid introducing excessive computational overhead in fine tuning the dimensionality of the subspace, we add rows in batches, with each batch being 10\% of the total number of rows in the original least-squares problem. 

The performance of the linear solvers is evaluated through \textit{performance profiles} \cite{Dolan}. Let us refer to $\mathcal{S}$ as the set of solvers and $\mathcal{P}$ as the test set. We assume that we have $n_s$ solvers and $n_p$ problems. In this paper, performance profiles are used to compare the computational times. To this end we introduce
\[
t_{p,s}=\text{computing time to solve problem}\; p \; \text{with solver}\; s.
\]
The comparison between the times taken by each solver is based on the performance ratio defined as
\[
r_{p,s} = \frac{t_{p,s}}{\min\{t_{p,s}:s\in\mathcal{S}\}}.
\]
The performance ratio allows one to compare the performance of solver $s$ on problem $p$ with the best performance by any solver to address the same problem $p$. In case a specific solver $s$ does not succeed in solving problem $p$, then a convention is adopted to set $r_{p,s}=r_{M}$ where $r_M$ is a maximal value. In our case we set $r_m=10,000$. 
The performance of one solver compared to the others' on the whole benchmark set is displayed by the cumulative distribution function $\rho_s(\tau)$ that is defined as follows:
\[
\rho_s(\tau)=\frac{1}{n_p}\text{size}\{p\in\mathcal{P}:r_{p,s}\le \tau \}.
\]
The value $\rho_s(\tau)$ represents the probability that solvers $s\in\mathcal{S}$ has a performance ratio $r_{p,s}$ less than or equal to the best possible ratio up to a scaling factor $\tau$. 
The convention adopted that prescribes $r_{p,s}=r_M$ if solver $s$ does not solve problem $p$ leads to the reasonable assumption that
\[
r_{p,s}\in[1,r_M].
\]
Therefore, $\rho_s(r_M)=1$ and 
\[
\rho_s^* = \lim_{\tau\rightarrow r_M^-}\rho_s(\tau)
\]
represents the probability that solver $s$ succeed in solving a generic problem from set $\mathcal{P}$. 
To improve the clarity of the graphs, we replace $\rho_s(\tau)$ with
\[
\tau \mapsto \frac{1}{n_p}\text{size}\{p\in\mathcal{P}:\log_2(r_{p,s})\le \tau\}
\] 
and we use this quantity for each performance profile displayed.\newline

The numerical results with ILU(0) and ILUT($\tau$) are shown in Figures \ref{profile_ilu0} and \ref{profile_ilut}, respectively. Randomized Alternating AA outperforms Subselected Alternating AA by converging on a larger set of problems, and outperforms standard AA and GMRES for time-to-solution. Therefore, an adaptive reduction of the dimensionality of the least-squares problem is shown to effectively reduce the total computational cost of AAR to iteratively solve a broad class of sparse linear systems while still maintaining convergence. 

\begin{figure}
  \centering
  \subfloat[]{\label{profile_ilu0}\includegraphics[width=0.5\textwidth]{./figures/ilu0_preconditioner.png}}
  \subfloat[]{\label{profile_ilut}\includegraphics[width=0.5\textwidth]{./figures/ilutp_preconditioner.png}}  
  \caption{\label{profiles}Performance profiles in a $\log_2$ scale.}
\end{figure}

\subsection{Non-linear deterministic fixed-point iteration: Picard iteration for non-linear time-dependent Boltzmann equation}
\newcommand{\vect}[1]{\boldsymbol{#1}}
\def\mathbb{R}{\mathbb{R}}
\def\mathbb{S}{\mathbb{S}}
\def\textsf{c}{\textsf{c}}
\newcommand{\mbox{\tiny EA}}{\mbox{\tiny EA}}
\newcommand{\mbox{\tiny ES}}{\mbox{\tiny ES}}
\newcommand{\mbox{\tiny IS}}{\mbox{\tiny IS}}
\newcommand{\mbox{\tiny IN}}{\mbox{\tiny IN}}
\newcommand{\mbox{\tiny OUT}}{\mbox{\tiny OUT}}
\newcommand{\TOTAL}{\mbox{\tiny{Tot}}}       

To test Alternating AA with projected least-squares on deterministic non-linear problems, we consider the non-linear systems that arise from implicit time discretization of the non-relativistic Boltzmann equation, which is a widely used model for neutrino transport in nuclear astrophysics applications \cite{MEZZACAPPA1999281}. 
The non-relativistic Boltzmann equation is given by
\begin{equation}\label{eq:Boltzmann}
  \partial_t{f} + \mathcal{T}(f) = \mathcal{C}(f)\:,
\end{equation}
where $f(\vect{x},\omega,\epsilon,t)$ denotes the neutrino distribution function that describes the density of neutrino particles at position $\vect{x}\in\mathbb{R}^3$
traveling along direction $\omega\in\mathbb{S}^2$ with energy $\epsilon\in\mathbb{R}^{+}$
at time $t\in\mathbb{R}^+$.
Here the advection and collision operators are denoted by $\mathcal{T}$  and $\mathcal{C}$, respectively. 
Implicit-explicit (IMEX) time integration schemes \cite{pareschiRusso_2005} are often used to solve \eqref{eq:Boltzmann}, in which the collision term is treated implicitly to relax the excessive explicit time-step restriction, while an explicit advection term is used to avoid spatially coupled non-linear systems.
With a simple backward Euler method, the implicit stage of an IMEX scheme applied to \eqref{eq:Boltzmann} at time $t^n$ takes the form
\begin{equation}\label{eq:Implicit}
    f(\vect{x},\omega,\epsilon,t^{n+1}) = f(\vect{x},\omega,\epsilon,t^{n}) + \Delta t \,\mathcal{C}(f(\vect{x},\cdot,\cdot,t^{n+1}))\:, \quad\forall\vect{x}\in\mathbb{R}^3\:,
\end{equation}
where we follow the approach in \cite{Paul_Laiu_2021} and write the (space-time homogeneous) neutrino collision operator $\mathcal{C}$ as 
\begin{equation}
    \mathcal{C}(f) = {\eta_{\TOTAL}}(f) -  {\chi_{\TOTAL}}(f)\,f\:,
\end{equation}
where ${\eta_{\TOTAL}}$ denotes the total emissivity and ${\chi_{\TOTAL}}$ denotes the total opacity.
As in \cite{Paul_Laiu_2021}, we formulate \eqref{eq:Implicit} at each time step into a non-linear fixed-point problem
\begin{equation}\label{eq:fixed_point_Boltzmann}
    f^{n+1} = G(f^{n+1})\, \text{ with }\,
    G(f^{n+1}):=({f^{n}+\Delta t \eta_{\TOTAL}(f^{n+1})})/({1+\Delta t \chi_{\TOTAL}(f^{n+1})})\:.
\end{equation}
This formulation ensures that the fixed-point operator $G$ is a contraction and allows for the standard Picard iteration as well as variants of AA to solve \eqref{eq:Implicit}.

In the numerical tests, we solve the spatially-homogeneous implicit system \eqref{eq:fixed_point_Boltzmann} at one time step, in which we apply a discrete ordinate angular discretization with a 110-point Lebedev quadrature rule on $\mathbb{S}^2$ and 64 geometrically progressing energy nodes in $[0,300]$ (MeV). 
Here we compare the performance of five iterative solvers: Picard iteration, standard AA, Alternating AA, and two versions of Reduced Alternating AA -- Subselected Alternating AA and Randomized Alternating AA. 
The history length $m=3$ is used for all four AA variants, and the alternating ones perform Anderson mixing every $p=3$ iterations. The Subselected Alternating AA and Randomized Alternating AA solvers are implemented as presented in Section~\ref{section:reduced_aar} with the reduced dimension $s$ chosen to satisfy the bound \eqref{error_bound2}, where $\epsilon=10^{-8}$ is used.
Each solver is tested on problems at six different matter density values. In general, higher matter density corresponds to stronger collision effect, which makes the implicit system more stiff and require more iteration to converge. 
In the test, the total emissivity and opacity are computed using opacity kernels provided in
\cite{bruenn_1985}.

Figure~\ref{fig:Boltzmann} reports the iteration count and computation time required for each compared solver at various matter density. These reported results are averaged over 30 runs. Figure~\ref{fig:Boltzmann_Iter} confirms that AA and Alternating AA require significantly fewer iterations to reach the prescribed $10^{-10}$ relative tolerance than the Picard iteration, especially for problems at higher densities (more stiff). Randomized Alternating AA essentially resembles the iteration counts of Alternating AA, while the Subselected Alternating AA requires much higher iteration counts than other AA variants. 
Figure~\ref{fig:Boltzmann_Time} illustrates that AA actually takes more computation time than Picard iteration, which indicates that the additional cost of solving least squares problems outweighs the reduction of iteration counts. Alternating AA lowers the computation time of AA by reducing the frequency of least square solves, and Randomized Alternating AA further improves the computation time by approximating the least square solve with solution to a sparsely projected smaller problem.



\begin{figure}[h]
  \centering
  \subfloat[]{\label{fig:Boltzmann_Iter}\includegraphics[width=0.475\textwidth]{./figures/Boltzmann_Iter.pdf}}~~
  \subfloat[]{\label{fig:Boltzmann_Time}\includegraphics[width=0.5\textwidth]{./figures/Boltzmann_Time.pdf}}  
  \caption{\label{fig:Boltzmann}Iteration counts and computation time (in $\log$ scale) for various fixed-point solvers for solving the non-linear system \eqref{eq:fixed_point_Boltzmann} at different matter density, which corresponds to the stiffness of the system.}
\end{figure}


\begin{comment}
\subsection{Non-linear stochastic fixed-point iteration: data-batched stochastic optimization for DL training}

As an example of nonlinear stochastic fixed-point iteration, we train a graph convolutional neural network (GCNN) model with batched stochastic optimization to predict the total mixing enthalpy, atomic charge transfer, and atomic magnetic moment for the ferromagnetic solid solution alloy iron-platinum as a function of the atomic configuration.  
The iterative training of a DL model can be recast as a stochastic non-linear fixed-point scheme \cite{Nonlinear_fixed_point}, where one iteration consists in a forward-backward pass of a stochastic first-order optimizer to update the parameters of a DL model.

The dataset used is open-source and provides density functional theory (DFT) calculations of the target properties for 32,000 randomized configurations, each described by 32 atoms arranged in a body centered tetragonal (BCT) crystal structure with fixed volume that spans the entire compositional range of the binary alloy. The data pre-processing normalizes every target property to restrain the values within the range $[0,1]$.
The training set represents 70\% of the total dataset; the validation and test sets each represent half of the remaining data. The compositional stratified splitting ensures that all the compositions are equally represented across in the training, validation, and testing portions of the dataset. The training loss function is the mean squared error (MSE).

We use Adam method \cite{adam} as stochastic first-order optimizer to train the GCNN model. 
The first 30 epochs are used warm-up the training of the GCNN model using only the Nesterov optimizer to enable a quick preliminary exploration of the loss function landscape, and AA starts being used one the preliminary 30 epochs are completed. 
The frequency at which AA is performed is $p=5$, the history depth is $m=20$, and the relaxation parameter is $\omega=0.1$. 
The numerical simulations compare the unaccelerated Adam, Adam accelerated with standard AA, Subselected AA and Randomized AA where the size of the projection subspace is allowed to vary between 1\% and 20\% of the original size. 

The numerical results in Table \ref{GCNN_table} show that AA reduces the root mean-squared error (RMSE) of the trained HydraGNN model for total mixing enthalpy, atomic charge transfer and atomic magnetic moment. Moreover, the projection of the least-squares for each AA step onto a projection subspace allows to significantly reduce the training time without affecting the accuracy. 

\begin{table}
\begin{tabular}{ |p{4.5cm}||p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|  }
 \hline
 & \multicolumn{3}{|c|}{\textbf{RMSE}} & & \\
 \hline
 \textbf{Training scheme} & \textbf{H} & \textbf{C} & \textbf{M} & \textbf{T} (sec.)& \textbf{LS-T} (sec.)\\
 \hline
 Adam   &  0.0095  & 0.0085 &  0.0136 & 717 & - \\
 Adam + AA & 0.0082 & 0.0072 & 0.0102 & 2,272 & 0.263\\
 Adam + Subselected AA & 0.0085 & 0.0070 & 0.0102 & 1,629 & 0.028\\
 Adam + Randomized AA & 0.0087 & 0.0072 & 0.0102 & 1,623 & 0.026\\
 \hline
\end{tabular}
\caption{RMSE of the normalized total mixing enthalpy (H), atomic charge transfer (C), and atomic magnetic moment (M), total wall-clock time for training (T) in seconds, and the averages wall-clock time in seconds to solve one least-square (LS-T) to train a GCNN model on the FePt dataset with Adam, Adam combined with AA, and Adam combined with Reduced AA. }
 \label{GCNN_table}
\end{table}

\end{comment}

\section{Conclusion and future development}
\label{section:conclusions}
Although AA has been shown to significantly improve the convergence of fixed-point iterations in several scientific applications, effectively performing AA for large scale problems without introducing excessive computational burden still remains a challenge. One possibility to reduce the computational cost consists in projecting the least-squares to compute the Anderson mixing onto a projection subspace, but this may compromise the convergence of the scheme. 

We derived rigorous theoretical bounds for AA that allow for efficient approximate calculations of the residual, which reduce computational time and memory storage while maintaining convergence. The theoretical bound provide useful insight to judiciously inject inaccuracy in the calculations whenever this allows to alleviate the computational burden without compromising the final accuracy of the physics based solver. In particular, we consider situations where the inaccuracy arises from projecting the least-squares problem onto a subspace. 

Guided by the theoretical bounds, we constructed a heuristic that  dynamically adjusts the dimension of the projection subspace at each iteration. As the process converges (as it is guaranteed to in absence of stagnations), the backward error decreases, which allows for an increase of the inaccuracy of the least-squares calculations while still maintaining the residual norm within the prescribed bound. This reduces the need for very accurate calculations needed to converge. 
While previous theoretical results provide error estimates when the accuracy is fixed a priori throughout the fixed-point scheme, thus potentially compromising inevitably the final attainable accuracy, our study allows to adaptively adjust the accuracy to perform AA while still ensuring that the final residual norm drops below a prescribed threshold. 
The resulting numerical scheme, called Reduced Alternating AA, is clearly preferable over standard Alternating AA when the computational budget is limited.

Numerical results on linear systems and non-linear non-relativistic Boltzmann equation show that our heuristic can be used to effectively reduce the computational time of Reduced Alternating AA without affecting the final attainable accuracy compared to traditional methods to apply AA. 

Future work will be dedicated to apply Reduced Alternating AA to non-linear stochastic fixed point iterations arising from: (i) Monte Carlo methods for solving the Boltzmann equation and (ii) training of physics informed graph convolutional neural networks for the prediction of materials properties from atomic information. 

\section*{Acknowledgement}
Massimiliano Lupo Pasini thanks Dr. Vladimir Protopo-pescu for his valuable feedback in the preparation of this manuscript.
This work was supported in part by the Office of Science of the Department of Energy, by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration, and by the Artificial Intelligence Initiative as part of the Laboratory Directed Research and Development (LDRD) Program of Oak Ridge National Laboratory managed by UT-Battelle, LLC for the US Department of Energy under contract DE-AC05-00OR22725. 


\bibliographystyle{siamplain}
