Anderson acceleration

From Wikipedia, the free encyclopedia

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 often arising in the field of computational science.


Given a function , consider the problem of finding a fixed point of , which is a solution to the equation . A classical approach to the problem is to employ a fixed-point iteration scheme;[2] that is, given an initial guess for the solution, to compute the sequence 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 is computationally expensive.[2] Anderson acceleration is a method to accelerate the convergence of the fixed-point sequence.[2]

Define the residual , and denote and (where corresponds to the sequence of iterates from the previous paragraph). Given an initial guess and an integer parameter , the method can be formulated as follows:[3][note 1]

where the matrix–vector multiplication , and is the th element of . Conventional stopping criteria can be used to end the iterations of the method. For example, iterations can be stopped when falls under a prescribed tolerance, or when the residual 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 , we know that , which is equivalent to saying that . We can therefore rephrase the problem as an optimization problem where we want to minimize .

Instead of going directly from to by choosing as in fixed-point iteration, let's consider an intermediate point that we choose to be the linear combination , where the coefficient vector , and is the matrix containing the last points, and choose such that it minimizes . Since the elements in sum to one, we can make the first order approximation , and our problem becomes to find the that minimizes . After having found , we could in principle calculate .

However, since is designed to bring a point closer to , is probably closer to than what is, so it makes sense to choose rather than . Furthermore, since the elements in sum to one, we can make the first order approximation . We therefore choose


Solution of the minimization problem[edit]

At each iteration of the algorithm, the constrained optimization problem , subject to 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 and , solve , and set ;[3][4]
  • solve , then set .[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, ) causes the method to break down, due to the singularity of the least-squares problem. Similarly, near-stagnation () results in bad conditioning of the least squares problem. Moreover, the choice of the parameter 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) .[1][3][4] At each step, compute the new iterate as

The choice of is crucial to the convergence properties of the method; in principle, might vary at each iteration, although it is often chosen to be constant.[4]

Choice of m[edit]

The parameter determines how much information from previous iterations is used to compute the new iteration . On the one hand, if is chosen to be too small, too little information is used and convergence may be undesirably slow. On the other hand, if 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 affects the size of the optimization problem. A too large value of 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 parameter.[3]

Choice of mk[edit]

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

Relations to other classes of methods[edit]

Newton's method can be applied to the solution of to compute a fixed point of with quadratic convergence. However, such method requires the evaluation of the exact derivative of , which can be very costly.[4] Approximating the derivative by means of finite differences is a possible alternative, but it requires multiple evaluations of at each iteration, which again can become very costly. Anderson acceleration requires only one evaluation of the function 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 ; 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 for some square matrix ), 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.

Example MATLAB implementation[edit]

The following is an example implementation in MATLAB language of the Anderson acceleration scheme for finding the fixed-point of the function . Notice that:

  • the optimization problem was solved in the form using QR decomposition;
  • the computation of the QR decomposition is sub-optimal: indeed, at each iteration a single column is added to the matrix , 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 is not needed;
  • the code is straightforwardly generalized to the case of a vector-valued .
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);
    k = k + 1;

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

See also[edit]


  1. ^ This formulation is not the same as given by the original author;[1] it is an equivalent, more explicit formulation given by Walker and Ni.[3]


  1. ^ a b c d Anderson, Donald G. (October 1965). "Iterative Procedures for Nonlinear Integral Equations". Journal of the ACM. 12 (4): 547–560. doi:10.1145/321296.321305.
  2. ^ a b c d Quarteroni, Alfio; Sacco, Riccardo; Saleri, Fausto. Numerical mathematics (2nd ed.). Springer. ISBN 978-3-540-49809-4.
  3. ^ a b c d e f g h i j k l m n Walker, Homer F.; Ni, Peng (January 2011). "Anderson Acceleration for Fixed-Point Iterations". SIAM Journal on Numerical Analysis. 49 (4): 1715–1735. doi:10.1137/10078356X.
  4. ^ a b c d e f g h Fang, Haw-ren; Saad, Yousef (March 2009). "Two classes of multisecant methods for nonlinear acceleration". Numerical Linear Algebra with Applications. 16 (3): 197–221. doi:10.1002/nla.617.
  5. ^ Evans, Claire; Pollock, Sara; Rebholz, Leo G.; Xiao, Mengying (20 February 2020). "A Proof That Anderson Acceleration Improves the Convergence Rate in Linearly Converging Fixed-Point Methods (But Not in Those Converging Quadratically)". SIAM Journal on Numerical Analysis. 58 (1): 788–810. arXiv:1810.08455. doi:10.1137/19M1245384.
  6. ^ Eyert, V. (March 1996). "A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences". Journal of Computational Physics. 124 (2): 271–285. doi:10.1006/jcph.1996.0059.
  7. ^ Broyden, C. G. (1965). "A class of methods for solving nonlinear simultaneous equations". Mathematics of Computation. 19 (92): 577–577. doi:10.1090/S0025-5718-1965-0198670-6.
  8. ^ Ni, Peng (November 2009). Anderson Acceleration of Fixed-point Iteration with Applications to Electronic Structure Computations (PhD).
  9. ^ a b Oosterlee, C. W.; Washio, T. (January 2000). "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows". SIAM Journal on Scientific Computing. 21 (5): 1670–1690. doi:10.1137/S1064827598338093.
  10. ^ Pulay, Péter (July 1980). "Convergence acceleration of iterative sequences. the case of scf iteration". Chemical Physics Letters. 73 (2): 393–398. doi:10.1016/0009-2614(80)80396-4.
  11. ^ Pulay, P. (1982). "ImprovedSCF convergence acceleration". Journal of Computational Chemistry. 3 (4): 556–560. doi:10.1002/jcc.540030413.
  12. ^ Carlson, Neil N.; Miller, Keith (May 1998). "Design and Application of a Gradient-Weighted Moving Finite Element Code I: in One Dimension". SIAM Journal on Scientific Computing. 19 (3): 728–765. doi:10.1137/S106482759426955X.
  13. ^ Miller, Keith (November 2005). "Nonlinear Krylov and moving nodes in the method of lines". Journal of Computational and Applied Mathematics. 183 (2): 275–287. doi:10.1016/
  14. ^ Daniel, J. W.; Gragg, W. B.; Kaufman, L.; Stewart, G. W. (October 1976). "Reorthogonalization and stable algorithms for updating the Gram-Schmidt $QR$ factorization". Mathematics of Computation. 30 (136): 772–772. doi:10.1090/S0025-5718-1976-0431641-8.