Chebyshev iteration method

An iterative algorithm for finding a solution to a linear equation

$$ \tag{1 } A u = f $$

that takes account of information about the inclusion of $ \mathop{\rm Tr} ( A) $— the spectrum of the operator $ A $— in a certain set $ \Omega $, and uses the properties and parameters of those polynomials that deviate least from zero on $ \Omega $ and are equal to 1 at 0.

The most well-developed Chebyshev iteration method is obtained when in (1), $ A $ is a linear self-adjoint operator and $ \mathop{\rm Tr} ( A) \in [ m , M ] $, where $ 0 < m < M $ are the boundary points of the spectrum; then the Chebyshev iteration method uses the properties of the Chebyshev polynomials of the first kind, $ T _ {n} ( x) $. For this case one considers two types of Chebyshev iteration methods:

$$ \tag{2 } u ^ {k+ 1} = u ^ {k} - \alpha _ {k+ 1} ( A u ^ {k} - f ) , $$

$$ \tag{3 } u ^ {k+ 1} = u ^ {k} - \alpha _ {k+ 1} ( A u ^ {k} - f ) - \beta _ {k+ 1} ( u ^ {k} - u ^ {k- 1} ) , $$

$$ \beta _ {1} = 0,\ k = 0 , $$

in which for a given $ u ^ {0} $ one obtains a sequence $ u ^ {k} \rightarrow u $ as $ k \rightarrow \infty $. In (2) and (3) $ \alpha _ {i} $ and $ \beta _ {i} $ are the numerical parameters of the method. If $ \epsilon ^ {k} = u - u ^ {k} $, then the initial error $ \epsilon ^ {0} $ and the error at the $ N $-th iteration $ \epsilon ^ {N} $ are related by the formula

$$ \epsilon ^ {N} = P _ {N} ( A) \epsilon ^ {0} , $$

where

$$ \tag{4 } P _ {N} ( t) = \prod _ { i= 1} ^ { N } ( 1 - \gamma _ {i} t ) ,\ \ P _ {N} ( 0) = 1 . $$

The polynomials $ P _ {N} ( t) $ are calculated using the parameters of each of the methods (2), (3): for method (2)

$$ \tag{5 } \alpha _ {k} = \gamma _ {j _ {k} } ,\ \ k = 1, \dots, N , $$

where $ 1 \leq j _ {k} \leq N $ are the elements of the permutation $ \kappa _ {N} = ( j _ {1}, \dots, j _ {N} ) $, while for method (3) they are calculated from the recurrence relations

$$ \tag{6 } P _ {k+ 1} ( t) = ( 1 - \beta _ {k+ 1} - \alpha _ {k+ 1} t ) P _ {k} ( t) + \beta _ {k+ 1} P _ {k- 1} ( t) , $$

$$ P _ {0} = 1 ,\ \beta _ {1} = 0 ,\ k = 0, \dots, N - 1 . $$

Here

$$ \| \epsilon ^ {N} \| \leq \sup _ {t \in [ m , M ] } \ | P _ {N} ( t) | \cdot \| \epsilon ^ {0} \| . $$

The methods (2) and (3) can be optimized on the class of problems for which $ \mathop{\rm Tr} ( A) \in [ m , M ] $ by choosing the parameters such that $ P _ {N} ( t) $ in (4) is the polynomial least deviating from zero on $ [ m , M ] $. It was proved in 1881 by P.L. Chebyshev that this is the polynomial

$$ \tag{7 } P _ {N} ( t) = \ \frac{T _ {N} ( {( M + m - 2 t) } / ( M- m) ) }{T _ {N} ( \theta ) } , $$

where $ \theta = ( M + m ) / ( M - m ) $. Then

$$ \tag{8 } \| \epsilon ^ {N} \| \leq \ \frac{2 \tau ^ {N} }{( 1 + \tau ^ {2N} ) } \| \epsilon ^ {0} \| , $$

where

$$ \tau = \frac{1 - \sqrt {m / M } }{1 + \sqrt {m / M } } . $$

