A Polynomial Preconditioner for the CMRH Algorithm

Abstract

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).

1. Introduction

In many practical applications, we have to solve a large and sparse, nonsymmetric linear system of equations 𝐴𝑥=𝑏, (1.1) where 𝐴∈ℛ^𝑛×𝑛 is nonsingular, and 𝑥,𝑏∈ℛ^𝑛. The linear system (1.1) can be solved efficiently by iterative methods, especially those based on the Krylov subspace𝐾_𝑘𝑟₀,𝐴=span𝑟₀,𝐴𝑟₀,…,𝐴^𝑘−1𝑟₀. (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 𝑘 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 𝑝_𝑚(𝑥) satisfying 𝑝_𝑚(𝐴)≈𝐴⁻¹ and transferring the original linear system of (1.1) into𝑝_𝑚(𝐴)𝐴𝑥=𝑝_𝑚(𝐴)𝑏, (1.3) then applying the iterative method based on the following Krylov subspace:𝐾^𝑘𝑝_𝑚(𝐴)𝐴;𝑟₀=span𝑟₀,𝑝_𝑚(𝐴)𝐴𝑟₀,…,𝑝_𝑚(𝐴)𝐴^𝑘−1𝑟₀ (1.4) to the new linear system (1.3), obviously, 𝐾^𝑘(𝑝_𝑚(𝐴)𝐴;𝑟₀)∈𝐾^{(𝑘+1)(𝑚−1)+1}(𝐴;𝑟₀). 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 𝑣, ‖𝑣‖ denotes the Euclidean norm ‖𝑣‖=𝑣^𝑇𝑣, and ‖𝑣‖_∞ denotes the maximum norm ‖𝑣‖_∞=max_{𝑖=1,…𝑛}|𝑣_𝑖|, where 𝑣_𝑖 is the 𝑖th component of a vector 𝑣, and |𝑣_𝑖| denotes the absolute value(modulus) of a real(complex) number 𝑣_𝑖. For a matrix 𝐴, ‖𝐴‖ denotes the 2-norm. Furthermore, we denote by 𝑒_𝑘^(𝑛) the 𝑘th canonical vector of 𝑅^𝑛,𝑒_𝑘^(𝑛)=(0,…,1,0,…,0)^𝑇, (1.5) and 𝐼_𝑘^(𝑛) the first 𝑘th columns of an 𝑛-dimensional identity matrix, and we denote by 𝐻^† the pseudoinverse of a matrix 𝐻.

2. PCMRH(m) Algorithm

2.1. 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 𝐾_𝑘(𝐴,𝑟₀): 𝑙₁,𝑙₂,…,𝑙_𝑘. Precisely, the Hessenberg process produces an 𝑛×𝑘 unit lower trapezoidal matrix 𝐿_𝑘=(𝑙₁,…,𝑙_𝑘) and a (𝑘+1)×𝑘 upper Hessenberg matrix 𝐻_𝑘∈ℛ^{(𝑘+1)×𝑘}, such that𝐴𝐿_𝑘=𝐿_𝑘+1𝐻_𝑘. (2.1) Then, we find a correction vector 𝑧∈𝐾_𝑘(𝐴,𝑟₀), such thatmin_{𝑧∈𝐾𝑘(𝐴,𝑟0)}𝑏−𝐴𝑥₀+𝑧=min_{𝑧∈𝐾𝑘(𝐴,𝑟0)}𝑟₀−𝐴𝑧. (2.2) Since 𝐾_𝑘(𝐴,𝑟₀)=span{𝑙₁,𝑙₂,…,𝑙_𝑘}, we let 𝑧=𝐿_𝑘𝑦_𝑘 and havemin_{𝑧∈𝐾𝑘(𝐴,𝑟0)}𝑏−𝐴𝑥₀+𝑧_𝑘=min_{𝑦𝑘∈𝑅𝑘}𝑟₀−𝐴𝐿_𝑘𝑦_𝑘=min_{𝑦𝑘∈𝑅𝑘}𝑟₀−𝐿_𝑘+1𝐻_𝑘𝑦_𝑘=min_{𝑦𝑘∈𝑅𝑘}𝑟₀𝑒₁^(𝑘+1)−𝐻_𝑘𝑦_𝑘. (2.3) Therefore, a least-squares solution could be obtained from the Krylov subspace by minimizing ‖‖r₀‖𝑒₁^(𝑘+1)−𝐻_𝑘𝑦_𝑘‖ over 𝑦, 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 (𝑙_𝑗)_𝑗 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 𝑥_𝑘=𝑥₀+𝐿_𝑘𝐻^†𝑟_0𝑖0𝑒₁^𝑘+1, (2.4) where 𝑖₀ satisfies |(𝑟₀)_𝑖0|=‖𝑟₀‖_∞, 𝐿_𝑘=(𝑙₁,…,𝑙_𝑘) is an 𝑛×𝑘 unit lower trapezoidal matrix, and 𝐻^† denotes the pseudoinverse of a matrix 𝐻 for more details, see (Algorithm 1) [2].

set 𝑝=(1,…,𝑛)𝑇. Determine 𝑖0 so |(𝑣)𝑖0|=‖𝑣‖∞; 𝑙1=𝑣0(𝑣)𝑖0; 𝑝(1)↔𝑝(𝑖0) (interchange 𝑝(1)  and 𝑝(𝑖0)); iterate: for  𝑘=1,2,…,𝑛      𝑢=𝐴𝑙𝑘      for  𝑗=1,2,…,𝑘   𝑐=𝑢𝑝(𝑗); ℎ𝑗,𝑘=𝑐; 𝑢𝑝(𝑗)=0   for  𝑙=𝑗+1,…,𝑛    𝑢𝑝(𝑙)=𝑢𝑝(𝑙)−𝑐∗(𝑙𝑗)𝑝(𝑙)   end    end    if  𝑘<𝑛   determine 𝑖0  so  |(𝑢)𝑖0|=‖𝑢‖∞; 𝑝(𝑘+1)↔𝑝(𝑖0)    ℎ𝑘+1,𝑘=𝑢𝑖0; 𝑙𝑘+1=𝑢/(𝑢)𝑖0    end end

Algorithm 1 

Hessenberg process with pivoting [2].

The notation “↔” means as “swap contents”: 𝛼↔𝛽⇔𝛾=𝛼;𝛼=𝛽;𝛽=𝛾. As 𝑘 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].

