Anderson acceleration explained

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

f(x)=x

often arising in the field of computational science.

Definition

Given a function

f:Rn\toRn

, consider the problem of finding a fixed point of

f

, which is a solution to the equation

f(x)=x

. A classical approach to the problem is to employ a fixed-point iteration scheme;[2] that is, given an initial guess

x0

for the solution, to compute the sequence

xi+1=f(xi)

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

f

is computationally expensive. Anderson acceleration is a method to accelerate the convergence of the fixed-point sequence.

Define the residual

g(x)=f(x)-x

, and denote

fk=f(xk)

and

gk=g(xk)

(where

xk

corresponds to the sequence of iterates from the previous paragraph). Given an initial guess

x0

and an integer parameter

m\geq1

, the method can be formulated as follows:[3]

x1=f(x0)

\forallk=1,2,...

mk=min\{m,k\}

Gk=\begin{bmatrix}

g
k-mk

&...&gk\end{bmatrix}

\alphak=

\operatorname{argmin}
\alpha\inAk

\|Gk\alpha\|2,whereAk=\{\alpha=(\alpha0,...,

\alpha
mk

)\in

mk+1
R

:

mk
\sum
i=0

\alphai=1\}

xk+1=

mk
\sum
i=0

(\alphak)i

f
k-mk+i

where the matrix–vector multiplication

Gk\alpha=

mk
\sum
i=0

(\alpha)i

g
k-mk+i
, and

(\alpha)i

is the

i

th element of

\alpha

. Conventional stopping criteria can be used to end the iterations of the method. For example, iterations can be stopped when

\|xk+1-xk\|

falls under a prescribed tolerance, or when the residual

g(xk)

falls under a prescribed tolerance.

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.

Derivation

For the solution

x*

, we know that

f(x*)=x*

, which is equivalent to saying that

g(x*)=\vec{0}

. We can therefore rephrase the problem as an optimization problem where we want to minimize

\|g(x)\|2

.

Instead of going directly from

xk

to

xk+1

by choosing

xk+1=f(xk)

as in fixed-point iteration, let's consider an intermediate point

x'k+1

that we choose to be the linear combination

x'k+1=Xk\alphak

, where the coefficient vector

\alphak\inAk

, and

Xk=\begin{bmatrix}

x
k-mk

&...&xk\end{bmatrix}

is the matrix containing the last

mk+1

points, and choose

x'k+1

such that it minimizes

\|g(x'k+1)\|2

. Since the elements in

\alphak

sum to one, we can make the first order approximation

g(Xk\alphak)=

mk
g\left(\sum
i=0

(\alphak)i

x
k-mk+i

\right)

mk
\sum
i=0

(\alphak)i

g(x
k-mk+i

)=Gk\alphak

, and our problem becomes to find the

\alpha

that minimizes

\|Gk\alpha\|2

. After having found

\alphak

, we could in principle calculate

x'k+1

.

However, since

f

is designed to bring a point closer to

x*

,

f(x'k+1)

is probably closer to

x*

than what

x'k+1

is, so it makes sense to choose

xk+1=f(x'k+1)

rather than

xk+1=x'k+1

. Furthermore, since the elements in

\alphak

sum to one, we can make the first order approximation

f(x'k+1)=

mk
f\left(\sum
i=0

(\alphak)i

x
k-mk+i

\right)

mk
\sum
i=0

(\alphak)i

f(x
k-mk+i

)=

mk
\sum
i=0

(\alphak)i

f
k-mk+i
. We therefore choose

xk+1=

mk
\sum
i=0

(\alphak)i

f
k-mk+i
.

Solution of the minimization problem

At each iteration of the algorithm, the constrained optimization problem

\operatorname{argmin}\|Gk\alpha\|2

, subject to

\alpha\inAk

needs to be solved. The problem can be recast in several equivalent formulations, yielding different solution methods which may result in a more convenient implementation:

l{G}k=\begin{bmatrix}

g
k-mk+1

-

g
k-mk

&...&gk-gk-1\end{bmatrix}

and

l{X}k=\begin{bmatrix}

x
k-mk+1

-

x
k-mk

&...&xk-xk-1\end{bmatrix}

, solve

\gammak=

\operatorname{argmin}
mk
\gamma\inR

\|gk-l{G}k\gamma\|2

, and set

xk+1=xk+gk-(l{X}k+l{G}k)\gammak

;[4]

\thetak=\{(\thetak)i\}

mk
i=1

=

\operatorname{argmin}
mk
\theta\inR

\left\|gk+

mk
\sum
i=1

\thetai(gk-i-gk)\right\|2

, then set

xk+1=xk+gk+

mk
\sum
j=1

(\thetak)j(xk-j-xk+gk-j-gk)

.

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 and singular value decomposition, 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.

Stagnation in the method (i.e. subsequent iterations with the same value,

xk+1=xk

) causes the method to break down, due to the singularity of the least-squares problem. Similarly, near-stagnation (

xk+1xk

) results in bad conditioning of the least squares problem. Moreover, the choice of the parameter

m

might be relevant in determining the conditioning of the least-squares problem, as discussed below.

Relaxation

The algorithm can be modified introducing a variable relaxation parameter (or mixing parameter)

\betak>0

. At each step, compute the new iterate as x_ = (1 - \beta_k)\sum_^(\alpha_k)_i x_ + \beta_k \sum_^(\alpha_k)_i f(x_)\;.The choice of

\betak

is crucial to the convergence properties of the method; in principle,

\betak

might vary at each iteration, although it is often chosen to be constant.

Choice of

The parameter

m

determines how much information from previous iterations is used to compute the new iteration

xk+1

. On the one hand, if

m

is chosen to be too small, too little information is used and convergence may be undesirably slow. On the other hand, if

m

is too large, information from old iterations may be retained for too many subsequent iterations, so that again convergence may be slow. Moreover, the choice of

m

affects the size of the optimization problem. A too large value of

m

may worsen the conditioning of the least squares problem and the cost of its solution. In general, the particular problem to be solved determines the best choice of the

m

parameter.

Choice of

With respect to the algorithm described above, the choice of

mk

at each iteration can be modified. One possibility is to choose

mk=k

for each iteration

k

(sometimes referred to as Anderson acceleration without truncation). This way, every new iteration

xk+1

is computed using all the previously computed iterations. A more sophisticated technique is based on choosing

mk

so as to maintain a small enough conditioning for the least-squares problem.

Relations to other classes of methods

Newton's method can be applied to the solution of

f(x)-x=0

to compute a fixed point of

f(x)

with quadratic convergence. However, such method requires the evaluation of the exact derivative of

f(x)

, which can be very costly. Approximating the derivative by means of finite differences is a possible alternative, but it requires multiple evaluations of

f(x)

at each iteration, which again can become very costly. Anderson acceleration requires only one evaluation of the function

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:

g(x)=0

; they also showed how the scheme can be seen as a method in the Broyden class;[7]

Ax=x

for some square matrix

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,[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

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

f(x)=\sin(x)+\arctan(x)

. Notice that:

\gammak=

\operatorname{argmin}
mk
\gamma\inR

\|gk-l{G}k\gamma\|2

using QR decomposition;

l{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]

xk

is not needed;

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 iterationsfprintf("Computed fixed point %f after %d iterations\n", x(end), k);

See also

Notes and References

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