Fourier Analysis of GMRES (m) preconditioned by multigrid
Fourier Analysis of GMRES(m) preconditioned by multigrid R. Wienands
Abstract:
This paper deals with convergence estimations of a preconditioned GMRES(m) method [6], where multigrid ([1], [8]) is used as the preconditioner. Fourier analysis is a well-known useful tool in the multigrid community for the prediction of two-grid convergence rates ([1], [8]). This analysis is generalized to the situation where multigrid is a preconditioner.
1 Introduction Nowadays, it has become popular to study the convergence of multilevel methods and use them in combination with a Krylov subspace acceleration method. Recently, also standard (multiplicative) multigrid methods are used as a preconditioner. This application of standard multigrid is bene cial in situations where standard multigrid alone does not converge fully satisfactorily, because certain error frequencies are not reduced well enough, which mainly occurs when complicated PDEs are solved. In general, it is dicult to construct robust multigrid solvers for large classes of problems. Among many other papers, this approach has been presented in [5] for linear singularly perturbed problems. Here, we analyze one combined solution method, the multigrid preconditioned GMRES(m) solution method, theoretically. We derive quantitative convergence results in order to understand the convergence and the bene ts of the combination. The basis for the theoretical convergence estimations is Fourier analysis, which is a well-known tool for obtaining two-grid convergence factors [8]. In this paper, we restrict ourselves to standard problems coming from scalar linear PDEs discretized with nite di erences on uniform Cartesian grids and we keep the multigrid method simple. For example, we restrict ourselves to point smoothers and standard grid coarsening and analyze the e ect of Krylov subspace acceleration for anisotropic di usion problems and for problems with mixed derivatives, although it is known that other (more expensive) smoothers ore other coarsening strategies are appropriate remedies for such equations. However, in this way we try to gain insight in the behavior of the combination of multigrid and GMRES(m). The typical use of the Fourier analysis in a multigrid context is that the spectral radius is obtained theoretically, which is the basis for asymptotic multigrid convergence estimations [8]. Two-grid Fourier analysis transforms a multigrid iteration matrix into a block diagonal matrix from which more information about the spectrum is easily obtained. The analysis makes use of unitary basis transformations to simplify the representation of the multigrid iteration matrix. It is therefore appropriate for the estimation of the multigrid preconditioned GMRES(m) convergence since unitary transformations do not a ect the convergence behavior of GMRES [6]. A preconditioned Krylov subspace acceleration method like GMRES(m) implicitly builds up a minimal residual polynomial. The determination of the polynomial coecients is easily possible with the help of the Fourier analysis and can be done explicitly, since a block diagonal matrix results. Based on the GMRES(m) polynomial and the Fourier analysis, sharp theoretical convergence estimates can be obtained. In section 2 the multigrid preconditioner is presented. The two-grid Fourier analysis results are used in section 3, where di erent ways to analyze restarted GMRES(m) preconditioned by multigrid are discussed. There are actually two types of Fourier analysis available in the multigrid community. The rst type is the analysis called \rigorous analysis" here, since it explicitly takes boundary conditions into account and discrete eigenfunctions are related to the discrete mesh size. The second is known as local mode analysis [1]. Exponential functions are the basis for the analysis, a GMD, Institute for Algorithms and Scienti c Computing, D-53754 Sankt Augustin, Germany
(email: wienands @ gmd.de)
1
discretization at domain boundaries cannot be taken into account. This analysis can be used for a wider range of PDEs and multigrid components. The generalization of Fourier analysis to multigrid as a preconditioner is applied for both types of analysis. By several numerical tests the theoretical convergence estimations are compared to the actual numerical convergence. The multigrid purist may not like the approach presented here with multigrid viewed as a preconditioner. It is, of course, also possible to try to construct an optimal multigrid method for each separate problem. For the purist, the analysis might give an indication how far his present method is from an optimal multigrid solution method. In other words, if an additional acceleration of the convergence is found by the Fourier analysis, it must also be possible to tune one of the multigrid components, so that multigrid as a stand alone solver is improved.
2 Rigorous Fourier Analysis of Multigrid The 3D Poisson equation, a simple test problem, is chosen to explain the possibilities of convergence estimations by Fourier analysis for multigrid as a solver. We apply the standard 7-point discretization on a cube with Dirichlet boundary conditions and use a uniform mesh with meshsize h = 1=n:
Ah uh(x) = ?h uh(x) = bh (x) on h = (0; 1)3 \ Gh ; (1) uh(x) = 0 on ?h = ([0; 1]3 n h ) \ Gh with Gh := fx = hj ; j 2 Z3g. We combine operator and matrix notation. Multigrid is
commonly explained with operators (see, for example, [8], [1]), whereas Krylov subspace methods are often described with matrices. A well-known ecient 3D multigrid method [9] for problem (1) consists of a red-black GaussSeidel smoothing method (GS-RB), Full Weighting of residuals to obtain the right-hand side on the coarse grids, tri-linear interpolation of corrections from coarse to ne grids and the coarse grid discretization of the PDE on a grid with meshsize 2h in each direction. We call this method the red-black multigrid solver (MG-RB) here. For problem (1) and MG-RB, it is possible to apply rigorous two-grid Fourier analysis, as it was presented for 2D problems in [8] and used for 3D (Poisson-type) problems in [9]. The existence of the orthonormal basis of discrete eigenfunctions
'k;`;m h (x; y; z ) = sin kx sin `y sin mz with k; `; m = 1; : : :; n ? 1 ((x; y; z ) 2 h )
(2)
for the operator (1) is crucial for the 3D discrete Fourier analysis. They generate the space of all grid functions, F ( h ). The error vj ?1 is transformed by a two-grid cycle as follows: (3) vj = Mh2h vj?1 with Mh2h = Sh2 Kh2h Sh1 = Sh2 (Ih ? I2hh (2h)?1Ih2h h )Sh1 ; where Sh is the smoothing operator, 1 and 2 indicate the number o f pre- and post-smoothing iterations, 2h is the discretization of the operator on a 2h coarse grid, I2hh and Ih2h are transfer
operators from coarse to ne grids and reversed. In the case of standard multigrid coarsening (H = 2h) in 3D, it is convenient to divide F ( h ) into a direct sum of (at most) eight-dimensional subspaces, the 2h-harmonics Ehk;`;m, which leaves Mh2h invariant (see [9], [8]). Hence, Mh2h is f consisting of (at most) 8 8 blocks ([8], [9]): orthogonally equivalent to a block matrix M h
i
^ M 2h j c(k; `; m) = f := M c(k; `; m) =^ Mh2h with M M h Ehk;`;m : k;`;m=1;::: ; n 2
(4)
The tilde e indicates the Fourier matrix representation of the related operator. By M we denote the matrix representation of the multigrid operator Mh2h and by A the matrix representation of f = UMU ?1 with a unitary matrix U , which h with respect to the Euclidian base. We have M 2
rows are given by the eigenvectors of A (for the corresponding eigenfunctions see (2)). From the c(k; `; m), we obtain the two-grid convergence factor by block matrices M f) = c(k; `; m)) : F := (M max (M 1k;`;m n 2
(5)
It is the asymptotic two-grid convergence which is a theoretical estimation of the multigrid convergence. It is furthermore possible to determine the whole spectrum and thus the eigenvalue f by calculating the eigenvalues of the block matrices M c(k; `; m). We use this distribution of M distribution in section 3 for the analysis of the convergence of the multigrid preconditioned GMRES(m) method.
3 Fourier analysis of multigrid preconditioned GMRES(m) In this section we discuss two approaches for the analysis of multigrid preconditioned GMRES(m). For both approaches it is crucial that the discrete Fourier analysis makes use of the unitary basis transformations which simplify the representation of the multigrid iteration matrix to the block matrix (4). We consider the linear system related to (1), Au = b. A two-grid (or more general a multigrid) cycle can be described by the following matrix splitting: Cuj + (A ? C )uj?1 = b ; (6) where uj represents a new and uj ?1 a previous approximation. The multigrid method (6) considered here is described in the previous section. It is used as a right preconditioner for GMRES(m). The parameter m represents the number of vectors stored after which the GMRES(m) algorithm will restart. GMRES(m) searches for a solution uj in the following subspace: C (uj ? uj?m ) 2 span[rj?m ; (AC ?1)rj?m ; :::; (AC ?1)m?1 rj?m ] =: K m(AC ?1 ; rj?m ) ; (7) where K m (AC ?1; rj ?m ) is the Krylov subspace. It selects the new solution vector uj by minimizing the corresponding residual rj in the L2-norm minm ?1 jjb ? Auj jj2 : (8) C (uj ?uj?m )2K (AC ;rj?m ) the ane subspace uj ?m + C ?1 K m (AC ?1 ; rj ?m) can be
Any element w in represented by w = uj?m + C ?1 ( 1rj?m + 2 AC ?1 rj?m + + m (AC ?1 )m?1 rj?m ): (9) By substituting (9) into the residual equation, we obtain rw = b ? Aw = Pm (AC ?1)rj?m ; (10) P
where Pm is a m-th order polynomial de ned by Pm () = 1 ? mk=1 k k . Comparing (8) and (10), we see that the L2-norm of rj satis es the following minimization property: jjrj jj2 = minfjjPm(AC ?1)rj?mjj2 j Pm 2 Pmg; (11) where Pm denotes the set of all polynomials of degree at most m with Pm (0) = 1. Since we are interested in the convergence behavior of multigrid preconditioned GMRES with a restart after m iterations we consider the residuals rm , r2m, ..., rim . This leads to the following de nition of the reduction factors i for a complete iteration, consisting of m multigrid preconditioned GMRES(m) steps: 1=mi1=i h : i := jjjjrrimjjjj2 0 2
3
(12)
In section 3.1 we review estimations for i which are based on an analysis of the spectrum of AC ?1 (see [6], [7]). In general, these estimations exhibit a qualitative character and concrete values are rarely found in the literature because the spectrum is usually not available. However, with multigrid as the preconditioner these estimations can be calculated explicitly since the spectrum is easily obtained. In section 3.2 we will introduce a sharper estimation, which is based on an analysis of the GMRES(m) polynomial.
3.1 Analysis with the spectrum of the iteration matrix
The common way to analyze the convergence of GMRES is to exploit information about the spectrum of the iteration matrix AC ?1 . If we compare the error transformation (3) and the f (4) is related to the two-grid residual transformation by one two-grid cycle, the Fourier matrix M ?1 iteration matrix I ? AC via the identity f = UA?1 (I ? AC ?1 )AU ?1 () M f = I ? UC ?1 AU ?1 ; M with the unitary matrix U from (4). With (13) we have
(13)
f(UAU ?1 )?1 = I ? A eM fA e?1 : UAC ?1 U ?1 = I ? UAU ?1 M
(14)
Ae and Ae?1 are diagonal matrices since the rows of U are eigenvectors of A. This means that UAC ?1 U ?1 is a block-diagonal matrix which consists of (at most) 8 8 blocks in the 3D case and of (at most) 4 4 blocks in the 2D case. The eigenvalues of these block matrices are easily calculated. First, we discuss the considerations on the convergence of GMRES(m) from [6] and [7], which are based on the spectrum of the iteration matrix. From (11) it follows that
jjrimjj2 jjPm(AC ?1)jj2 jjr(i?1)mjj2 for all Pm 2 Pm : (15) Suppose that all eigenvalues of AC ?1 are located in an ellipse E (c; d; a) which excludes the origin. c denotes the center, d the focal distance and a the major semi-axis. Then, it is known [7] that the absolute value of the polynomial tm (z) := Tm ( dc ? 1d z)=Tm( dc ) = Tm(zb)=Tm( dc ) with z; zb := ( dc ? d1 z) 2 C
(16)
is small on the spectrum of AC ?1 . Here, Tm represents the Chebychev polynomial of degree m of the rst kind which is given by the three-term recurrence relation
T0(Ab) = I ; T1 (Ab) = Ab ; Tm+1 (Ab) = 2Tm(Ab) ? Tm?1 (Ab) for m 1 : If AC ?1 is diagonalizable by AC ?1 = XDX ?1 then (15) yields
jjrimjj2 jjtm(AC ?1)jj2 jjr(i?1)mjj2 ii h 2(X ) Tm( da )=Tm( dc ) jjr0jj2
h
ii
jjtm(AC ?1)jj2 jjr0jj2
(17) (18) (19)
where 2 (X ) denotes the spectral condition number of the transformation matrix X [7]. Inequality (19) uses the fact that the maximum modulus of a complex analytical function is reached on the boundary of the domain which is due to Liouville's theorem [7]. Using the above relations (18) and (19) we obtain the following well known estimations ([7]) for i
i NmE (2 (X ))1=m TmE 1=m : with NmE := (jjtm(AC ?1 )jj2)1=m and TmE := Tm ( da )=Tm ( dc ) 4
(20) (21)
TmE can be further approximated by TmE
=
p
+ ( ad )2 ? 1 m + da + ( ad )2 ? 1 ?m 1=m a + a2 ? d2 p2 2: t p c 2 ? 1m + ? c + p( c )2 ? 1?m c + c ?d + ( ) d d d d
?a d ?c
p
?
p
(22)
The approximation from (22) is sharp if a d and c d or m 1 holds. The rst inequality is ful lled if the ellipse E (c; d; a) tends to a circle which means that the focal distance d becomes very small. In all cases considered, the multigrid preconditioner ensures that the spectrum is clustered around 1, which means that the second inequality is satis ed. If (22) is sharp due to a circular shape of or due to a large m, TmE is independent of the restart parameter m. In this case, it is dicult to accelerate the multigrid preconditioner. A large Krylov subspace (7) does not lead to additional acceleration. The main drawback of the estimations (20) for the reduction factors i is their i-independence. Thus, they can not expected to be sharp approximations for every i. Furthermore the constant 2 (X ) can be very large if the iteration matrix AC ?1 is far from normal but it is found in [4] that the condition number 2 (X ) is not always relevant for the convergence of GMRES. In [7] it is stated that (20) is an asymptotic result and that the actual residual norm should behave like TmE without 2 (X ). Therefore, explicit values also for TmE are quoted in section 3.3 and compared with asymptotic numerical convergence results. However, TmE is a heuristic estimation for i 1 and not an upper bound like NmE .
3.2 Analysis with the GMRES-polynomial
To obtain a quantitative result for the convergence of GMRES(m), we present a second approach to estimate i, which explicitly depends on the iteration index i in contrast to the previous section. An initial residual r0 is prescribed and the GMRES(m)-polynomials Pmi are explicitly determined for every i. The choice of r0 does not in uence the average reduction factor for i >> 1. For any unitary matrix U , in particular for U from (4), we have jjrimjj2 = jjUrimjj2. In order to nd the coecients ik (k = 1; :::; m) of the GMRES(m)-polynomials Pmi , the following function f is minimized:
f ( i1 ; : : :; im) := UPmi (AC ?1 )r(i?1)m; UPmi (AC ?1 )r(i?1)m =
m X
(23)
l;k=1
+2
il ik (UAC ?1 U ?1 )lUr(i?1)m; (UAC ?1U ?1 )k Ur(i?1)m
m X l=1
il (UAC ?1 U ?1 )lUr(i?1)m ; Ur(i?1)m :
(24)
The ik are obtained by solving the following linear system of equations: m @f =2 X ?1U ?1 )k Ur i (UAC ?1 U ?1 )l Ur ; ( UAC (i?1)m (i?1)m @ il k=1 k + 2 (UAC ?1 U ?1 )l Ur(i?1)m; Ur(i?1)m = 0 for l = 1; : : :; m :
(25)
fA e?1 of the GMRES(m) polynomials In (14) it is seen that the argument UAC ?1 U ?1 = I ? AeM Pmi has a simple block structure. Therefore, the solution of the linear system (25) can be easily computed for every complete iteration i due to the sparse structure of (UAC ?1U ?1 )l (l = 1; :::; m), if the previous transformed residual Ur(i?1)m is given. We prescribe right-hand side b = 0 and a randomly chosen initial approximation Uu0, which yield the transformed initial residual
5
P
n?1 (k; `; m)'k;`;m e 0 . More precisely, we select an initial solution u0 = Ur0 = ?AUu k;`;m=1 where the amplitudes (k; `; m) of the eigenvectors are randomly chosen. Hence we get Uu0 = ( (1; 1; 1); (n ? 1; n ? 1; n ? 1); : : :; ( n2 ; n2 ; n2 ) )T , see (2). From Ur0 it is possible to evaluate the coecients 1k by solving (25) for i = 1. This gives Urm = Pm1 (UAC ?1U ?1 )Ur0. It follows that Urim can be easily calculated for every i. This leads to the de nition of an average reduction factor Ui , which can be obtained by the Fourier analysis: h jjUr jj 1=m i1=i U : (26) i := jjUrimjj 2 0 2 In most of the tests, presented in the next subsection, Ui becomes constant for i 20. Thus, U20
is expected to match very well with numerical reference values and in particular to be sharper than the estimation from section 3.1.
3.3 Fourier results of multigrid preconditioned GMRES(m)
Here, we investigate the quality of the theoretical considerations by applying the multigrid preconditioned GMRES(m) method to the 3D Poisson problem. We compare the average accelerated multigrid convergence rate acc h from numerical calculations with the theoretical prediction U E i20, the heuristic estimation Tm and the upper bound NmE . Insight in the behavior of the acceleration method is gained by considering di erent combinations of the relaxation parameter ! and the GMRES(m) restart parameter m. In [10] it is presented that an overrelaxation parameter (! > 1) improves the smoothing properties of GS-RB (leading to ! -GS-RB) and therefore improves the convergence of MG-RB for d-dimensional Poisson-type equations. We investigate if it is possible to further accelerate the multigrid methods with these relaxation parameters by GMRES(m). The results are presented in Table 1. We observe that Ui20 and TmE provide accurate predictions of the multigrid preconditioned GMRES(m) convergence acc h for di erent E overrelaxation parameters ! . The upper bound Nm becomes sharper with increasing m as TmE and NmE coincide for large m, because the condition number 2 (X ) cancels out (see (20)). Table 1: Multigrid convergence mg convergence acc h , predicted converh , acceleratedE multigrid U E gence F , i20 and convergence estimations Tm , Nm from Chebychev polynomials for the 3D Poisson equation; W(1,1), h = 1=32. Relaxation MG-RB acc. MG-RB, m = 2 acc. MG-RB, m = 5 acc U E E acc parameter mg T N Ui20 TmE NmE F m m i20 h h h ! = 1:0 0.192 0.194 0.085 0.086 0.088 0.311 0.070 0.074 0.085 0.149 ! = 1:1 0.089 0.097 0.045 0.056 0.059 0.272 0.042 0.049 0.059 0.102 ! = 1:15 [10] 0.070 0.077 0.050 0.057 0.060 0.234 0.045 0.048 0.060 0.102 In Table 1, we observe a satisfactory convergence improvement for ! = 1. If overrelaxation is selected (! = 1:1; 1:15) the bene ts of the GMRES(m) acceleration are small and the convergence cannot be improved by using a larger Krylov space (7). This behavior is expected from the spectral picture in Figure 1. The ellipses become orbital (d a) with the overrelaxation so that (22) holds by which a m-independent acceleration is predicted.
4 Local mode analysis of multigrid preconditioned GMRES(m) The previous Fourier analysis is restricted to model operators Ah , for which the discrete functions (2) form a basis of eigenfunctions. Furthermore, it is not possible to use the analysis for multigrid methods based on Gauss-Seidel relaxations with a lexicographical ordering of the grid points (GSLEX). The eigenfunctions (2) are no longer invariant under the GS-LEX operators. Instead they 6
ellipse E(c=0.903,d=0.078,a=0.097) spectrum
0.04
0.04
0.02
0.02
0
0
-0.02
-0.02
-0.04
-0.04
-0.06
-0.06
0.8
0.85
0.9 REAL PART
0.95
ellipse E(c=0.959,d=0.005,a=0.056) spectrum
0.06
IMAGINARY PART
IMAGINARY PART
0.06
1
0.9
0.92
0.94
0.96 REAL PART
0.98
1
1.02
Figure 2-a: ! = 1:00 Figure 2-b: ! = 1:10 Figure 1: Eigenvalue spectra of the Red-Black Two-Grid 3D Poisson solver from rigorous Fourier analysis.(Note the di erent scaling on the real- and imaginary-axis in Fig. 2-a.) are intermixed, which leads to a full matrix Se. For a general linear operator Ah with constant or frozen coecients it is, however, possible to use the local mode Fourier analysis [1] for analyzing multigrid methods based on GS-LEX as the smoother. The local mode (or in nite grid) analysis gives satisfactory estimations of the actual multigrid convergence. Its theoretical justi cation is, however, not as straightforward as for the previous analysis [2]. Instead of considering the nite domain with discrete eigenfunctions (2) of Ah an in nite domain h := fx = (jx h; jy h; jz h) : jx ; jy; jz 2 Zg is considered with continuous eigenfunctions h (; x) = eix=h = eij = ei(jx x +jy y +jz z ) ; (27) where the Fourier frequencies = (x ; y ; z ) vary continuously in = (?; ]3. The Fourier space "h = spanfeij : 2 g can be divided into eight-dimensional subspaces, the 2hharmonics "h , which leaves Mh2h invariant [8]. In order to get a well-de ned two-grid operator we replace "h by a slightly shrunk subspace "~h = spanfeij : 2 n g such that (Ah )?1 exists and also (A2h )?1 can be reasonably de ned on the coarser grid, like in [8]. The spectral radius f) of the corresponding block diagonal Fourier two-grid matrix is determined similarly as in (M the previous Fourier analysis (5): ^ M 2h j : f) = sup (M c()) with ~ l = (?=2; =2]3n and M c() = F := (M (28) h "h ~l 2
The whole spectrum is again obtained similarly as in section 2. We, however, always observe a spectrum from a discrete mesh, whereas we have a continuously de ned set of eigenfunctions (27) in order to get h-independent upper bounds for the convergence rates. With the local mode analysis the analytical values Ui20 and TmE , NmE for the accelerated multigrid method can be obtained by straightforward modi cations of the de nitions (26) and (21) taken from the rigorous case.
4.1 Results for Poisson-type equations
The rst example for the local mode Fourier analysis of multigrid preconditioned GMRES(m) deals with the anisotropic di usion equation. We analyze the di erence in convergence between 7
performing two smoothing iterations (1 = 2 = 1) of the forward lexicographical point GaussSeidel (GS-fLEX) method and one iteration of forward lexicographical point Gauss-Seidel and one iteration of backward lexicographical point Gauss-Seidel (GS-bLEX). (The other multigrid components are identical to the previously introduced multigrid method MG-RB.) We call these multigrid methods MG-FF and MG-FB, respectively, in the following. Local mode analysis yields sharp quantitative predictions for the multigrid convergence of MG-FF and MG-FB for the 2D anisotropic Poisson equation, as can be seen in Table 2. In the anisotropic case it is bene cial to introduce an overrelaxation parameter also for these lexicographical smoothers leading to ! -MG-FF and ! -MG-FB. (The improvement is less impressive than for ! -MG-RB in [10], even if the optimal overrelaxation parameter (!opt = 1:40 for " = 0:1, !opt = 1:75 or " = 0:01) is selected.) Table 2: Multigrid convergence mg h and predicted convergence F for the 2D anisotropic Poisson equation; h = 1=128. solution " = 0:1, ! = 1:00 " = 0:1, ! = 1:40 " = 0:01, ! = 1:00 " = 0:01, ! = 1:75 method mg F mg F mg F mg F h h h h MG-FF 0.693 0.696 0.410 0.433 0.957 0.961 0.755 0.758 MG-FB 0.693 0.697 0.437 0.440 0.957 0.962 0.750 0.759 However, in contrast to !opt -MG-FF (and !opt-MG-RB), it is possible to further accelerate !opt-MG-FB signi cantly with GMRES(m). This can be observed by comparing Table 2 and Table 3 and can be explained as follows. MG-FB results in a symmetric iteration matrix due to U Table 3: Accelerated multigrid convergence acc h , predicted convergence i20 and estimations TmE , NmE from Chebychev polynomials for the 2D anisotropic Poisson equation; h = 1=128. solution " = 0:1, ! = 1:00 " = 0:1, ! = 1:40 acc U E E acc method h i20 Tm Nm h Ui20 TmE NmE acc. MG-FF, m = 2 0.475 0.515 0.529 0.911 0.410 0.420 0.430 0.524 acc. MG-FF, m = 5 0.450 0.469 0.519 0.907 0.410 0.420 0.430 0.441 acc. MG-FB, m = 2 0.370 0.402 0.408 0.718 0.185 0.185 0.203 0.445 acc. MG-FB, m = 5 0.305 0.331 0.331 0.487 0.158 0.160 0.165 0.202 solution " = 0:01, ! = 1:00 " = 0:01, ! = 1:75 U E E acc U method acc T N TmE NmE m m i20 i20 h h acc. MG-FF, m = 2 0.850 0.895 0.907 0.952 0.753 0.750 0.757 0.923 acc. MG-FF, m = 5 0.800 0.837 0.896 0.897 0.753 0.750 0.757 0.813 acc. MG-FB, m = 2 0.784 0.850 0.865 1.403 0.397 0.380 0.479 0.717 acc. MG-FB, m = 5 0.723 0.755 0.768 1.102 0.362 0.370 0.392 0.457
the construction of the coarse grid correction (FW as the restriction and bilinear interpolation as the prolongation operator) and the two relaxations which are adjoined to each other. This means that the spectrum of the iteration matrix is real-valued, and uniformly distributed on the interval [min ; max]. This feature is maintained even if !opt is introduced, whereas for !opt -MGFF an acceleration by GMRES(m) is not possible, because the spectrum becomes orbital. It is found in [3] that a symmetrization is not bene cial if Bi-CGstab is used as the Krylov subspace acceleration method. The analytical prediction Ui20 shows a sharp quantitative character again and also the estimation TmE gives a reliable indication for the average accelerated convergence rate presented in Table 3. NmE is an interesting quantity for large m. Remark: For symmetric matrices it is possible to use the Conjugate Gradient (CG) iteration as acceleration method. The asymptotic convergence rate cg for the preconditioned CG iteration 8
is well-known. If we consider the de nition TmE (20) for a real-valued spectrum [min ; max] then the ellipse E (c; d; a) deteriorates to a line. In particular, we have a = d = (max ? min )=2 and c = min + a, which leads to TmE = (Tm ( aa )=Tm ( ac ))1=m = cg (see (21), (22)). Thus it can be expected heuristically that GMRES(m) has the same asymptotic convergence as CG. Remark: Also for the 3D Poisson equation the multigrid convergence of the two methods MG-FF and MG-FB is similar. The theoretically obtained convergence factor F is 0.321 for both methods, although the spectra are completely di erent (as in 2D). Very satisfactory accelerated convergence rates are obtained by the Krylov subspace acceleration of MG-FB. Due to the symmetry we can predict the convergence rate of the accelerated MG-FB accurately from cg .
4.2 Results for a problem with mixed derivative
Finally, we discuss an equation where reasonable smoothing properties can be achieved with point smoothers, but the coarse grid correction turns out to be problematic in combination with standard coarsening [8]. Again an acceleration with the help of GMRES(m) can be applied as a remedy to overcome this diculty. We investigate the 2D di erential equation ?u ? uxy = b with a mixed derivative discretized via the O(h2 ) 9-point operator 2 3 4 ?1 ? 4 1 Ah uh (x) = h2 4 ?1 4 ?1 5 uh (x) = bh (x) on h : (29) ? 4 ?1 4 The di erential equation is no longer elliptic if j j = 2 and the eciency of all previously considered multigrid methods (MG-RB, MG-FF, MG-FB) deteriorates for j j ! 2 (see for example Table 4 and [8]). This behavior can be explained if we perform a simpli ed two-grid analysis where only the very low Fourier frequencies along the characteristic direction of the di erential operator (First Di erential Approximation (FDA)) are considered. The transfer operators act almost as identities for very low-frequency harmonics and no reduction of amplitudes can be achieved by the Gauss-Seidel relaxation. Thus we obtain the following two-grid ampli cation factor () (see (3)), where (; h) denotes the eigenvalue which corresponds to the eigenfunction h (; x) from (27): (30) () := 1 ? (2(;; 2hh) ) = 0:75 for ! 2 and 2 = 1 ! 0 : For = ?1:99 and MG-FB, the limiting value of 0.75 is found by the predicted two-grid factor F and also by the numerically obtained tgh , see Table 4. The multigrid convergence mg h increases further since the coarse grid problem occurs on all coarser grids (actually 6 grids are used). Ways to handle the coarse grid problem successfully is to change the smoother to the more tg Table 4: Multi- and two-grid convergence mg h and h , predicted convergence F and smoothing factor F for the -problem, MG-FB, h = 1=64. = ?1:90 = ?1:99
mg h
tgh
F
mg h
F
tgh
F
F
0.747 0.654 0.662 0.499 0.839 0.730 0.745 0.518 expensive ILU-relaxation, so that the coarse grid problem is removed by a powerful smoother on ner grids [8], or to change the grid coarsening to nonstandard coarsening which takes also care of the problematic low-frequency error components. We keep the point Gauss-Seidel smoother and standard coarsening and solve this problem by the Krylov subspace acceleration of MGFB. In Table 5, we observe a signi cant improvement of the multigrid convergence due to the acceleration even with a small Krylov subspace. The accelerated two-grid convergence acc h (tg ) U E is very well predicted by the analytical estimations i20 and Tm . 9
Comparing the accelerated two- and multigrid convergence from Table 5 it is expected that it is possible to further improve acc h by incorporating the Krylov acceleration into the multigrid cycle and apply it also on the coarse grids. For example, if GMRES(2) is used on all grids the convergence rate improves from 0:533 to 0:470 for = ?1:99. acc Table 5: Accelerated multi- and two-grid convergence acc h and h (tg ), predicted convergence factor Ui20 and estimation TmE for the -problem, h = 1=64. solution = ?1:90 = ?1:99 acc acc U E acc acc method h h (tg) i20 Tm h h (tg) Ui20 TmE acc. MG-FB, m = 2 0.420 0.333 0.372 0.373 0.533 0.404 0.415 0.463 acc. MG-FB, m = 5 0.370 0.287 0.302 0.304 0.472 0.348 0.354 0.378
5 Conclusions In this paper we have presented a way to obtain sharp quantitative convergence estimations for GMRES(m) preconditioned by multigrid on the basis of Fourier analysis. For all the cases considered the estimations are accurate compared to measured numerical convergence of the multigrid preconditioned GMRES(m) method. It has been found that it is not easily possible to further accelerate multigrid methods which are optimally tuned with, for example, overrelaxation parameters. In other situations, however, very satisfactory convergence improvement is achieved with the Krylov subspace acceleration. For a fair comparison one should take the additional work into account, which is not the main issue here. The proposed acceleration is, of course, also applicable to situations in which it is not easy to tune a multigrid method.
References [1] A. Brandt, Multi{level adaptive solutions to boundary{value problems. Math. Comp., 31, 333{390 (1977). [2] A. Brandt, Rigorous quantitive analysis of multigrid, I: Constant coecients two-level cycle with L2 -norm. SIAM J. Numer. Anal., 31, 1695{1730 (1994). [3] Michael Holst and Stefan Vandewalle, Schwarz methods: To symmetrize or not to symmetrize. SIAM J. Numer. Anal., 34(2), 699-722 (1997). [4] Noel M. Nachtigal, Satish G. Reddy and Lloyd N. Trefethen, How fast are nonnsymmetric matrix iterations ? SIAM J. Matrix Anal. Appl., 13, 778{795 (1992). [5] C.W. Oosterlee and T. Washio, An evaluation of parallel multigrid as a solver and a preconditioner for singular perturbed problems. SIAM J. Sci. Comput., 19, 87{110 (1998) [6] Y. Saad and M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Comput., 7, 856{869 (1986). [7] Y. Saad, Iterative Methods for Sparse Linear Systems. PWS Publishing, Boston, USA (1996). [8] K. Stuben and U. Trottenberg, Multigrid methods: fundamental algorithms, model problem analysis and applications. In: W.Hackbusch, U. Trottenberg (eds.), Multigrid Methods, Lecture Notes in Math. 960, 1{176, Springer, Berlin Germany (1982). [9] C.A. Thole, U. Trottenberg, Basic smoothing procedures for the multigrid treatment of elliptic 3-D operators. Appl. Math. Comp., 19, 333{345 (1986). [10] I. Yavneh, On Red-Black smoothing in multigrid. SIAM J. Sci. Comput., 17, 180{192 (1996).
10
Abstract:
This paper deals with convergence estimations of a preconditioned GMRES(m) method [6], where multigrid ([1], [8]) is used as the preconditioner. Fourier analysis is a well-known useful tool in the multigrid community for the prediction of two-grid convergence rates ([1], [8]). This analysis is generalized to the situation where multigrid is a preconditioner.
1 Introduction Nowadays, it has become popular to study the convergence of multilevel methods and use them in combination with a Krylov subspace acceleration method. Recently, also standard (multiplicative) multigrid methods are used as a preconditioner. This application of standard multigrid is bene cial in situations where standard multigrid alone does not converge fully satisfactorily, because certain error frequencies are not reduced well enough, which mainly occurs when complicated PDEs are solved. In general, it is dicult to construct robust multigrid solvers for large classes of problems. Among many other papers, this approach has been presented in [5] for linear singularly perturbed problems. Here, we analyze one combined solution method, the multigrid preconditioned GMRES(m) solution method, theoretically. We derive quantitative convergence results in order to understand the convergence and the bene ts of the combination. The basis for the theoretical convergence estimations is Fourier analysis, which is a well-known tool for obtaining two-grid convergence factors [8]. In this paper, we restrict ourselves to standard problems coming from scalar linear PDEs discretized with nite di erences on uniform Cartesian grids and we keep the multigrid method simple. For example, we restrict ourselves to point smoothers and standard grid coarsening and analyze the e ect of Krylov subspace acceleration for anisotropic di usion problems and for problems with mixed derivatives, although it is known that other (more expensive) smoothers ore other coarsening strategies are appropriate remedies for such equations. However, in this way we try to gain insight in the behavior of the combination of multigrid and GMRES(m). The typical use of the Fourier analysis in a multigrid context is that the spectral radius is obtained theoretically, which is the basis for asymptotic multigrid convergence estimations [8]. Two-grid Fourier analysis transforms a multigrid iteration matrix into a block diagonal matrix from which more information about the spectrum is easily obtained. The analysis makes use of unitary basis transformations to simplify the representation of the multigrid iteration matrix. It is therefore appropriate for the estimation of the multigrid preconditioned GMRES(m) convergence since unitary transformations do not a ect the convergence behavior of GMRES [6]. A preconditioned Krylov subspace acceleration method like GMRES(m) implicitly builds up a minimal residual polynomial. The determination of the polynomial coecients is easily possible with the help of the Fourier analysis and can be done explicitly, since a block diagonal matrix results. Based on the GMRES(m) polynomial and the Fourier analysis, sharp theoretical convergence estimates can be obtained. In section 2 the multigrid preconditioner is presented. The two-grid Fourier analysis results are used in section 3, where di erent ways to analyze restarted GMRES(m) preconditioned by multigrid are discussed. There are actually two types of Fourier analysis available in the multigrid community. The rst type is the analysis called \rigorous analysis" here, since it explicitly takes boundary conditions into account and discrete eigenfunctions are related to the discrete mesh size. The second is known as local mode analysis [1]. Exponential functions are the basis for the analysis, a GMD, Institute for Algorithms and Scienti c Computing, D-53754 Sankt Augustin, Germany
(email: wienands @ gmd.de)
1
discretization at domain boundaries cannot be taken into account. This analysis can be used for a wider range of PDEs and multigrid components. The generalization of Fourier analysis to multigrid as a preconditioner is applied for both types of analysis. By several numerical tests the theoretical convergence estimations are compared to the actual numerical convergence. The multigrid purist may not like the approach presented here with multigrid viewed as a preconditioner. It is, of course, also possible to try to construct an optimal multigrid method for each separate problem. For the purist, the analysis might give an indication how far his present method is from an optimal multigrid solution method. In other words, if an additional acceleration of the convergence is found by the Fourier analysis, it must also be possible to tune one of the multigrid components, so that multigrid as a stand alone solver is improved.
2 Rigorous Fourier Analysis of Multigrid The 3D Poisson equation, a simple test problem, is chosen to explain the possibilities of convergence estimations by Fourier analysis for multigrid as a solver. We apply the standard 7-point discretization on a cube with Dirichlet boundary conditions and use a uniform mesh with meshsize h = 1=n:
Ah uh(x) = ?h uh(x) = bh (x) on h = (0; 1)3 \ Gh ; (1) uh(x) = 0 on ?h = ([0; 1]3 n h ) \ Gh with Gh := fx = hj ; j 2 Z3g. We combine operator and matrix notation. Multigrid is
commonly explained with operators (see, for example, [8], [1]), whereas Krylov subspace methods are often described with matrices. A well-known ecient 3D multigrid method [9] for problem (1) consists of a red-black GaussSeidel smoothing method (GS-RB), Full Weighting of residuals to obtain the right-hand side on the coarse grids, tri-linear interpolation of corrections from coarse to ne grids and the coarse grid discretization of the PDE on a grid with meshsize 2h in each direction. We call this method the red-black multigrid solver (MG-RB) here. For problem (1) and MG-RB, it is possible to apply rigorous two-grid Fourier analysis, as it was presented for 2D problems in [8] and used for 3D (Poisson-type) problems in [9]. The existence of the orthonormal basis of discrete eigenfunctions
'k;`;m h (x; y; z ) = sin kx sin `y sin mz with k; `; m = 1; : : :; n ? 1 ((x; y; z ) 2 h )
(2)
for the operator (1) is crucial for the 3D discrete Fourier analysis. They generate the space of all grid functions, F ( h ). The error vj ?1 is transformed by a two-grid cycle as follows: (3) vj = Mh2h vj?1 with Mh2h = Sh2 Kh2h Sh1 = Sh2 (Ih ? I2hh (2h)?1Ih2h h )Sh1 ; where Sh is the smoothing operator, 1 and 2 indicate the number o f pre- and post-smoothing iterations, 2h is the discretization of the operator on a 2h coarse grid, I2hh and Ih2h are transfer
operators from coarse to ne grids and reversed. In the case of standard multigrid coarsening (H = 2h) in 3D, it is convenient to divide F ( h ) into a direct sum of (at most) eight-dimensional subspaces, the 2h-harmonics Ehk;`;m, which leaves Mh2h invariant (see [9], [8]). Hence, Mh2h is f consisting of (at most) 8 8 blocks ([8], [9]): orthogonally equivalent to a block matrix M h
i
^ M 2h j c(k; `; m) = f := M c(k; `; m) =^ Mh2h with M M h Ehk;`;m : k;`;m=1;::: ; n 2
(4)
The tilde e indicates the Fourier matrix representation of the related operator. By M we denote the matrix representation of the multigrid operator Mh2h and by A the matrix representation of f = UMU ?1 with a unitary matrix U , which h with respect to the Euclidian base. We have M 2
rows are given by the eigenvectors of A (for the corresponding eigenfunctions see (2)). From the c(k; `; m), we obtain the two-grid convergence factor by block matrices M f) = c(k; `; m)) : F := (M max (M 1k;`;m n 2
(5)
It is the asymptotic two-grid convergence which is a theoretical estimation of the multigrid convergence. It is furthermore possible to determine the whole spectrum and thus the eigenvalue f by calculating the eigenvalues of the block matrices M c(k; `; m). We use this distribution of M distribution in section 3 for the analysis of the convergence of the multigrid preconditioned GMRES(m) method.
3 Fourier analysis of multigrid preconditioned GMRES(m) In this section we discuss two approaches for the analysis of multigrid preconditioned GMRES(m). For both approaches it is crucial that the discrete Fourier analysis makes use of the unitary basis transformations which simplify the representation of the multigrid iteration matrix to the block matrix (4). We consider the linear system related to (1), Au = b. A two-grid (or more general a multigrid) cycle can be described by the following matrix splitting: Cuj + (A ? C )uj?1 = b ; (6) where uj represents a new and uj ?1 a previous approximation. The multigrid method (6) considered here is described in the previous section. It is used as a right preconditioner for GMRES(m). The parameter m represents the number of vectors stored after which the GMRES(m) algorithm will restart. GMRES(m) searches for a solution uj in the following subspace: C (uj ? uj?m ) 2 span[rj?m ; (AC ?1)rj?m ; :::; (AC ?1)m?1 rj?m ] =: K m(AC ?1 ; rj?m ) ; (7) where K m (AC ?1; rj ?m ) is the Krylov subspace. It selects the new solution vector uj by minimizing the corresponding residual rj in the L2-norm minm ?1 jjb ? Auj jj2 : (8) C (uj ?uj?m )2K (AC ;rj?m ) the ane subspace uj ?m + C ?1 K m (AC ?1 ; rj ?m) can be
Any element w in represented by w = uj?m + C ?1 ( 1rj?m + 2 AC ?1 rj?m + + m (AC ?1 )m?1 rj?m ): (9) By substituting (9) into the residual equation, we obtain rw = b ? Aw = Pm (AC ?1)rj?m ; (10) P
where Pm is a m-th order polynomial de ned by Pm () = 1 ? mk=1 k k . Comparing (8) and (10), we see that the L2-norm of rj satis es the following minimization property: jjrj jj2 = minfjjPm(AC ?1)rj?mjj2 j Pm 2 Pmg; (11) where Pm denotes the set of all polynomials of degree at most m with Pm (0) = 1. Since we are interested in the convergence behavior of multigrid preconditioned GMRES with a restart after m iterations we consider the residuals rm , r2m, ..., rim . This leads to the following de nition of the reduction factors i for a complete iteration, consisting of m multigrid preconditioned GMRES(m) steps: 1=mi1=i h : i := jjjjrrimjjjj2 0 2
3
(12)
In section 3.1 we review estimations for i which are based on an analysis of the spectrum of AC ?1 (see [6], [7]). In general, these estimations exhibit a qualitative character and concrete values are rarely found in the literature because the spectrum is usually not available. However, with multigrid as the preconditioner these estimations can be calculated explicitly since the spectrum is easily obtained. In section 3.2 we will introduce a sharper estimation, which is based on an analysis of the GMRES(m) polynomial.
3.1 Analysis with the spectrum of the iteration matrix
The common way to analyze the convergence of GMRES is to exploit information about the spectrum of the iteration matrix AC ?1 . If we compare the error transformation (3) and the f (4) is related to the two-grid residual transformation by one two-grid cycle, the Fourier matrix M ?1 iteration matrix I ? AC via the identity f = UA?1 (I ? AC ?1 )AU ?1 () M f = I ? UC ?1 AU ?1 ; M with the unitary matrix U from (4). With (13) we have
(13)
f(UAU ?1 )?1 = I ? A eM fA e?1 : UAC ?1 U ?1 = I ? UAU ?1 M
(14)
Ae and Ae?1 are diagonal matrices since the rows of U are eigenvectors of A. This means that UAC ?1 U ?1 is a block-diagonal matrix which consists of (at most) 8 8 blocks in the 3D case and of (at most) 4 4 blocks in the 2D case. The eigenvalues of these block matrices are easily calculated. First, we discuss the considerations on the convergence of GMRES(m) from [6] and [7], which are based on the spectrum of the iteration matrix. From (11) it follows that
jjrimjj2 jjPm(AC ?1)jj2 jjr(i?1)mjj2 for all Pm 2 Pm : (15) Suppose that all eigenvalues of AC ?1 are located in an ellipse E (c; d; a) which excludes the origin. c denotes the center, d the focal distance and a the major semi-axis. Then, it is known [7] that the absolute value of the polynomial tm (z) := Tm ( dc ? 1d z)=Tm( dc ) = Tm(zb)=Tm( dc ) with z; zb := ( dc ? d1 z) 2 C
(16)
is small on the spectrum of AC ?1 . Here, Tm represents the Chebychev polynomial of degree m of the rst kind which is given by the three-term recurrence relation
T0(Ab) = I ; T1 (Ab) = Ab ; Tm+1 (Ab) = 2Tm(Ab) ? Tm?1 (Ab) for m 1 : If AC ?1 is diagonalizable by AC ?1 = XDX ?1 then (15) yields
jjrimjj2 jjtm(AC ?1)jj2 jjr(i?1)mjj2 ii h 2(X ) Tm( da )=Tm( dc ) jjr0jj2
h
ii
jjtm(AC ?1)jj2 jjr0jj2
(17) (18) (19)
where 2 (X ) denotes the spectral condition number of the transformation matrix X [7]. Inequality (19) uses the fact that the maximum modulus of a complex analytical function is reached on the boundary of the domain which is due to Liouville's theorem [7]. Using the above relations (18) and (19) we obtain the following well known estimations ([7]) for i
i NmE (2 (X ))1=m TmE 1=m : with NmE := (jjtm(AC ?1 )jj2)1=m and TmE := Tm ( da )=Tm ( dc ) 4
(20) (21)
TmE can be further approximated by TmE
=
p
+ ( ad )2 ? 1 m + da + ( ad )2 ? 1 ?m 1=m a + a2 ? d2 p2 2: t p c 2 ? 1m + ? c + p( c )2 ? 1?m c + c ?d + ( ) d d d d
?a d ?c
p
?
p
(22)
The approximation from (22) is sharp if a d and c d or m 1 holds. The rst inequality is ful lled if the ellipse E (c; d; a) tends to a circle which means that the focal distance d becomes very small. In all cases considered, the multigrid preconditioner ensures that the spectrum is clustered around 1, which means that the second inequality is satis ed. If (22) is sharp due to a circular shape of or due to a large m, TmE is independent of the restart parameter m. In this case, it is dicult to accelerate the multigrid preconditioner. A large Krylov subspace (7) does not lead to additional acceleration. The main drawback of the estimations (20) for the reduction factors i is their i-independence. Thus, they can not expected to be sharp approximations for every i. Furthermore the constant 2 (X ) can be very large if the iteration matrix AC ?1 is far from normal but it is found in [4] that the condition number 2 (X ) is not always relevant for the convergence of GMRES. In [7] it is stated that (20) is an asymptotic result and that the actual residual norm should behave like TmE without 2 (X ). Therefore, explicit values also for TmE are quoted in section 3.3 and compared with asymptotic numerical convergence results. However, TmE is a heuristic estimation for i 1 and not an upper bound like NmE .
3.2 Analysis with the GMRES-polynomial
To obtain a quantitative result for the convergence of GMRES(m), we present a second approach to estimate i, which explicitly depends on the iteration index i in contrast to the previous section. An initial residual r0 is prescribed and the GMRES(m)-polynomials Pmi are explicitly determined for every i. The choice of r0 does not in uence the average reduction factor for i >> 1. For any unitary matrix U , in particular for U from (4), we have jjrimjj2 = jjUrimjj2. In order to nd the coecients ik (k = 1; :::; m) of the GMRES(m)-polynomials Pmi , the following function f is minimized:
f ( i1 ; : : :; im) := UPmi (AC ?1 )r(i?1)m; UPmi (AC ?1 )r(i?1)m =
m X
(23)
l;k=1
+2
il ik (UAC ?1 U ?1 )lUr(i?1)m; (UAC ?1U ?1 )k Ur(i?1)m
m X l=1
il (UAC ?1 U ?1 )lUr(i?1)m ; Ur(i?1)m :
(24)
The ik are obtained by solving the following linear system of equations: m @f =2 X ?1U ?1 )k Ur i (UAC ?1 U ?1 )l Ur ; ( UAC (i?1)m (i?1)m @ il k=1 k + 2 (UAC ?1 U ?1 )l Ur(i?1)m; Ur(i?1)m = 0 for l = 1; : : :; m :
(25)
fA e?1 of the GMRES(m) polynomials In (14) it is seen that the argument UAC ?1 U ?1 = I ? AeM Pmi has a simple block structure. Therefore, the solution of the linear system (25) can be easily computed for every complete iteration i due to the sparse structure of (UAC ?1U ?1 )l (l = 1; :::; m), if the previous transformed residual Ur(i?1)m is given. We prescribe right-hand side b = 0 and a randomly chosen initial approximation Uu0, which yield the transformed initial residual
5
P
n?1 (k; `; m)'k;`;m e 0 . More precisely, we select an initial solution u0 = Ur0 = ?AUu k;`;m=1 where the amplitudes (k; `; m) of the eigenvectors are randomly chosen. Hence we get Uu0 = ( (1; 1; 1); (n ? 1; n ? 1; n ? 1); : : :; ( n2 ; n2 ; n2 ) )T , see (2). From Ur0 it is possible to evaluate the coecients 1k by solving (25) for i = 1. This gives Urm = Pm1 (UAC ?1U ?1 )Ur0. It follows that Urim can be easily calculated for every i. This leads to the de nition of an average reduction factor Ui , which can be obtained by the Fourier analysis: h jjUr jj 1=m i1=i U : (26) i := jjUrimjj 2 0 2 In most of the tests, presented in the next subsection, Ui becomes constant for i 20. Thus, U20
is expected to match very well with numerical reference values and in particular to be sharper than the estimation from section 3.1.
3.3 Fourier results of multigrid preconditioned GMRES(m)
Here, we investigate the quality of the theoretical considerations by applying the multigrid preconditioned GMRES(m) method to the 3D Poisson problem. We compare the average accelerated multigrid convergence rate acc h from numerical calculations with the theoretical prediction U E i20, the heuristic estimation Tm and the upper bound NmE . Insight in the behavior of the acceleration method is gained by considering di erent combinations of the relaxation parameter ! and the GMRES(m) restart parameter m. In [10] it is presented that an overrelaxation parameter (! > 1) improves the smoothing properties of GS-RB (leading to ! -GS-RB) and therefore improves the convergence of MG-RB for d-dimensional Poisson-type equations. We investigate if it is possible to further accelerate the multigrid methods with these relaxation parameters by GMRES(m). The results are presented in Table 1. We observe that Ui20 and TmE provide accurate predictions of the multigrid preconditioned GMRES(m) convergence acc h for di erent E overrelaxation parameters ! . The upper bound Nm becomes sharper with increasing m as TmE and NmE coincide for large m, because the condition number 2 (X ) cancels out (see (20)). Table 1: Multigrid convergence mg convergence acc h , predicted converh , acceleratedE multigrid U E gence F , i20 and convergence estimations Tm , Nm from Chebychev polynomials for the 3D Poisson equation; W(1,1), h = 1=32. Relaxation MG-RB acc. MG-RB, m = 2 acc. MG-RB, m = 5 acc U E E acc parameter mg T N Ui20 TmE NmE F m m i20 h h h ! = 1:0 0.192 0.194 0.085 0.086 0.088 0.311 0.070 0.074 0.085 0.149 ! = 1:1 0.089 0.097 0.045 0.056 0.059 0.272 0.042 0.049 0.059 0.102 ! = 1:15 [10] 0.070 0.077 0.050 0.057 0.060 0.234 0.045 0.048 0.060 0.102 In Table 1, we observe a satisfactory convergence improvement for ! = 1. If overrelaxation is selected (! = 1:1; 1:15) the bene ts of the GMRES(m) acceleration are small and the convergence cannot be improved by using a larger Krylov space (7). This behavior is expected from the spectral picture in Figure 1. The ellipses become orbital (d a) with the overrelaxation so that (22) holds by which a m-independent acceleration is predicted.
4 Local mode analysis of multigrid preconditioned GMRES(m) The previous Fourier analysis is restricted to model operators Ah , for which the discrete functions (2) form a basis of eigenfunctions. Furthermore, it is not possible to use the analysis for multigrid methods based on Gauss-Seidel relaxations with a lexicographical ordering of the grid points (GSLEX). The eigenfunctions (2) are no longer invariant under the GS-LEX operators. Instead they 6
ellipse E(c=0.903,d=0.078,a=0.097) spectrum
0.04
0.04
0.02
0.02
0
0
-0.02
-0.02
-0.04
-0.04
-0.06
-0.06
0.8
0.85
0.9 REAL PART
0.95
ellipse E(c=0.959,d=0.005,a=0.056) spectrum
0.06
IMAGINARY PART
IMAGINARY PART
0.06
1
0.9
0.92
0.94
0.96 REAL PART
0.98
1
1.02
Figure 2-a: ! = 1:00 Figure 2-b: ! = 1:10 Figure 1: Eigenvalue spectra of the Red-Black Two-Grid 3D Poisson solver from rigorous Fourier analysis.(Note the di erent scaling on the real- and imaginary-axis in Fig. 2-a.) are intermixed, which leads to a full matrix Se. For a general linear operator Ah with constant or frozen coecients it is, however, possible to use the local mode Fourier analysis [1] for analyzing multigrid methods based on GS-LEX as the smoother. The local mode (or in nite grid) analysis gives satisfactory estimations of the actual multigrid convergence. Its theoretical justi cation is, however, not as straightforward as for the previous analysis [2]. Instead of considering the nite domain with discrete eigenfunctions (2) of Ah an in nite domain h := fx = (jx h; jy h; jz h) : jx ; jy; jz 2 Zg is considered with continuous eigenfunctions h (; x) = eix=h = eij = ei(jx x +jy y +jz z ) ; (27) where the Fourier frequencies = (x ; y ; z ) vary continuously in = (?; ]3. The Fourier space "h = spanfeij : 2 g can be divided into eight-dimensional subspaces, the 2hharmonics "h , which leaves Mh2h invariant [8]. In order to get a well-de ned two-grid operator we replace "h by a slightly shrunk subspace "~h = spanfeij : 2 n g such that (Ah )?1 exists and also (A2h )?1 can be reasonably de ned on the coarser grid, like in [8]. The spectral radius f) of the corresponding block diagonal Fourier two-grid matrix is determined similarly as in (M the previous Fourier analysis (5): ^ M 2h j : f) = sup (M c()) with ~ l = (?=2; =2]3n and M c() = F := (M (28) h "h ~l 2
The whole spectrum is again obtained similarly as in section 2. We, however, always observe a spectrum from a discrete mesh, whereas we have a continuously de ned set of eigenfunctions (27) in order to get h-independent upper bounds for the convergence rates. With the local mode analysis the analytical values Ui20 and TmE , NmE for the accelerated multigrid method can be obtained by straightforward modi cations of the de nitions (26) and (21) taken from the rigorous case.
4.1 Results for Poisson-type equations
The rst example for the local mode Fourier analysis of multigrid preconditioned GMRES(m) deals with the anisotropic di usion equation. We analyze the di erence in convergence between 7
performing two smoothing iterations (1 = 2 = 1) of the forward lexicographical point GaussSeidel (GS-fLEX) method and one iteration of forward lexicographical point Gauss-Seidel and one iteration of backward lexicographical point Gauss-Seidel (GS-bLEX). (The other multigrid components are identical to the previously introduced multigrid method MG-RB.) We call these multigrid methods MG-FF and MG-FB, respectively, in the following. Local mode analysis yields sharp quantitative predictions for the multigrid convergence of MG-FF and MG-FB for the 2D anisotropic Poisson equation, as can be seen in Table 2. In the anisotropic case it is bene cial to introduce an overrelaxation parameter also for these lexicographical smoothers leading to ! -MG-FF and ! -MG-FB. (The improvement is less impressive than for ! -MG-RB in [10], even if the optimal overrelaxation parameter (!opt = 1:40 for " = 0:1, !opt = 1:75 or " = 0:01) is selected.) Table 2: Multigrid convergence mg h and predicted convergence F for the 2D anisotropic Poisson equation; h = 1=128. solution " = 0:1, ! = 1:00 " = 0:1, ! = 1:40 " = 0:01, ! = 1:00 " = 0:01, ! = 1:75 method mg F mg F mg F mg F h h h h MG-FF 0.693 0.696 0.410 0.433 0.957 0.961 0.755 0.758 MG-FB 0.693 0.697 0.437 0.440 0.957 0.962 0.750 0.759 However, in contrast to !opt -MG-FF (and !opt-MG-RB), it is possible to further accelerate !opt-MG-FB signi cantly with GMRES(m). This can be observed by comparing Table 2 and Table 3 and can be explained as follows. MG-FB results in a symmetric iteration matrix due to U Table 3: Accelerated multigrid convergence acc h , predicted convergence i20 and estimations TmE , NmE from Chebychev polynomials for the 2D anisotropic Poisson equation; h = 1=128. solution " = 0:1, ! = 1:00 " = 0:1, ! = 1:40 acc U E E acc method h i20 Tm Nm h Ui20 TmE NmE acc. MG-FF, m = 2 0.475 0.515 0.529 0.911 0.410 0.420 0.430 0.524 acc. MG-FF, m = 5 0.450 0.469 0.519 0.907 0.410 0.420 0.430 0.441 acc. MG-FB, m = 2 0.370 0.402 0.408 0.718 0.185 0.185 0.203 0.445 acc. MG-FB, m = 5 0.305 0.331 0.331 0.487 0.158 0.160 0.165 0.202 solution " = 0:01, ! = 1:00 " = 0:01, ! = 1:75 U E E acc U method acc T N TmE NmE m m i20 i20 h h acc. MG-FF, m = 2 0.850 0.895 0.907 0.952 0.753 0.750 0.757 0.923 acc. MG-FF, m = 5 0.800 0.837 0.896 0.897 0.753 0.750 0.757 0.813 acc. MG-FB, m = 2 0.784 0.850 0.865 1.403 0.397 0.380 0.479 0.717 acc. MG-FB, m = 5 0.723 0.755 0.768 1.102 0.362 0.370 0.392 0.457
the construction of the coarse grid correction (FW as the restriction and bilinear interpolation as the prolongation operator) and the two relaxations which are adjoined to each other. This means that the spectrum of the iteration matrix is real-valued, and uniformly distributed on the interval [min ; max]. This feature is maintained even if !opt is introduced, whereas for !opt -MGFF an acceleration by GMRES(m) is not possible, because the spectrum becomes orbital. It is found in [3] that a symmetrization is not bene cial if Bi-CGstab is used as the Krylov subspace acceleration method. The analytical prediction Ui20 shows a sharp quantitative character again and also the estimation TmE gives a reliable indication for the average accelerated convergence rate presented in Table 3. NmE is an interesting quantity for large m. Remark: For symmetric matrices it is possible to use the Conjugate Gradient (CG) iteration as acceleration method. The asymptotic convergence rate cg for the preconditioned CG iteration 8
is well-known. If we consider the de nition TmE (20) for a real-valued spectrum [min ; max] then the ellipse E (c; d; a) deteriorates to a line. In particular, we have a = d = (max ? min )=2 and c = min + a, which leads to TmE = (Tm ( aa )=Tm ( ac ))1=m = cg (see (21), (22)). Thus it can be expected heuristically that GMRES(m) has the same asymptotic convergence as CG. Remark: Also for the 3D Poisson equation the multigrid convergence of the two methods MG-FF and MG-FB is similar. The theoretically obtained convergence factor F is 0.321 for both methods, although the spectra are completely di erent (as in 2D). Very satisfactory accelerated convergence rates are obtained by the Krylov subspace acceleration of MG-FB. Due to the symmetry we can predict the convergence rate of the accelerated MG-FB accurately from cg .
4.2 Results for a problem with mixed derivative
Finally, we discuss an equation where reasonable smoothing properties can be achieved with point smoothers, but the coarse grid correction turns out to be problematic in combination with standard coarsening [8]. Again an acceleration with the help of GMRES(m) can be applied as a remedy to overcome this diculty. We investigate the 2D di erential equation ?u ? uxy = b with a mixed derivative discretized via the O(h2 ) 9-point operator 2 3 4 ?1 ? 4 1 Ah uh (x) = h2 4 ?1 4 ?1 5 uh (x) = bh (x) on h : (29) ? 4 ?1 4 The di erential equation is no longer elliptic if j j = 2 and the eciency of all previously considered multigrid methods (MG-RB, MG-FF, MG-FB) deteriorates for j j ! 2 (see for example Table 4 and [8]). This behavior can be explained if we perform a simpli ed two-grid analysis where only the very low Fourier frequencies along the characteristic direction of the di erential operator (First Di erential Approximation (FDA)) are considered. The transfer operators act almost as identities for very low-frequency harmonics and no reduction of amplitudes can be achieved by the Gauss-Seidel relaxation. Thus we obtain the following two-grid ampli cation factor () (see (3)), where (; h) denotes the eigenvalue which corresponds to the eigenfunction h (; x) from (27): (30) () := 1 ? (2(;; 2hh) ) = 0:75 for ! 2 and 2 = 1 ! 0 : For = ?1:99 and MG-FB, the limiting value of 0.75 is found by the predicted two-grid factor F and also by the numerically obtained tgh , see Table 4. The multigrid convergence mg h increases further since the coarse grid problem occurs on all coarser grids (actually 6 grids are used). Ways to handle the coarse grid problem successfully is to change the smoother to the more tg Table 4: Multi- and two-grid convergence mg h and h , predicted convergence F and smoothing factor F for the -problem, MG-FB, h = 1=64. = ?1:90 = ?1:99
mg h
tgh
F
mg h
F
tgh
F
F
0.747 0.654 0.662 0.499 0.839 0.730 0.745 0.518 expensive ILU-relaxation, so that the coarse grid problem is removed by a powerful smoother on ner grids [8], or to change the grid coarsening to nonstandard coarsening which takes also care of the problematic low-frequency error components. We keep the point Gauss-Seidel smoother and standard coarsening and solve this problem by the Krylov subspace acceleration of MGFB. In Table 5, we observe a signi cant improvement of the multigrid convergence due to the acceleration even with a small Krylov subspace. The accelerated two-grid convergence acc h (tg ) U E is very well predicted by the analytical estimations i20 and Tm . 9
Comparing the accelerated two- and multigrid convergence from Table 5 it is expected that it is possible to further improve acc h by incorporating the Krylov acceleration into the multigrid cycle and apply it also on the coarse grids. For example, if GMRES(2) is used on all grids the convergence rate improves from 0:533 to 0:470 for = ?1:99. acc Table 5: Accelerated multi- and two-grid convergence acc h and h (tg ), predicted convergence factor Ui20 and estimation TmE for the -problem, h = 1=64. solution = ?1:90 = ?1:99 acc acc U E acc acc method h h (tg) i20 Tm h h (tg) Ui20 TmE acc. MG-FB, m = 2 0.420 0.333 0.372 0.373 0.533 0.404 0.415 0.463 acc. MG-FB, m = 5 0.370 0.287 0.302 0.304 0.472 0.348 0.354 0.378
5 Conclusions In this paper we have presented a way to obtain sharp quantitative convergence estimations for GMRES(m) preconditioned by multigrid on the basis of Fourier analysis. For all the cases considered the estimations are accurate compared to measured numerical convergence of the multigrid preconditioned GMRES(m) method. It has been found that it is not easily possible to further accelerate multigrid methods which are optimally tuned with, for example, overrelaxation parameters. In other situations, however, very satisfactory convergence improvement is achieved with the Krylov subspace acceleration. For a fair comparison one should take the additional work into account, which is not the main issue here. The proposed acceleration is, of course, also applicable to situations in which it is not easy to tune a multigrid method.
References [1] A. Brandt, Multi{level adaptive solutions to boundary{value problems. Math. Comp., 31, 333{390 (1977). [2] A. Brandt, Rigorous quantitive analysis of multigrid, I: Constant coecients two-level cycle with L2 -norm. SIAM J. Numer. Anal., 31, 1695{1730 (1994). [3] Michael Holst and Stefan Vandewalle, Schwarz methods: To symmetrize or not to symmetrize. SIAM J. Numer. Anal., 34(2), 699-722 (1997). [4] Noel M. Nachtigal, Satish G. Reddy and Lloyd N. Trefethen, How fast are nonnsymmetric matrix iterations ? SIAM J. Matrix Anal. Appl., 13, 778{795 (1992). [5] C.W. Oosterlee and T. Washio, An evaluation of parallel multigrid as a solver and a preconditioner for singular perturbed problems. SIAM J. Sci. Comput., 19, 87{110 (1998) [6] Y. Saad and M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Comput., 7, 856{869 (1986). [7] Y. Saad, Iterative Methods for Sparse Linear Systems. PWS Publishing, Boston, USA (1996). [8] K. Stuben and U. Trottenberg, Multigrid methods: fundamental algorithms, model problem analysis and applications. In: W.Hackbusch, U. Trottenberg (eds.), Multigrid Methods, Lecture Notes in Math. 960, 1{176, Springer, Berlin Germany (1982). [9] C.A. Thole, U. Trottenberg, Basic smoothing procedures for the multigrid treatment of elliptic 3-D operators. Appl. Math. Comp., 19, 333{345 (1986). [10] I. Yavneh, On Red-Black smoothing in multigrid. SIAM J. Sci. Comput., 17, 180{192 (1996).
10
