NUMERICAL STABILITY OF THE CHEBYSHEV METHOD FOR THE SOLUTION OF LARGE LINEAR SYSTEMS

This paper contains the rounding error analysis for the Chebyshev method for the solution of large linear systems Ax+g = 0 where A = A is positive definite. We prove that the Chebyshev method in float ing point arithmetic is numerically stable, which means that the computed sequence {x^} approximates the solution a such that lim|k is of order C ||a||.||A~ \\.\\y\\ where £ is the relative computer precision, k We also point out that in general the Chebyshev method is not well-behaved, which means that x^, k large, is not the exact solution for a slightly perturbed A or equivalently that the computed residuals r k = Ax k +g are of order C ||a|| II*" 1 || ||or||.


INTRODUCTION
Direct methods of numerical interest for the solution of linear systems Ax+g = 0 are numerically stable. This means that they produce an approximation y of the exact solution a such that |(y-a|| is of order £||a|| ||a" || ||CY|| where Q is the relative computer precision. It might seem that the numerical accuracy of iterations for solving large linear systems can be better or even not depend on the condition number of A, k (A) = ||a|| ||a ^ || . In this paper we consider the Chebyshev method which is one of the most effective iterations for the solution of large linear systems. We show that this method is stable and that the condition number of A is crucial for this iteration.
Moreover direct methods are also well-behaved which means that the computed y is the exact solution for a slightly perturbed A, i.e., (A+E)y+g = 0 where ||e|| is of order C||a|| . Unfortunately this does not hold for the Chebyshev method. Thus, from the numerical accuracy point of view direct methods seem to be better than Chebyshev.
In Section 2 we briefly recall the main properties of the Chebyshev method T [a,b] for the solution of large linear systems Ax+g = 0 where A = A is positive definite, shortly denoted by A = A > 0. In the classical case, the interval [a,b] contains all eigenvalues of A. We consider the case where b £ ||a|| and a is an arbitrary positive number. We also propose an extension of the Chebyshev method for singular matrices A = A 2: 0.
Section 3 deals with a perturbed Chebyshev method which generates a sequence { x^} such that for suitable p^ ^ and q^. We express the solution of (1.1) in terms of §^ and prove some asymptotic results.
In Section 4 we present an algorithm for the computation of p^ ^ and q^. We prove that this algorithm in floating point arithmetic computes p^_^ and with high relative precision.
Section 5 deals with the proof of numerical stability of the Chebyshev method. We prove that T[a,b] generates {x^} such that lim||x k -o/|| is of order C||A|| ||A -1 || whenever b/a is of order ||A|| |[A 1 11.
In Section 6 we discuss well-behavior of the Chebyshev method. In general, the residual vectors in the Chebyshev method r^ = Ax^+S are °* order £||A|| |JA ^ || ||o/|| which contradicts well-behavior. However, sometimes r^ can be or order CLMI IMI* Such a case yields well-behavior.

CHEBYSHEV METHOD
Let us consider the numerical solution of a large linear system and without loss of generality we can assume c^ ^ 0, for 1 £ j ^ m and X^ < < ... < X ffl , for m ^ n.
We consider a class of iterative methods which generate the sequences {x^} of the approximation of a such that where W fc is a polynomial of degree £ k. Since we only do know Ao+g • 0 than to eliminate or from (2.3) we have to assume (2.4) W k (0) -1.

Remark
Another motivation of (2.3) and (2.4) is to consider a class of iterative methods such that where and are arbitrary polynomials of degree ^ k. Assume that if x^ = a then x^ = a for any a.
Then W k (x) -U k (x)-x + 1 and Let P k (0,l) denote a class of polynomials P of degree £ k such that P(0) -1.
In the Chebyshev method T[a,b], the W k are defined as the polynomials of the smallest possible norms (2.5), i.e., (2.6) Ik || » inf ||P|| , P € P k «U) and the solution of (2.6) is given by (See, for instance, Stiefel (1958) and Rutishauser and Stiefel (1959).) Usually, the eigenvalues X, and X from (2.2) are equal to the smallest eigenvalue X . , and to the I m min largest eigenvalue, X , of A. Hence, the best convergence in T[a,b] is for a = X . and b *» X 6 & ' max* mm max However, in numerical practice X . and X are known only for a few problems. In many cases we can ' mm max easily find b ^ ^m ax (setting for instance b • ||A|| where ||« || is any matrix norm, see Young (1971), page 32). A much harder problem is to find a suitable approximation of X . . Without knowledge of X . one r mm ° mm can use the Chebyshev method T [a,b] for any values of a > 0 and b ^ ^m ax ' Then instead of (2.8) we get (2.12) l^-oll * 2q(X min ) k |^0-< ,|l 2 where v^X + V(a-X) r _ j-(2.13) q(X) -; ± ^ ^ for (a-X) + • a-X if a-X ^ 0 and zero otherwise.
Note that (i) if X € (0,a> then q(X) < 1 which means the convergence of T[a,b], however for X -> 0 + ,

Jb + Ja
(iii) if X < 0 then q(X) > 1. This implies that T[a,b] is divergence whenever X^ from (2.2) is negative.

*
One can also consider the Chebyshev method for a singular matrix A = A ^ 0. In such case by a we mean the normal solution of Ax+g 83 0, i.e., the vector of the minimal spectral norm which minimizes the spectral norm of the residual. Let g a g-j + g 2 where Ag^ • 0 and g^ is orthogonal to g 2 . Note that Aa+g 2 -0. It is straightforward to verify that {x^} defined by (2.9) in T[a,b] for a > 0 and b £ ^m ax > satisfies (2.14) x k -cy -W k (A)(x 0 -or) + \(0)g } for W k from (2.7) and where Av. • X.v., X, -0, 0 < X. < X. < ... < X , (v., v.) -6. .. Note that the normal solution or is orthogonal to v^. Let us discuss the two cases.
Thus, if c.j 88 0 (which holds for instance if x^ = 0) the Chebyshev method is convergent and the best possible speed of convergence is for a = \^, i.e., k K-lk * 2(--<) ik 0 -it-Case II. Let g 1 / 0. In that case the iterative process is divergent, although lim r^ -g^. This suggests constructing y k = x^ -w^( 0) r k' Then which once more implies the convergence of the Chebyshev method.

PERTURBED CHEBYSHEV METHOD
Recall we consider a large linear system where A • A > 0. We want to solve it by the Chebyshev method T [a,b] where it is only assumed that b £ ^m ax and a > °-The Chebyshev method generates a sequence [x^) defined by (2.9), (2.10) and (2.11).
However a sequence computed in floating point arithmetic can at best satisfy a perturbed relation (2.9), i.e.
for suitable vectors ? k » A form of § k will be discussed in Section 5. In order to analyze the Chebyshev method in fl arithmetic we start to solve (3.1) for an arbitrary {? k } an d find some asymptotical properties of the perturbed sequence t^l-Let e k • ^e the error of the kth approximant. Then from (3.1) we get Let T5K) be an arbitrary sequence and let {x k } be a perturbed sequence generated by T[a,b] defined by (3.1), (2.10) and (2.11). Then and T^, denote the Chebyshev polynomials of the first and second kind of degree k, respectively. Proof Assume now that (3.3) holds for all i £ k. Let B k = ^1 --A, W fc -W fe (A) and R^ =» R^A). Note q k that (3.4) easily follows from (2.10) and (2.11) and it is easy to verify that (3.7) W k+1 =B k W k + (1-fyW k _ 1 , (3.8) B k -P k W r From (3.2), (3.3), (3.4), (3.7) and (3.8) we get Since K K ^ -1 and R 1 -2W 1 , (3.9) follows from (3.8). To prove (3.10) we use (3.7) which holds for W k-i a n d ^t-i* B y com P arin 8 t h e coefficients at AW^^, W^^., W^^, AR^^., \ _ u ± and \_ 2 _ ± we get three equations on K K ^,

3). TO PROVE THE LIMIT O F ^ I NOTE THAT
t r 2(I+2) 2(K-I)    To prove conclusion (ii) we rewrite (3.2) as follows:

