Anderson acceleration

In mathematics, Anderson acceleration, also called Anderson mixing, is a method for the acceleration of the convergence rate of fixed-point iterations. Introduced by Donald G. Anderson,^[1] this technique can be used to find the solution to fixed point equations ${\displaystyle f(x)=x}$ often arising in the field of computational science.

Given a function ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} ^{n}}$, consider the problem of finding a fixed point of ${\displaystyle f}$, which is a solution to the equation ${\displaystyle f(x)=x}$. A classical approach to the problem is to employ a fixed-point iteration scheme;^[2] that is, given an initial guess ${\displaystyle x_{0}}$ for the solution, to compute the sequence ${\displaystyle x_{i+1}=f(x_{i})}$ until some convergence criterion is met. However, the convergence of such a scheme is not guaranteed in general; moreover, the rate of convergence is usually linear, which can become too slow if the evaluation of the function ${\displaystyle f}$ is computationally expensive.^[2] Anderson acceleration is a method to accelerate the convergence of the fixed-point sequence.^[2]

Define the residual ${\displaystyle g(x)=f(x)-x}$, and denote ${\displaystyle f_{k}=f(x_{k})}$ and ${\displaystyle g_{k}=g(x_{k})}$ (where ${\displaystyle x_{k}}$ corresponds to the sequence of iterates from the previous paragraph). Given an initial guess ${\displaystyle x_{0}}$ and an integer parameter ${\displaystyle m\geq 1}$, the method can be formulated as follows:^{[3][note 1]}

${\displaystyle x_{1}=f(x_{0})}$

${\displaystyle \forall k=1,2,\dots }$

${\displaystyle m_{k}=\min\{m,k\}}$

${\displaystyle G_{k}={\begin{bmatrix}g_{k-m_{k}}&\dots &g_{k}\end{bmatrix}}}$

${\displaystyle \alpha _{k}=\operatorname {argmin} _{\alpha \in A_{k}}\|G_{k}\alpha \|_{2},\quad {\text{where}}\;A_{k}=\{\alpha =(\alpha _{0},\dots ,\alpha _{m_{k}})\in \mathbb {R} ^{m_{k}+1}:\sum _{i=0}^{m_{k}}\alpha _{i}=1\}}$