Input: m: the dimension of a Krylov subspace and the approximation precision 𝜀 start: 𝑥0=0,𝑟0=𝑏−𝐴𝑥0, and set 𝑝(𝑖)=𝑖, 𝑖=1,…,𝑛 Determine 𝑖0  so |(𝑟0)𝑖0|=‖𝑟0‖∞; 𝑙1=𝑟0(𝑟0)𝑖0; 𝑝(1)↔𝑝(𝑖0) (interchange 𝑝(1) and 𝑝(𝑖0)); iterate: for  𝑘=1,2,…,𝑚      𝑢=𝐴𝑙𝑘      for  𝑗=1,2,…,𝑘      𝑐=𝑢𝑝(𝑗); ℎ𝑗,𝑘=𝑐; 𝑢𝑝(𝑗)=0      for  𝑙=𝑗+1,…,𝑛      𝑢𝑝(𝑙)=𝑢𝑝(𝑙)−𝑐∗(𝑙𝑗)𝑝(𝑙)      end      end      if  𝑘<𝑛      determine 𝑖0so |(𝑢0)𝑖0|=‖𝑢0‖∞; 𝑝(𝑘+1)↔𝑝(𝑖0)      ℎ𝑘+1,𝑘=𝑢𝑖0; 𝑙𝑘+1=𝑢/(𝑢)𝑖0      end      if (an estimate of)‖𝑏−𝐴𝑥𝑘‖ is small enough or 𝑘=𝑛 then      𝑥𝑘=𝑥0+𝐿𝑘𝑦𝑘, where 𝑦𝑘 minimize ‖𝐻𝑘𝑦−(𝑟0)𝑖0𝑒1(𝑘+1)‖      stop iteration      end end 𝑥𝑚=𝑥0+𝐿𝑚𝑦𝑚, where 𝑦𝑚∈ℛ𝑚 minimizes ‖𝐻𝑚𝑦−(𝑟0)𝑖0𝑒1(𝑚+1)‖. 𝑥0∶=𝑥𝑚, go to start