Substituting (7) for $ N = k - 1 , k , k + 1 $ in (6), the parameters $ \alpha _ {k+ 1} , \beta _ {k+ 1} $ of the method (3) are determined:

$$ \tag{9 } \alpha _ {k+ 1} = \frac{4 \delta _ {k+ 1} }{M - m } ,\ \beta _ {k+ 1} = - \delta _ {k} \delta _ {k+ 1} , $$

where

$$ \tag{10 } \delta _ {0} = 0 ,\ \delta _ {1} = \theta ^ {- 1} ,\ \ \delta _ {k+ 1} = ( 2 \theta - \delta _ {k} ) ^ {- 1} , $$

$$ k = 1, \dots, N - 1 . $$

Thus, computing $ \alpha _ {k+ 1} $ and $ \beta _ {k+ 1} $ by the formulas (9) and (10), one obtains the Chebyshev iteration method (3) for which $ \| \epsilon ^ {N} \| $ is optimally small for each $ N \geq 1 $.

To optimize (2) for a given $ N $, the parameters $ \alpha _ {k+ 1} $ are chosen corresponding to the permutation $ \kappa _ {N} $ in formula (5) in such a way that (7) holds, that is,

$$ \tag{11 } \gamma _ {j} = 2 ( M + m - ( M - m ) \cos \pi \psi _ {j} )^{-1} , $$

$$ \psi _ {j} = \frac{2 j - 1 }{2N} ,\ j = 1, \dots, N . $$

Then after $ N $ iterations, inequality (8) holds for $ \| \epsilon ^ {N} \| $.

An important problem for small $ m / M $ is the question of the stability of the method (2), (5), (11). An imprudent choice of $ \kappa _ {N} $ may lead to a catastrophic increase in $ \| u ^ {k} \| $ for some $ 1 \leq k \leq N $, to the loss of significant figures, or to an increase in the rounding-off errors allowed on intermediate iteration. There exist algorithms that mix the parameters in (11) and guarantee the stability of the calculations: for $ N = 2 ^ {p} $ see Iteration algorithm; and for $ N = 3 ^ {p} $ one of the algorithms for constructing $ \kappa _ {N} $ is as follows. Let $ \kappa _ {1} = ( 1) $, and suppose that $ \kappa _ {3 ^ {r- 1} } = ( j _ {1}, \dots, j _ {3 ^ {r- 1} } ) $ has been constructed, then

$$ \tag{12 } \kappa _ {3 ^ {r} } = ( j _ {1} , 2 \cdot 3 ^ {r- 1} + j _ {1} ,\ 2 \cdot 3 ^ {r- 1} + 1 - j _ {1}, \dots $$

$$ \dots, {} 2 \cdot 3 ^ {r- 1} + 1 - j _ {3 ^ {r- 1} } ) ,\ r = 1, \dots, p . $$

There exists a class of methods (2) — the stable infinitely repeated optimal Chebyshev iteration methods — that allows one to repeat the method (2), (5), (11) after $ N $ iterations in such a way that it is stable and such that it becomes optimal again for some sequence $ N _ {i} \rightarrow \infty $. For the case $ N _ {i} = 3 ^ {i} N $, it is clear from the formula

$$ \tag{13 } T _ {3N} ( x) = T _ {N} ( x) ( 2 T _ {N} ( x) + \sqrt 3 ) ( 2 T _ {N} ( x) - \sqrt 3 ) $$

that $ P _ {3N} ( t) $ agrees with (11). If after $ N $ iterations one repeats the iteration (2), (5), (11) further, taking for $ \psi _ {j} $ in (11) the $ 2N $ values

$$ \tag{14 } \widetilde \psi _ {j} = \frac{2 j - 1 }{6 N } ,\ 2j \not\equiv 1 ( \mathop{\rm mod} 3 ) , $$