THIS COMPLETES THE PROOF OF THEOREM
Pk-iW^q k (e k-fr e k>-Since V V l -V V l and lim P k « P* -(J^f^) then lim sup ||Ae k -q k ? k || ^ (q*+p*)n -n(b+a) k which completes the proof of Corollary 3.2. •

ALGORITHM OF p^ AND q fe
In this section we deal with the computation of p k ^ and q k which appear in the Chebyshev method This suggests the following"algorithm for the computation of p k ^ and q k » We also assume that for x » rd(x), flCv/x) = Jx(l-e), | e| £ £. (See Wilkinson (1963).) To simplify the further estimations of roundoff errors we shall use the relation ~, i.e., if a(t) and b(t) are bounded functions of t, t ^ t^ > 0 then a(t) ~b(t) iff there exists K independent of t such that a(t) « b(t)(1+€(t)2" t ) where |e(t)| £ K for t £ t Q .
Next a(t) < b(t) iff a(t) £ b(t) or a(t) ~b(t). (For more details see e.g. Wozniakowski (1974).) Let us denote any computed value x in Algorithm 4.1 by x and let x • xO+1^). Thus 1^ is the relative error of x.
where 0 £ 1^ ^ 15.5 + 64K for K -Ja/b/O+Jajb) and lim 1^-3.5 • k Theorem 4.1 means that we comput«wD k and q R with high relative precision for all values of a and b.
We propose the following algorithm of the Chebyshev method (see (2.3) and Rutishauser, Stiefel and others (1959)).