Algorithm 2 

CMRH(m) [2].

The main body of Algorithm 2 is the Hessenberg process with pivoting [2], in the last, we get 𝑥_𝑚=𝑥₀+𝐿_𝑚𝑦_𝑚, where 𝑦_𝑚∈ℛ^𝑚 minimizes ‖𝐻_𝑚𝑦−(𝑟₀)_𝑖0𝑒₁^(𝑚+1)‖; this is a little different from the GMRES method because of the pivot strategy used in the Hessenberg process.

2.2. Construction of the Polynomial Preconditioner

The main idea of the polynomial preconditioned method is to construct a polynomial 𝑠(𝑡) satisfying 𝑠(𝐴)≈𝐴⁻¹ and then solve the linear system of equations 𝑠(𝐴)𝐴𝑥=𝑠(𝐴)𝑏 instead of solving the linear system (1.1). Generally, we do not need to compute 𝑠(𝐴)𝐴 really, since what we need is only that the matrix 𝑠(𝐴)𝐴 approaches the unit matrix. In practice, if a polynomial 𝑞_𝑛 is constructed, such that 𝑥_𝑘+1=𝑥₀+𝑞_𝑘(𝐴)𝑟₀ is an approximation of the solution of the linear system (1.1), then when 𝑥₀=0, the following relationship holds: 𝐴⁻¹𝑏≈𝑥_𝑘+1=𝑞_𝑘(𝐴)𝑏. (2.5) That is, 𝑞_𝑘(𝐴) is a polynomial preconditioner.

In the following, we utilize the residual polynomial produced from the CMRH process to construct a polynomial preconditioner.

Let 𝐾_𝑘 be an 𝑛×𝑘 matrix formed by Krylov vectors, that is, 𝐾_𝑘=𝑟₀,𝐴𝑟₀,…,𝐴^𝑘−1𝑟₀. (2.6) From (2.1), we get𝑙_𝑘+1=ℎ⁻¹_𝑘+1,𝑘𝐴𝑙_𝑘−𝐿_𝑘ℎ_𝑘, (2.7) where ℎ_𝑘=(ℎ_1𝑘,…,ℎ_𝑘𝑘)^𝑇 is the 𝑘th column of 𝐻_𝑘. Since the columns of 𝐿_𝑘 span the same space as the columns of 𝐾_𝑘, the following relationship holds:𝐿_𝑘=𝐾_𝑘𝐶_𝑘, (2.8) where 𝐶_𝑘 is an upper triangular matrix. Thus, 𝐿_𝑘ℎ_𝑘=𝐾_𝑘𝐶_𝑘ℎ_𝑘=𝐾_𝑘𝐴^𝑘𝑟₀C_𝑘ℎ_𝑘0,𝐴𝑙_𝑘=𝐴𝐾_𝑘𝑐_1𝑘⋮𝑐_𝑘𝑘=𝐴𝑟₀𝐴²𝑟₀⋯𝐴^𝑘𝑟₀𝑐_1𝑘⋮𝑐_𝑘𝑘=𝐾_𝑘+10𝑐_1𝑘⋮𝑐_𝑘𝑘. (2.9) Hence, from (2.7), we have𝑙_𝑘+1=ℎ⁻¹_𝑘+1,𝑘𝐾_𝑘+10𝑐_1𝑘⋮𝑐_𝑘𝑘−𝐶_𝑘ℎ_𝑘0. (2.10) On the other hand, from 𝐿_𝑘+1=𝐾_𝑘+1𝐶_𝑘+1, we have 𝑙_𝑘+1=𝐾_𝑘+1^{𝑐1,𝑘+1⋮}_{𝑐𝑘+1,𝑘+1}. Since the columns of 𝐾_𝑘+1 are linear independent, we obtain𝑐_1,𝑘+1𝑐_2,𝑘+1⋮𝑐_{𝑘+1,𝑘+1}=ℎ⁻¹_𝑘+1,𝑘0𝑐_1𝑘⋮𝑐_𝑘𝑘−ℎ⁻¹_𝑘+1,𝑘𝐶_𝑘ℎ_𝑘0. (2.11) Thus, from the CMRH algorithm, we can obtain the approximate solution 𝑥_𝑘+1=𝑥₀+𝐿_𝑘+1𝑦=𝑥₀+𝐾_𝑘+1𝐶_𝑘+1𝑦. (2.12) Since 𝐿_𝑘+1𝑦=𝑥₀+𝐾_𝑘+1𝐶_𝑘+1𝑦, we know that the entries of 𝐶_𝑘+1𝑦 are the coefficients of 𝑞_𝑘(𝑧), that is, if 𝐶_𝑘+1𝑦=(𝛼₀,𝛼₁,…,𝛼_𝑘)^𝑇, then 𝑞_𝑘(𝑧)=𝛼₀+𝛼₁𝑧+⋯+𝛼_𝑘𝑧^𝑘. Note that 𝑝_𝑘+1(𝑧)=1−𝑧𝑞_𝑘(𝑧) is the residual polynomial, so we have 𝑞_𝑘(𝑧)=(1−𝑝_𝑘+1(𝑧))/𝑧 and 𝑞_𝑘(𝐴)≈𝐴⁻¹ [4]; therefore, 𝑞_𝑘(𝐴) becomes a preconditioner.

