View source for Acceleration methods

The theory and techniques for increasing the speed of convergence in the iterative solution of linear systems <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101601.png" />. See [[Acceleration of convergence|Acceleration of convergence]] and the references given there.

The fundamental problem of [[Linear-algebra(2)|linear algebra]] is the solution of finite-dimensional linear algebraic systems

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101602.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a1)</td></tr></table>

in which <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101603.png" /> is a given <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101604.png" />-matrix. For simplicity, let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101605.png" /> be an <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101606.png" />-matrix with real entries <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101607.png" />. The numerical solution of very large systems (e.g., <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101608.png" />) is not uncommon today (1996). Commonly, the saving grace of such large systems is their sparseness: most <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a1101609.png" /> are zero (cf. [[Sparse matrix|Sparse matrix]]). The classical direct method of solution, Gauss elimination, has at least three disadvantages: the pivoting operation fills in zeros, the operation count is <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016010.png" />, and the storage requirements when building a classical triangular factorization of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016011.png" /> increases superlinearly with <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016012.png" />.

For these reasons, indirect iterative methods for the solution of (a1) have recently come to the fore. Iterative methods may be roughly broken into two classes: splitting methods and gradient methods. Splitting methods decompose <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016013.png" /> according to

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016014.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a2)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016015.png" /> is non-singular (cf. [[Non-singular matrix|Non-singular matrix]]). The solution to (a1) is then sought as the fixed point <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016016.png" /> of the equation

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016017.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a3)</td></tr></table>

Accordingly, one picks an initial guess <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016018.png" /> and forms the iterates

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016019.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a4)</td></tr></table>

Provided that the [[Spectral radius|spectral radius]] <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016020.png" /> is less than one, the iteration (a4) converges to the solution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016021.png" /> of (a1) for any initial guess <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016022.png" />.

Gradient methods, on the other hand, form a variational problem from (a1). For example, minimizing the quadratic form

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016023.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a5)</td></tr></table>

for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016024.png" /> a symmetric positive-definite matrix, is equivalent to the problem (a1).

Just as for the classical direct methods of solution of (a1), the indirect iterative methods in their original formulations are too slow for ever larger matrices <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016025.png" />. Therefore the original Jacobi iterative method (cf. [[Jacobi method|Jacobi method]])

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016026.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a6)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016027.png" /> has been written in terms of its lower-triangular, diagonal, and upper-triangular parts, has been accelerated by the Gauss–Seidel method

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016028.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a7)</td></tr></table>

which in turn has been accelerated by the successive overrelaxation method

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016029.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a8)</td></tr></table>

(Cf. also [[Hyperrelaxation method|Hyperrelaxation method]].) In (a8), <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016030.png" /> is a free parameter which may be sought optimally. Generally, the error of a splitting iterative method is reduced by each iteration by a factor bounded by <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016031.png" />.

In like vein, the original gradient iterative method, steepest descent, in which the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016032.png" />-th residual error vector

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016033.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a9)</td></tr></table>

being the negative of the gradient of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016034.png" /> hence defines the direction of steepest descent, converges more slowly as one nears the bottom of the quadratic surface (cf. also [[Steepest descent, method of|Steepest descent, method of]]). The convergence also becomes slow when the condition number

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016035.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a10)</td></tr></table>

of the matrix <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016036.png" /> is large. Here <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016037.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016038.png" /> denote the largest and smallest eigenvalues of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016039.png" /> for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016040.png" /> symmetric positive-definite, and may be replaced in (a10) by their absolute values or, more generally, by <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016041.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016042.png" /> for general matrices <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016043.png" />. Therefore, the steepest descent method has been accelerated by the conjugate-gradient method. Both steepest descent and conjugate gradient may be written as iterations

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016044.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a11)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016045.png" /> denotes an optimal magnitude of change in the direction <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016046.png" />. In steepest descent, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016047.png" /> is the residual <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016048.png" /> of (a9) and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016049.png" /> is then chosen to minimize <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016050.png" />. In conjugate gradient one takes <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016051.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016052.png" /> is chosen so that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016053.png" /> is conjugate to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016054.png" />, i.e.,

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016055.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a12)</td></tr></table>

Then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016056.png" /> is chosen to minimize <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016057.png" /> as before. Because the conjugate-gradient method may be seen to choose a new dimension each iteration, in theory (i.e., in the absence of roundoff errors) it must converge in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016058.png" /> steps. In practice it usually converges even faster than that.

The above acceleration methods, i.e., those based upon splittings and those based upon gradients, respectively, can be referred to as the classical acceleration methods: relaxation methods and direction methods, respectively. Recently (1996) these methods have been sharpened in a number of ways. Two of the most important sharpenings are Chebyshev acceleration and preconditioning acceleration. Both methods often involve polynomial iterations

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016059.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a13)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016060.png" /> is some polynomial chosen from some knowledge or estimation of the eigenvalues of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016061.png" />. These methods therefore make use of results and techniques from both [[Approximation theory|approximation theory]] and [[Spectral theory|spectral theory]].