then once again one obtains a Chebyshev iteration method after $ 3 N $ iterations. To ensure stability, the set (14) is decomposed into two sets: in the $ i $- th set, $ i = 1 , 2 $, one puts the $ \widetilde \psi _ {j} $ for which $ \pi \cos \widetilde \psi _ {j} $ is a root of the $ i $-th bracket in (13); within each of the subsets the $ \widetilde \psi _ {j} $ are permuted according to the permutation $ \kappa _ {N} $. For $ N < k < 2 N $ one substitutes elements of the first set in (5), (11), and for $ 2 N < k \leq 3 N $ one uses the second subset; the permutation $ \kappa _ {3N} $ is defined in the same way. Continuing in an analogous way the process of forming parameters, one obtains an infinite sequence $ \{ \omega _ {j} \} _ {1} ^ \infty $, uniformly distributed on $ [ 0 , 1 ] $, called a $ T $- sequence, for which the method (2) becomes optimal with $ N _ {i} = 3 ^ {i} N $ and

$$ \tag{15 } \alpha _ {k+ 1} = 2 ( M + m - ( M - m ) \cos \pi \omega _ {k+ 1} ) ^ {- 1} , $$

$$ k = 0 , 1 , . . . . $$

The theory of the Chebyshev iteration methods (2), (3) can be extended to partial eigen value problems. Generalizations also exist to a certain class of non-self-adjoint operators, when $ \mathop{\rm Tr} ( A) $ lies in a certain interval or within a certain domain of special shape (in particular, an ellipse); when information is known about the distribution of the initial error; or when the Chebyshev iteration method is combined with the method of conjugate gradients.

One of the effective methods of speeding up to the convergence of the iterations (2), (3) is a preliminary transformation of equation (1) to an equivalent equation of the form

$$ B A u = B f , $$

and the application of the Chebyshev iteration method to this equation. The operator $ B $ is defined by taking account of two facts: 1) the algorithm for computing a quantity of the form $ Bv $ should not be laborious; and 2) $ \mathop{\rm Tr} ( B A ) $ should lie in a set that ensures the fast convergence of the Chebyshev iteration method.

References

[1]	G.I. Marchuk, V.I. Lebedev, "Numerical methods in the theory of neutron transport" , Harwood (1986) (Translated from Russian)
[2]	N.S. Bakhvalov, "Numerical methods: analysis, algebra, ordinary differential equations" , MIR (1977) (Translated from Russian)
[3]	G.I. Marchuk, "Methods of numerical mathematics" , Springer (1982) (Translated from Russian)
[4]	A.A. Samarskii, "Theorie der Differenzverfahren" , Akad. Verlagsgesell. Geest u. Portig K.-D. (1984) (Translated from Russian)
[5a]	V.I. Lebedev, S.A. Finogenov, "The order of choices of the iteration parameters in the cyclic Chebyshev iteration method" Zh. Vychisl. Mat. i Mat. Fiz. , 11 : 2 (1971) pp. 425–438 (In Russian)
[5b]	V.I. Lebedev, S.A. Finogenov, "Solution of the problem of parameter ordering in Chebyshev iteration methods" Zh. Vychisl. Mat. i Mat. Fiz , 13 : 1 (1973) pp. 18–33 (In Russian)
[5c]	V.I. Lebedev, S.A. Finogenov, "The use of ordered Chebyshev parameters in iteration methods" Zh. Vychisl. Mat. i Mat. Fiz. , 16 : 4 (1976) pp. 895–907 (In Russian)
[6a]	V.I. Lebedev, "Iterative methods for solving operator equations with spectrum located on several segments" Zh. Vychisl. Mat. i Mat. Fiz. , 9 : 6 (1969) pp. 1247–1252 (In Russian)
[6b]	V.I. Lebedev, "Iteration methods for solving linear operator equations, and polynomials deviating least from zero" , Mathematical analysis and related problems in mathematics , Novosibirsk (1978) pp. 89–108 (In Russian)

Comments

In the Western literature the method (2), (5), (11) is known as the Richardson method of first degree [a2] or, more widely used, the Chebyshev semi-iterative method of first degree. The method goes back to an early paper of L.F. Richardson , where the method (2), (5) was already proposed. However, Richardson did not identify the zeros $ 1/ \alpha _ {k} $ of $ P _ {N} ( t) $ with the zeros of (shifted) Chebyshev polynomials as done in (11), but (less sophisticatedly) sprinkled them uniformly over the interval $ [ m, M] $. The use of Chebyshev polynomials seems to be proposed for the first time in [a1] and [a3].

