A Polynomial Preconditioner for the CMRH Algorithm

Many large and sparse linear systems can be solved efficiently by restarted GMRES and CMRH methods Sadok 1999. The CMRH m method is less expensive and requires slightly less storage than GMRES m . But like GMRES, the restarted CMRH method may not converge. In order to remedy this defect, this paper presents a polynomial preconditioner for CMRH-based algorithm. Numerical experiments are given to show that the polynomial preconditioner is quite simple and easily constructed and the preconditioned CMRH m with the polynomial preconditioner has better performance than CMRH m .


Introduction
In many practical applications, we have to solve a large and sparse, nonsymmetric linear system of equations Ax b, 1.1 where A ∈ R n×n is nonsingular, and x, b ∈ R n .The linear system 1.1 can be solved efficiently by iterative methods, especially those based on the Krylov subspace K k r 0 , A span r 0 , Ar 0 , . . ., A k−1 r 0 .

1.2
Although the GMRES method 1 and CMRH method 2 are Krylov-type algorithms, they construct different Krylov subspace bases in different ways; the GMRES method uses the Arnoldi process, while CMRH uses the Hessenberg process.As the number of iterations increases, the methods become impractical because of growth of memory and computational requirement as k increases, so they must be remedied by restarting; thus, the GMRES m and CMRH m algorithms were developed.It is concluded in 2 that the CMRH m algorithm only requires half the arithmetical operations per iteration and slightly less storage than GMRES m .However, like GMRES m , the CMRH m method may not converge in some cases.In order to remedy this defect, van Gijzen 3 presented a polynomial preconditioner to improve the GMRES m method.The main idea is finding the low-degree polynomials p m x satisfying p m A ≈ A −1 and transferring the original linear system of 1.1 into then applying the iterative method based on the following Krylov subspace: to the new linear system 1.3 , obviously, K k p m A A; r 0 ∈ K k 1 m−1 1 A; r 0 .This method becomes very effective when large computation of inner product is needed.In this paper, we apply the same technique to the CMRH m algorithm.Our method is efficient since the CMRH process itself is utilized to construct the polynomial preconditioner.Numerical experiments are given to show that the CMRH m method with a polynomial preconditioner has better performance, in the sense of fewer iterations, than the original CMRH m method.The rest of the paper is arranged as follows.The first part of Section 2 reviews the CMRH m method.The construction of the polynomial preconditioner is described in the second part of Section 2, and in the third part, the PCMRH m method is developed.We present some numerical experiments to show that the polynomial preconditioner is effective for the CMRH m method in Section 3. Section 4 is a simple conclusion.
Throughout the paper, all vectors and matrices are assumed to be real.For a vector v, v denotes the Euclidean norm v √ v T v, and v ∞ denotes the maximum norm v ∞ max i 1,...n |v i |, where v i is the ith component of a vector v, and |v i | denotes the absolute value modulus of a real complex number v i .For a matrix A, A denotes the 2-norm.Furthermore, we denote by e n k the kth canonical vector of R n , e n k 0, . . ., 1, 0, . . ., 0 T , 1.5 and I n k the first kth columns of an n-dimensional identity matrix, and we denote by H † the pseudoinverse of a matrix H.

CMRH(m) Algorithm
The principle of CMRH 2 method is to use the Hessenberg process with pivoting to construct a basis of the Krylov subspace K k A, r 0 : l 1 , l 2 , . . ., l k .Precisely, the Hessenberg process produces an n × k unit lower trapezoidal matrix L k l 1 , . . ., l k and a k 1 × k upper Hessenberg matrix H k ∈ R k 1 ×k , such that 2.1 Then, we find a correction vector z ∈ K k A, r 0 , such that min Since K k A, r 0 span{l 1 , l 2 , . . ., l k }, we let z L k y k and have min − H k y k .

2.3
Therefore, a least-squares solution could be obtained from the Krylov subspace by minimizing r 0 e k 1 1 − H k y k over y, just as in the GMRES method, except that the Hessenberg process is replaced by the Arnoldi process; this is the basic idea of the CMRH method.Since the Hessenberg process will be breakdown if l j j is zero, Sadok 2 uses a pivoting strategy such as in the Gaussian elimination method to avoid such a breakdown, and then the approximating solution can be computing by where . ., l k is an n × k unit lower trapezoidal matrix, and H † denotes the pseudoinverse of a matrix H for more details, see Algorithm 1 2 .
The notation "↔" means as "swap contents": α ↔ β ⇔ γ α; α β; β γ .As k increases, the CMRH method Hessenberg process with pivot becomes impractical because of the growth of memory and computational requirement, so it must be remedied by restarting, thus the CMRH m algorithm is developed, which is described in Algorithm 2 2 .
The main body of Algorithm 2 is the Hessenberg process with pivoting 2 , in the last, we get x m x 0 L m y m , where ; this is a little different from the GMRES method because of the pivot strategy used in the Hessenberg process.