Chebyshev methods depend upon the property of the [[Chebyshev polynomials|Chebyshev polynomials]] of minimizing the maximum error <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016062.png" /> for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016063.png" /> in some interval <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016064.png" />. For example, for the conjugate-gradient method it may be shown that (see [[#References|[a5]]])

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016065.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a14)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016066.png" /> is the Chebyshev polynomial of the first kind of degree <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016067.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016068.png" />, and where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016069.png" />. A disadvantage of Chebyshev methods in practice is their dependence upon knowledge or good estimates of the eigenvalues of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016070.png" />.

Preconditioning has recently become the most important acceleration method. There are several preconditioning strategies, but the basic idea is to left or right multiply <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016071.png" /> in (a1) by some preconditioning matrix so that the transformed system <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016072.png" /> is easier to solve. For example, if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016073.png" /> approximates <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016074.png" /> in some correct sense of representing <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016075.png" />'s main invertibility features, i.e., so that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016076.png" /> approximates <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016077.png" /> in some good way, then (a1) is transformed by right preconditioning to the equivalent system

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016078.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a15)</td></tr></table>

Left preconditioning yields

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016079.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a16)</td></tr></table>

Two-sided preconditioning by multiplicatively splitting <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016080.png" /> into <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016081.png" /> becomes

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016082.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a17)</td></tr></table>

Because the convergence rate of iterative methods generally depends on the condition number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016083.png" /> of (a10), a rule of thumb has been to seek preconditioners <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016084.png" /> such that the transformed system <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016085.png" /> has lower, hopefully much lower, condition number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016086.png" />.

There are roughly three generally successful theories for finding preconditioners: incomplete factorizations, Frobenius-norm methods and relaxation methods. The incomplete factorizations approximate <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016087.png" /> by its approximate factorization into lower and upper triangular factors, but insisting on sparse approximations by dropping small elements. Then a conjugate-gradient method or Krylov-space variant of it (e.g., schemes such as ORTHOMIN or GMRES) may converge well. The Frobenius-norm methods, roughly, seek a weight matrix <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016088.png" /> such that the trace <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016089.png" /> is minimized. Often, an incomplete Cholesky factorization enters into applications of this method. Relaxation methods have been described above. Often, these are combined (see [[#References|[a1]]]) with block Schur complement decompositions of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016090.png" />. Relaxation methods such as (a6), (a7), (a8) are used today chiefly as preconditioners. Even the lowly Jacobi method (a6) can be a practical preconditioner on highly parallel computers when appropriately implemented.

Acceleration methods are also used for the iterative solution of the eigenvalue-eigenvector problem

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016091.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a18)</td></tr></table>

(a18) may be regarded as a special case of (a1) in which <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016092.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016093.png" />. However, one must approximate the eigenvalue <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016094.png" /> as one is approximating the eigenvector <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016095.png" />. Both Chebyshev and preconditioning acceleration methods enter into these (e.g., Arnoldi and Krylov) methods for eigenvalue-eigenvector computation. An important preconditioning method somewhat special to these problems is the so-called shift and invert method, using inverted shifts <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016096.png" /> to more widely disperse the eigenvalues (see [[#References|[a4]]]).

There are other useful variations on preconditioning for accelerating convergence. Conjugate gradient can be applied as a preconditioner for any of the relaxation schemes. Other methods, such as multigrid schemes, have been successfully used as preconditioners (see [[#References|[a3]]]). Preconditioning strategies can be cast in terms of a new theory of anti-eigenvalues <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016097.png" /> rather than the eigenvalue condition number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/a/a110/a110160/a11016098.png" /> (see [[#References|[a2]]] and [[Anti-eigenvalue|Anti-eigenvalue]]). Finding optimal preconditioners for an important specific application often usefully involves specific physical details of that application. For this reason, in practice preconditioning is sometimes labeled  "black magic" .

====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  O. Axelsson,   "Iterative solution methods" , Cambridge , New York  (1994)</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  K. Gustafson,   "Lectures on computational fluid dynamics, mathematical physics, and linear algebra" , Kaigai &amp; World Sci.  (1996/7)</TD></TR><TR><TD valign="top">[a3]</TD> <TD valign="top">  W. Hackbusch,   "Iterative solution of large sparse systems of equations" , Springer  (1994)</TD></TR><TR><TD valign="top">[a4]</TD> <TD valign="top">  Y. Saad,   "Numerical methods for large eigenvalue problems" , Halsted  (1992)</TD></TR><TR><TD valign="top">[a5]</TD> <TD valign="top">  Y. Saad,   "Iterative methods for sparse linear systems" , PWS Publishing  (1996)</TD></TR></table>