The "stable infinitely repeated optimal Chebyshev iteration methods" outlined above are based on the identity $ T _ {pq} ( x) \equiv T _ {p} ( T _ {q} ( x)) $, which immediately leads to the factorization

$$ T _ {pq} ( x) = \ \prod _ {j = 1 } ^ { p } \frac{T _ {q} ( x) - \cos (( 2j - 1) \pi / 2p) }{1 - \cos (( 2j - 1) \pi / 2p) } . $$

This formula has already been used in [a1] in the numerical determination of fundamental modes.

The method (3), (9) is known as Richardson's method or Chebyshev's semi-iterative method of second degree. It was suggested in [a9] and turns out to be completely stable; thus, at the cost of an extra storage array the instability problems associated with the first-degree process are avoided.

As to the choice of the transformation operator $ B $( called "preconditioningpreconditioning" ), an often used "preconditionerpreconditioner" is the so-called SSOR matrix (Symmetric Successive Over-Relaxation matrix) proposed in [a8].

Introductions to the theory of Chebyshev semi-iterative methods are provided by [a2] and [a3]. An extensive analysis can be found in [a10], Chapt. 5 and in [a4]. In this work the spectrum of the operator $ A $ is assumed to be real. An analysis of the case where the spectrum is not real can be found in [a5].

Instead of using minimax polynomials, one may consider integral measures for "minimizing" $ P _ {N} ( t) $ on $ [ m, M] $. This leads to the theory of kernel polynomials introduced in [a9] and extended in [a11], Chapt. 5.

Iterative methods as opposed to direct methods (cf. Direct method) only make sense when the matrix is sparse (cf. Sparse matrix). Moreover, their versatility depends on how large an error $ ( \epsilon _ {N} ) $ is tolerated; often other errors, e.g., truncation errors in discretized systems of partial differential equations, are more dominant.

When no information about the eigen structure of $ A $ is available, or in the non-self-adjoint case, it is often preferable to use the method of conjugate gradients (cf. Conjugate gradients, method of). Numerical algorithms based on the latter method combined with incomplete factorization have proven to be one of the most efficient ways to solve linear problems up to now (1987).

References

[a1]	D.A. Flanders, G. Shortley, "Numerical determination of fundamental modes" J. Appl. Physics , 21 (1950) pp. 1326–1332
[a2]	G.E. Forsythe, W.R. Wasow, "Finite difference methods for partial differential equations" , Wiley (1960)
[a3]	G.H. Golub, C.F. van Loan, "Matrix computations" , North Oxford Acad. (1983)
[a4]	G.H. Golub, R.S. Varga, "Chebyshev semi-iterative methods, successive over-relaxation methods and second-order Richardson iterative methods I, II" Num. Math. , 3 (1961) pp. 147–156; 157–168
[a5]	T.A. Manteuffel, "The Tchebychev iteration for nonsymmetric linear systems" Num. Math. , 28 (1977) pp. 307–327
[a6a]	L.F. Richardson, "The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam" Philos. Trans. Roy. Soc. London Ser. A , 210 (1910) pp. 307–357
[a6b]	L.F. Richardson, "The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam" Proc. Roy. Soc. London Ser. A , 83 (1910) pp. 335–336
[a7]	G. Shortley, "Use of Tchebycheff-polynomial operators in the numerical solution of boundary-value problems" J. Appl. Physics , 24 (1953) pp. 392–396
[a8]	J.W. Sheldon, "On the numerical solution of elliptic difference equations" Math. Tables Aids Comp. , 9 (1955) pp. 101–112
[a9]	E.L. Stiefel, "Kernel polynomials in linear algebra and their numerical applications" , Appl. Math. Ser. , 49 , Nat. Bur. Standards (1958)
[a10]	R.S. Varga, "Matrix iterative analysis" , Prentice-Hall (1962)
[a11]	E.L. Wachspress, "Iterative solution of elliptic systems, and applications to the neutron diffusion equations of nuclear physics" , Prentice-Hall (1966)