In the next subsection, we present the preconditioned CMRH(m) with the polynomial preconditioner described above.

2.3. 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).

Input: the dimension of Krylov subspace 𝑚; the prescribed stopping criterion 𝜀;   the degree of the polynomial preconditioner 𝑘𝑘; the initial approximate solution 𝑥0;       the residual vector 𝑟0=𝑏−𝐴𝑥0. 1. The procedure of constructing polynomial preconditioner 𝑞𝑘𝑘  (1)  Let 𝑝(𝑖)=𝑖,𝑖=1,2,…,𝑛;      determine 𝑖0 so |(𝑟0)|𝑖0=‖𝑟0‖∞; 𝑙1=𝑟0(𝑟0)𝑖0; 𝑝(1)↔𝑝(𝑖0),      1𝐶(1,1)=(𝑟0)𝑖0   (2)  for 𝑘=1,2,…,𝑘𝑘 do     𝑢=𝐴𝑙𝑘         for𝑗=1,2,…,𝑘 do    𝑐=𝑢𝑝(𝑗),ℎ𝑗,𝑘=𝑐,𝑢𝑝(𝑗)=0    for𝑙=𝑗+1,…,𝑛 do       𝑢𝑝(𝑙)=𝑢𝑝(𝑙)−𝑐×(𝑙𝑗)𝑝(𝑙)    end      if 𝑘<𝑛     compute 𝑖0 so |(𝑢0)|𝑖0=‖𝑢0‖∞; 𝑝(𝑘+1)↔𝑝(𝑖0)     ℎ𝑘+1,𝑘=𝑢𝑖0; 𝑙𝑘+1=𝑢(𝑢)𝑖0     let ℎ𝑗=ℎ1,𝑗⋮ℎ𝑗,𝑗,𝐶𝑗=𝑐11⋯𝑐1𝑗0⋱⋮00𝑐𝑗𝑗     compute     𝑐1,𝑗+1𝑐2,𝑗+1⋮𝑐𝑗+1,𝑗+1=ℎ−1𝑗+1,𝑗0𝑐1𝑗⋮𝑐𝑗𝑗−ℎ−1𝑗+1,𝑗𝐶𝑗ℎ𝑗0   end  end   (3) Form the approximate solution 𝑥𝑘=𝑥0+𝐿𝑘𝑦𝑘, where 𝑦𝑘∈𝐶𝑘 minimizes             𝐹(𝑦)=‖𝐻𝑘𝑦−(𝑟0)𝑖0𝑒1‖. (4) Form the polynomial preconditioner 𝑞𝑘𝑘(𝑧)             𝐶𝑘𝑦𝑘=(𝛼0⋯𝛼𝑘−1)𝑇;             𝑞𝑘𝑘(𝑥)=𝛼0+𝛼1𝑥+⋯+𝛼𝑘−1𝑥𝑘−1.         2. Using the CMRH(m) algorithm for 𝑞𝑘𝑘(𝐴)𝐴𝑥=𝑞𝑘𝑘(𝐴)𝑏.