${\displaystyle x_{k+1}=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$

where the matrix–vector multiplication ${\displaystyle G_{k}\alpha =\sum _{i=0}^{m_{k}}(\alpha )_{i}g_{k-m_{k}+i}}$, and ${\displaystyle (\alpha )_{i}}$ is the ${\displaystyle i}$th element of ${\displaystyle \alpha }$. Conventional stopping criteria can be used to end the iterations of the method. For example, iterations can be stopped when ${\displaystyle \|x_{k+1}-x_{k}\|}$ falls under a prescribed tolerance, or when the residual ${\displaystyle g(x_{k})}$ falls under a prescribed tolerance.^[2]

With respect to the standard fixed-point iteration, the method has been found to converge faster and be more robust, and in some cases avoid the divergence of the fixed-point sequence.^[3][4]

For the solution ${\displaystyle x^{*}}$, we know that ${\displaystyle f(x^{*})=x^{*}}$, which is equivalent to saying that ${\displaystyle g(x^{*})={\vec {0}}}$. We can therefore rephrase the problem as an optimization problem where we want to minimize ${\displaystyle \|g(x)\|_{2}}$.

Instead of going directly from ${\displaystyle x_{k}}$ to ${\displaystyle x_{k+1}}$ by choosing ${\displaystyle x_{k+1}=f(x_{k})}$ as in fixed-point iteration, let's consider an intermediate point ${\displaystyle x'_{k+1}}$ that we choose to be the linear combination ${\displaystyle x'_{k+1}=X_{k}\alpha _{k}}$, where the coefficient vector ${\displaystyle \alpha _{k}\in A_{k}}$, and ${\displaystyle X_{k}={\begin{bmatrix}x_{k-m_{k}}&\dots &x_{k}\end{bmatrix}}}$ is the matrix containing the last ${\displaystyle m_{k}+1}$ points, and choose ${\displaystyle x'_{k+1}}$ such that it minimizes ${\displaystyle \|g(x'_{k+1})\|_{2}}$. Since the elements in ${\displaystyle \alpha _{k}}$ sum to one, we can make the first order approximation ${\displaystyle g(X_{k}\alpha _{k})=g\left(\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}\right)\approx \sum _{i=0}^{m_{k}}(\alpha _{k})_{i}g(x_{k-m_{k}+i})=G_{k}\alpha _{k}}$, and our problem becomes to find the ${\displaystyle \alpha }$ that minimizes ${\displaystyle \|G_{k}\alpha \|_{2}}$. After having found ${\displaystyle \alpha _{k}}$, we could in principle calculate ${\displaystyle x'_{k+1}}$.

However, since ${\displaystyle f}$ is designed to bring a point closer to ${\displaystyle x^{*}}$, ${\displaystyle f(x'_{k+1})}$ is probably closer to ${\displaystyle x^{*}}$ than what ${\displaystyle x'_{k+1}}$ is, so it makes sense to choose ${\displaystyle x_{k+1}=f(x'_{k+1})}$ rather than ${\displaystyle x_{k+1}=x'_{k+1}}$. Furthermore, since the elements in ${\displaystyle \alpha _{k}}$ sum to one, we can make the first order approximation ${\displaystyle f(x'_{k+1})=f\left(\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}\right)\approx \sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f(x_{k-m_{k}+i})=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$. We therefore choose

${\displaystyle x_{k+1}=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$.

At each iteration of the algorithm, the constrained optimization problem ${\displaystyle \operatorname {argmin} \|G_{k}\alpha \|_{2}}$, subject to ${\displaystyle \alpha \in A_{k}}$ needs to be solved. The problem can be recast in several equivalent formulations,^[3] yielding different solution methods which may result in a more convenient implementation:

• defining the matrices ${\displaystyle {\mathcal {G}}_{k}={\begin{bmatrix}g_{k-m_{k}+1}-g_{k-m_{k}}&\dots &g_{k}-g_{k-1}\end{bmatrix}}}$ and ${\displaystyle {\mathcal {X}}_{k}={\begin{bmatrix}x_{k-m_{k}+1}-x_{k-m_{k}}&\dots &x_{k}-x_{k-1}\end{bmatrix}}}$, solve ${\displaystyle \gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}\|g_{k}-{\mathcal {G}}_{k}\gamma \|_{2}}$, and set ${\displaystyle x_{k+1}=x_{k}+g_{k}-({\mathcal {X}}_{k}+{\mathcal {G}}_{k})\gamma _{k}}$;^[3][4]
• solve ${\displaystyle \theta _{k}=\{(\theta _{k})_{i}\}_{i=1}^{m_{k}}=\operatorname {argmin} _{\theta \in \mathbb {R} ^{m_{k}}}\left\|g_{k}+\sum _{i=1}^{m_{k}}\theta _{i}(g_{k-i}-g_{k})\right\|_{2}}$, then set ${\displaystyle x_{k+1}=x_{k}+g_{k}+\sum _{j=1}^{m_{k}}(\theta _{k})_{j}(x_{k-j}-x_{k}+g_{k-j}-g_{k})}$.^[1]

For both choices, the optimization problem is in the form of an unconstrained linear least-squares problem, which can be solved by standard methods including QR decomposition^[3] and singular value decomposition,^[4] possibly including regularization techniques to deal with rank deficiencies and conditioning issues in the optimization problem. Solving the least-squares problem by solving the normal equations is generally not advisable due to potential numerical instabilities and generally high computational cost.^[4]

Stagnation in the method (i.e. subsequent iterations with the same value, ${\displaystyle x_{k+1}=x_{k}}$) causes the method to break down, due to the singularity of the least-squares problem. Similarly, near-stagnation (${\displaystyle x_{k+1}\approx x_{k}}$) results in bad conditioning of the least squares problem. Moreover, the choice of the parameter ${\displaystyle m}$ might be relevant in determining the conditioning of the least-squares problem, as discussed below.^[3]

The algorithm can be modified introducing a variable relaxation parameter (or mixing parameter) ${\displaystyle \beta _{k}>0}$.^[1][3][4] At each step, compute the new iterate as $${\displaystyle x_{k+1}=(1-\beta _{k})\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}+\beta _{k}\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f(x_{k-m_{k}+i})\;.}$$The choice of ${\displaystyle \beta _{k}}$ is crucial to the convergence properties of the method; in principle, ${\displaystyle \beta _{k}}$ might vary at each iteration, although it is often chosen to be constant.^[4]

The parameter ${\displaystyle m}$ determines how much information from previous iterations is used to compute the new iteration ${\displaystyle x_{k+1}}$. On the one hand, if ${\displaystyle m}$ is chosen to be too small, too little information is used and convergence may be undesirably slow. On the other hand, if ${\displaystyle m}$ is too large, information from old iterations may be retained for too many subsequent iterations, so that again convergence may be slow.^[3] Moreover, the choice of ${\displaystyle m}$ affects the size of the optimization problem. A too large value of ${\displaystyle m}$ may worsen the conditioning of the least squares problem and the cost of its solution.^[3] In general, the particular problem to be solved determines the best choice of the ${\displaystyle m}$ parameter.^[3]

With respect to the algorithm described above, the choice of ${\displaystyle m_{k}}$ at each iteration can be modified. One possibility is to choose ${\displaystyle m_{k}=k}$ for each iteration ${\displaystyle k}$ (sometimes referred to as Anderson acceleration without truncation).^[3] This way, every new iteration ${\displaystyle x_{k+1}}$ is computed using all the previously computed iterations. A more sophisticated technique is based on choosing ${\displaystyle m_{k}}$ so as to maintain a small enough conditioning for the least-squares problem.^[3]

Newton's method can be applied to the solution of ${\displaystyle f(x)-x=0}$ to compute a fixed point of ${\displaystyle f(x)}$ with quadratic convergence. However, such method requires the evaluation of the exact derivative of ${\displaystyle f(x)}$, which can be very costly.^[4] Approximating the derivative by means of finite differences is a possible alternative, but it requires multiple evaluations of ${\displaystyle f(x)}$ at each iteration, which again can become very costly. Anderson acceleration requires only one evaluation of the function ${\displaystyle f(x)}$ per iteration, and no evaluation of its derivative. On the other hand, the convergence of an Anderson-accelerated fixed point sequence is still linear in general.^[5]

Several authors have pointed out similarities between the Anderson acceleration scheme and other methods for the solution of non-linear equations. In particular:

• Eyert^[6] and Fang and Saad^[4] interpreted the algorithm within the class of quasi-Newton and multisecant methods, that generalize the well known secant method, for the solution of the non-linear equation ${\displaystyle g(x)=0}$; they also showed how the scheme can be seen as a method in the Broyden class;^[7]
• Walker and Ni^[3][8] showed that the Anderson acceleration scheme is equivalent to the GMRES method in the case of linear problems (i.e. the problem of finding a solution to ${\displaystyle A\mathbf {x} =\mathbf {x} }$ for some square matrix ${\displaystyle A}$), and can thus be seen as a generalization of GMRES to the non-linear case; a similar result was found by Washio and Oosterlee.^[9]

Moreover, several equivalent or nearly equivalent methods have been independently developed by other authors,^{[9][10][11][12][13]} although most often in the context of some specific application of interest rather than as a general method for fixed point equations.

The following is an example implementation in MATLAB language of the Anderson acceleration scheme for finding the fixed-point of the function ${\displaystyle f(x)=\sin(x)+\arctan(x)}$. Notice that:

• the optimization problem was solved in the form ${\displaystyle \gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}\|g_{k}-{\mathcal {G}}_{k}\gamma \|_{2}}$ using QR decomposition;
• the computation of the QR decomposition is sub-optimal: indeed, at each iteration a single column is added to the matrix ${\displaystyle {\mathcal {G}}_{k}}$, and possibly a single column is removed; this fact can be exploited to efficiently update the QR decomposition with less computational effort;^[14]
• the algorithm can be made more memory-efficient by storing only the latest few iterations and residuals, if the whole vector of iterations ${\displaystyle x_{k}}$ is not needed;
• the code is straightforwardly generalized to the case of a vector-valued ${\displaystyle f(x)}$.

f = @(x) sin(x) + atan(x); % Function whose fixed point is to be computed.
x0 = 1; % Initial guess.

k_max = 100; % Maximum number of iterations.
tol_res = 1e-6; % Tolerance on the residual.
m = 3; % Parameter m.

x = [x0, f(x0)]; % Vector of iterates x.
g = f(x) - x; % Vector of residuals.

G_k = g(2) - g(1); % Matrix of increments in residuals.
X_k = x(2) - x(1); % Matrix of increments in x.

k = 2;
while k < k_max && abs(g(k)) > tol_res
    m_k = min(k, m);
 
    % Solve the optimization problem by QR decomposition.
    [Q, R] = qr(G_k);
    gamma_k = R \ (Q' * g(k));
 
    % Compute new iterate and new residual.
    x(k + 1) = x(k) + g(k) - (X_k + G_k) * gamma_k;
    g(k + 1) = f(x(k + 1)) - x(k + 1);
 
    % Update increment matrices with new elements.
    X_k = [X_k, x(k + 1) - x(k)];
    G_k = [G_k, g(k + 1) - g(k)];
 
    n = size(X_k, 2);
    if n > m_k
        X_k = X_k(:, n - m_k + 1:end);
        G_k = G_k(:, n - m_k + 1:end);
    end
 
    k = k + 1;
end

% Prints result: Computed fixed point 2.013444 after 9 iterations
fprintf("Computed fixed point %f after %d iterations\n", x(end), k);

• Fixed-point
• Fixed-point iteration
• Root-finding algorithms