NUMERICAL STABILITY OF THE CHEBYSHEV METHOD
In this section we deal with numerical stability for the Chebyshev method. Let us briefly that an iterative method for the solution of the linear equation Ax+g -0 is numerically stab produces a sequence {x^} such that (5.1) lim sup l^-crll * C^IMI II*' 1 II IHI + k where K can only depend on the size n (see Wozniakowski (1975)).
We propose the following algorithm of the Chebyshev method (see (2.3) and Rutishauser, Stiefel and others (1959)).
Xq is a given initial approximation^, for k » 0,1,... compute q k and p j by Algorithm 4.1, (5.2) r k Aa^+g; If C is small then one can prove that the constant which appears in the "0" notation in (5.10) only -1 2 depends on (||A|| || ) , K and K 9 (see (5.4). Note that if b £ L*X . then for any a 2> X . , (5.9) • ^ max mm holds; however, for increasing a the convergence of T[a,b] is getting worse, see (2.13) and (3.15).
We want to show that without additional assumption on D^, E k and© k in (5.8), the estimate (5.10) is sharp which means that the condition number of A is crucial for the accuracy of the solution of linear equation solving by iteration.
Thus, without loss of generality, a method is well-behaved if the computed y is the exact solution of the problem with a slightly perturbed matrix A.
We wish to consider the well-behavior problem for the Chebyshev method T[a,b]. This means we must verify if the computed vectors r^ = fl(Ax^+g) satisfies condition (6.4) for large k. From (5.4) we get Hr k -r*||<K lC \M IfccJI where r^ = Ae^.
Thus the Chebyshev method is well-behaved iff r^ satisfies (6.4). Let us assume for simplicity that a = X . and b = X . Note that fr. } satisfies similar recurrence formula as fx, 1, see (5.7), i.e., mm max L k * L k. J ' ' (6 -5) Vi ' r k + K-i (r kr k-i } " + A V k Numerical tests of Algorithm 5.1 confirm that (6.6) is sharp which means that in general the Chebyshev method Is not well-behaved. Note that direct methods for small dense systems such as Gaussian elimination with pivoting, the Householder method and the Gram-Schmidt reorthogonalization method are wellbehaved (see Wilkinson (1965) for two first, Kielbasinski (1974), Kielbasinski and Jankowska (1974) for the last). The lack of well-behavior for the Chebyshev method makes the termination of iteration which is based on [r^l difficult. For instance, if we want to find x^ such that ||TJJ| ^ ell^H then, in general, we can guarantee the exitence of such x^ only if e is of order £||A|| |JA ^ || ||r^ ||.
However, it can happen that (6.4) holds. Let us mention only two examples (rather theoretical).
If {g ] from (6.5) is convergent to g, ||g|| is of order C|| a ll» then applying Corollary (3.2) we get from which the well-behavior holds.
ACKNOWLEDGMENT I would like to thank H. T. Kung and J. F. Traub for their comments on this paper.