Algorithm 3 

PCMRH(m).

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 𝑥₀=(0,0,…,0)^𝑇, 𝑏=(1,1,…,1)^𝑇, and 𝑚=20. We use the inequality norm(𝑟)/norm(𝑏)≤10⁻¹⁰ as stopping criteria. In the figures, solid lines(-) denote the CMRH(m) method, and dots(-.) denote the PCMRH(m) method. The 𝑥-axis denotes the number of iterations, and 𝑦-axis denotes the norm of the residual 𝑟: norm(𝑟).

Example 3.1. The first example is taken from Brown [5]. The 𝑛×𝑛 matrix𝐴=𝜀1−1𝜀1⋱⋱⋱−1𝜀1−1𝜀 (3.1) is used in three experiments, whose results are listed in Figures 1 and 2 and Table 1, respectively. In Figure 1, 𝜀=0.1, 𝑛=40, and the degree of the preconditioned polynomial: 𝑘𝑘=20. In Figure 2, 𝜀=0.01, 𝑛=40, and 𝑘𝑘=20. In Table 1, 𝜀=0.01, 𝑛=100 and 𝑘𝑘 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, 𝑛=40, the CMRH(20) method needs 107 restarts to make sure that the residual norm is 9.2925𝑒−011, while PCMRH(20) can get 1.1143𝑒−011 only in 3 restarts. When 𝜀=0.01, 𝑛=40, the CMRH(20) method needs 840 restarts to make sure that the residual norm is 9.0075𝑒−011, while PCMRH(20) can get 4.4771𝑒−011 only by 6 restarts. The polynomial preconditioner is well done.
Let 𝑘𝑘 range from 1 to 20 and list PCMRH(20) in Table 1. We see that only when 𝑘𝑘=1, the PCMRH(20) method does not converge, and when 𝑘𝑘=19, PCMRH(20) has the best performance.

Table 1 

Figure 1 

𝜀=0.1, 𝑛=40.

Figure 2 

𝜀=0.01, 𝑛=40.

Example 3.2. The matrix𝑎𝐴=111…11₁𝑎11…11₁𝑎₂𝑎1…11⋮⋮⋱⋱⋱⋮₁𝑎₂𝑎₃…𝑎_𝑛−11with𝑎_𝑖=1+𝑖𝜀 (3.2) is taken from [6]. In this experiment, we let 𝜀=10⁻², 𝑛=100, and 𝑘𝑘=2.
From Figure 3, we can also see that the PCMRH(20) method converges much faster than CMRH(20). The CMRH(20) method needs 317 restarts to make sure the residual norm is 9.9221𝑒−011, while PCMRH(20) reaches 9.0007𝑒−011 only in 37 restarts. So, we can conclude that the polynomial preconditioner makes CMRH(m) algorithm more efficient.

Figure 3 

𝜀=10−2, 𝑛=100.

Example 3.3. Let 𝐴=𝑆𝐷𝑆⁻¹∈ℛ^1000×1000, where 𝑆 is a bidiagonal matrix of the form1𝑆=10.910.9⋱⋱0⋱0.9∈ℛ^1000×1000, (3.3) and 𝐷=diag(−10,−9,…,−1,1,2,…,990)∈ℛ^1000×1000.
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 𝐴 are negative; it uses 1000 restarts or about 368 seconds to make the residual norm reach 1.8110𝑒−006. Meanwhile, the CMRH(20) method uses 883 restarts or about 100 seconds to get the residual norm be 8.4761𝑒−011; however, the PCMRH(20) method with 𝑘𝑘=3 can reach 8.2944𝑒−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 2 

Figure 4 

𝐴=𝑆𝐷𝑆−1.

Table 3 lists the performance of PCMRH(20). It is seen that the PCMRH(20) method is convergent when the value of 𝑘𝑘 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 𝑘𝑘 increases(decreases). In fact, when 𝑘𝑘=3, the PCMRH(20) method performs the best.

Table 3 

4. 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.

Acknowledgment

This work is supported by the National Natural Science Foundation of China no. 10961010.