Construction of the Polynomial Preconditioner
The main idea of the polynomial preconditioned method is to construct a polynomial s t satisfying s A ≈ A −1 and then solve the linear system of equations s A Ax s A b instead of solving the linear system 1.1 .Generally, we do not need to compute s A A really, since what we need is only that the matrix s A A approaches the unit matrix.In practice, if Algorithm 1: Hessenberg process with pivoting 2 .a polynomial q n is constructed, such that x k 1 x 0 q k A r 0 is an approximation of the solution of the linear system 1.1 , then when x 0 0, the following relationship holds: That is, q k A is a polynomial preconditioner.
In the following, we utilize the residual polynomial produced from the CMRH process to construct a polynomial preconditioner.
Let K k be an n × k matrix formed by Krylov vectors, that is, From 2.1 , we get where h k h 1k , . . ., h kk T is the kth column of H k .Since the columns of L k span the same space as the columns of K k , the following relationship holds:

8
Input: m: the dimension of a Krylov subspace and the approximation precision ε start: x 0 0, r 0 b − Ax 0 , and set p i i, i 1, . . ., n where C k is an upper triangular matrix.Thus, Hence, from 2.7 , we have

Mathematical Problems in Engineering
On the other hand, from . Since the columns of K k 1 are linear independent, we obtain Thus, from the CMRH algorithm, we can obtain the approximate solution is the residual polynomial, so we have q k z 1 − p k 1 z /z and q k A ≈ A −1 4 ; therefore, q k A becomes a preconditioner.
In the next subsection, we present the preconditioned CMRH m with the polynomial preconditioner described above.

PCMRH(m) Algorithm
The polynomial preconditioned CMRH m PCMRH m method is given in the following, which is composed of two parts, in the first part, the polynomial preconditioner is constructed, and in the second part, the CMRH m method is applied to the new linear systems Algorithm 3 .

Numerical Experiments
In this section, we give some numerical examples to compare the performance of PCMRH m with CMRH m .All codes were written in Matlab 6.5 and run on AMD Athlon tm 64 X2 Dual Core Processor 2.20 GHz equipped with 1.5 G of RAM.In these numerical experiments, all the matrices are from Sadok's paper 2 we always take the initial vector x 0 0, 0, . . ., 0 T , b 1, 1, . . ., 1 T , and m 20.We use the inequality norm r /norm b ≤ 10 −10 as stopping criteria.In the figures, solid lines -denote the CMRH m method, and dots -. denote the PCMRH m method.The x-axis denotes the number of iterations, and y-axis denotes the norm of the residual r: norm r .
Input: the dimension of Krylov subspace m; the prescribed stopping criterion ε; the degree of the polynomial preconditioner kk; the initial approximate solution x 0 ; the residual vector r 0 b − Ax 0 .1.The procedure of constructing polynomial preconditioner q kk 1 Let p i i, i 1, 2, . . ., n; Cj hj 0 end end 3 Form the approximate solution x k x 0 L k y k , where y k ∈ C k minimizes F y H k y − r 0 i0 e 1 .4 Form the polynomial preconditioner Using the CMRH m algorithm for q kk A Ax q kk A b.   is used in three experiments, whose results are listed in Figures 1 and 2 and Table 1, respectively.In Figure 1, ε 0.1, n 40, and the degree of the preconditioned polynomial: kk 20.In Figure 2, ε 0.01, n 40, and kk 20.In Table 1, ε 0.01, n 100 and kk ranges from 1 to 20.From Figures 1 and 2, we can see that the PCMRH m method converges much faster than CMRH m , when ε 0.1, n 40, the CMRH 20 method needs 107 restarts to make sure that the residual norm is 9.2925e − 011, while PCMRH 20 can get 1.1143e − 011 only in 3 restarts.When ε 0.01, n 40, the CMRH 20 method needs 840 restarts to make sure that the residual norm is 9.0075e − 011, while PCMRH 20 can get 4.4771e − 011 only by 6 restarts.The polynomial preconditioner is well done.
Let kk range from 1 to 20 and list PCMRH 20 in Table 1.We see that only when kk 1, the PCMRH 20 method does not converge, and when kk 19, PCMRH 20 has the best performance.
In this experiment, we compare the PCMRH m method with CMRH m and GMRES m .In Figure 4, we use the dot line to denote the GMRES m method.It can be seen from Figure 4 that the GMRES 20 method is stagnate since some eigenvalues of A are negative; it uses 1000 restarts or about 368 seconds to make the residual norm reach 1.8110e − 006.Meanwhile, the CMRH 20 method uses 883 restarts or about 100 seconds to get the residual norm be 8.4761e − 011; however, the PCMRH 20 method with kk 3 can reach 8.2944e − 011 only by 481 restarts or about 57 seconds.These results are also listed in Table 2 and show that the polynomial preconditioner is effective.Table 3 lists the performance of PCMRH 20 .It is seen that the PCMRH 20 method is convergent when the value of kk ranges from 1 to 20, but it can not conclude that the number of restarts and the time needed by PCMRH 20 decrease increase as kk increases decreases .In fact, when kk 3, the PCMRH 20 method performs the best.

Conclusion
To accelerate the convergence, we take advantage of the CMRH process to construct a polynomial preconditioner for CMRH m algorithm.Numerical experiments show that the preconditioned CMRH m algorithm is more efficient than the CMRH m algorithm without the polynomial preconditioner.However, the degree of the polynomial used as preconditioner should be small, otherwise it possibly reduces the preconditioner's effectiveness.
p 1 and p i 0 ;