[PPT] - AN INTRODUCTION TO MULTIGRID METHODS VIA SUBSPACE CORRECTION PowerPoint Presentation

SLIDE 1

AN INTRODUCTION TO MULTIGRID METHODS VIA SUBSPACE CORRECTION FRAMEWORK

LUDMIL ZIKATANOV, DEPARTMENT OF MATHEMATICS, PENN STATE

multigrid convergence.

CONTENTS 1. Introduction 2 2. An introduction to linear and product iterative methods 2 2.1. On the choice of B 3 2.2. Error equation and error transfer operator 3 2.3. Symmetrization of B and convergence 3 2.4. Gauss-Seidel method 4 2.5. Convergence rate 4 2.6. Convergence of Gauss-Seidel (cont) 4 2.7. Convergence – continued 4 2.8. Simple application 5 2.9. Preconditioning and convergent iterations 6 2.10. Some more on linear Iterative Methods 6 2.11. Errors in terms of product of projections 6 3. Subspace corrections 6 3.1. Some idea behind these methods 6 3.2. Variational problem and some notation 7 3.3. Subspace correction method: Algorithm 7 3.4. Error equations 7 3.5. Relationship with the Method of Alternating Projections 8 3.6. The method of alternating projections 8

SLIDE 2

4. Convergence of subspace correction methods 9 4.1. Convergence result (preliminaries) 9 4.2. Abstract theorem 10 4.3. Outline of the proof 10 5. How to estimate c0 12 5.1. A special decomposition 12 5.2. A multigrid method convergence proof 13 5.3. More general elliptic PDEs 15 5.4. Discretization and subspace splitting 16 5.5. Basic multigrid iteration 16 5.6. Basic facts from Finite element theory 16 5.7. Error estimates 17 5.8. On the addtive methods 17 5.9. Multiplicative overlapping DD-method 18 5.10. Convergence for finite element equations 18 5.11. Convergence via subspace correction framework 18 5.12. Further simplifications: Richardson smoother 19 5.13. Some known results 19 5.14. Stability and convergence 20 Basic references 20

(1) Au = f,

(2) uk+1 = uk + B(f − Auk), where B ∈ Rn×n is some given matrix (an approximate inverse of A).

SLIDE 3

2.1. On the choice of B.

(3) ek := u − uk, lim

(u, v)A := Au, vℓ2, u2

Let us take u to be the solution to (1). Then by f − Au = 0 we get that (4) u = u + B(f − Au), 2.2. Error equation and error transfer operator.

uk+1 = uk + B(f − Auk), u = u + B(f − Au) it follows that ek = u − uk = (I − BA)ek−1 = . . . = (I − BA)ke0.

ekA = (I − BA)ke0A ≤ I − BAk

(5) v2

BAv, vA, where (6) ¯ B := Bt + B − BtAB. 2.3. Symmetrization of B and convergence.

Gv2

BAv, vA, and we will have convergence, as long as ¯ BA, is positive definite (note that this matrix is symmetric in ·, ·A inner product. ( ¯ BAv, w)A = (A ¯ BAv, w) = (Av, ¯ BAw) = (v, ¯ BAw)A

SLIDE 4

2.4. Gauss-Seidel method.

A.

B for this choice of B: ¯ B := Bt(B−1 + B−t − A)B = Bt(D − L + D − Lt − D + L + Lt)B = BtDB.

B is symmetric and positive definite. and as excersize one can prove that GA < 1, just based on this fact. 2.5. Convergence rate.

G2

(I − BA)2

v2

B is invertible. Let us compute ¯ B−1: ¯ B−1 = B−1D−1B−T = (D − L)D−1(D − LT) = A + LD−1LT. 2.6. Convergence of Gauss-Seidel (cont).

B is below: ¯ B = (A + LD−1LT)−1.

(I − BA)v2

BAv, vA. 2.7. Convergence – continued.

(I − BA)2

¯ BAv, vA v, vA = 1 − inf

A ¯ BAv, v Av, v

B and change of variables w = A1/2v then can be used to obtain that: (I − BA)2

A1/2(A + LD−1LT)−1A1/2w, w w, w .

SLIDE 5

(I − BA)2

1 1 + supv=0

.

we get that G2

1 1 + c0 , and that is exactly the contraction number in A norm of the Gauss-Seidel iteration. In applications one needs to give an upper bound on c0 in order to get a convergence result. 2.8. Simple application.

LD−1LTv, v = 1 h2

v2

c0 = sup

N−1

h2v, vA .

N−1

h2v, vA

boundary point being shifted one h to the right. Hence: c0 = 1 4 sin2

.

c0 ∼ N 2 which gives the desired result.

SLIDE 6

2.9. Preconditioning and convergent iterations.

any eigelue λ(BA) we have |λ(I − BA)| < 1 ⇒ 1 − δ < |λ(BA)| < 1 + δ and hence B is a good preconditioner.

A =   2 1 1 1 2 1 1 1 2   , B−1 = diag(A) = 2I.

2.10. Some more on linear Iterative Methods. x ← x + Br; r := Ae = b − Ax.

2.11. Errors in terms of product of projections.

Pkv = αk(v)φk, αk = (Aφk, v)/φk2

then ⋆ Jacobi method is an example of additive preconditioner: B = Pk ⋆ Gaus-Seidel and SOR method are product iterative methods, with error transfer operators given by (for Gauss-Seidel ω = 1): ek = Eek, E = (I − ωP1) . . . (I − ωPn).

3.1. Some idea behind these methods.

⋆ Split the space into sum of subspaces. ⋆ Solve smaller problems on each of the subspaces. ⋆ “Glue up” the solutions of the subspace problems together to get an approximation of the solution in the entire space.

SLIDE 7

3.2. Variational problem and some notation.

a(u, v) = f(v), ∀v ∈ V.

ai(Tiv, vi) = a(v, vi), v ∈ V, vi ∈ Vi. with Ti = Pi if ai = a. 3.3. Subspace correction method: Algorithm.

Let u0 ∈ V be given. for ℓ = 1, 2, . . . uℓ−1 = uℓ−1. for i = 1 : J Let ei ∈ Vi solve ai(ei, vi) = f(vi) − a(uℓ−1

∀vi ∈ Vi. uℓ−1

= uℓ−1

endfor uℓ = uℓ−1

. endfor 3.4. Error equations.

= (I − Ti)(u − uℓ−1

u − uℓ = E(u − uℓ−1) = . . . = Eℓ(u − u0) where E = (I − TJ)(I − TJ−1) . . . (I − T1).

ai([Tiv], vi) = a(v, vi) ∀vi ∈ Vi.

SLIDE 8

3.5. Relationship with the Method of Alternating Projections.

E = (I − PJ)(I − PJ−1) . . . (I − P1).

3.6. The method of alternating projections. Let H be a Hilbert space and Mi ⊂ H be closed subspaces

PM2PM1 is a projection if and only if PM2PM1 = PM1PM2, furthermore PM2PM1 = PM1∩M2.

For general closed subspace M1 and M2 lim

✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✰ ✑✑✑✑✑✑✑✑✑✑✑✑✑✑✑ ✸

M1

✛ ✲

M2

✇

w0 = w

❏ ❏ ❏ ❏ ❏ ❏✇

w1

✇

w2

❏ ❏ ❏ ❏ ❏ ❏ ✇

w3

✇

w4

❏ ❏ ❏ ❏ ✇

w5

✇

w6

❏ ❏ ❏ ✇

w7

✇

w8

❏ ❏ ✇

w9

SLIDE 9

(PM2PM1)k − PM1∩M2 = [c(M1, M2)]2k−1 where c(M1, M2) is the cosine of the angle between M1 and M2: c(M1, M2) = sup |(u, v)| uv : u ∈ M1 ∩ (M1 ∩ M2)⊥, v ∈ M2 ∩ (M1 ∩ M2)⊥

⇒ uniform convergence.

4.1. Convergence result (preliminaries).

v ∈ V for some constant ω ∈ (0, 2). LEMMA. Assume that Ti satisfies the above assumption and that Range(Ti) = Vi. Then (1) I − Ti is nonexpansive. (2) Ti, T ∗

Ti have the same kernels: N( ¯ T i) = N(Ti) = N(T ∗

(3) Ti, T ∗

Ti have the same range: R( ¯ T i) = R(Ti) = R(T ∗

(4) The following inequality holds: 2 − ω ω Tiv2 ≤ v2 − (I − Ti)v2 = ( ¯ Tiv, v). (5) As operators restricted on Vi, the above (1)–(3) are still valid. (6) Ti, T ∗

Ti are all isomorphisms from Vi to itself. (7) ¯ Ti is nonnegative on V and symmetric positive definite on Vi. PROOF. It is obvious that (1) holds. By (1) I − Ti ≤ 1, namely Tiv2 ≤ 2(Tiv, v) = 2(v, T ∗

N(T ∗

Ti that N(Ti) ⊂ N( ¯ T i). Note that by the assumptions on the subspace solvers ( ¯ Tiv, v) = v2 − (I − Ti)v2 = 2(Tiv, v) − Tiv2 ≥ 2 − ω ω Tiv2 which implies N(Ti) ⊃ N( ¯ T i). This completes the proof of (2). (3) follows from (2) by the well-know relations between kernel and range of

SLIDE 10

4.2. Abstract theorem. Theorem (Xu and Z. 2003, J. of AMS) E2 ≡ (I − TJ)(I − TJ−1) . . . (I − T1)2 = c0 1 + c0 where ¯ Ti = T ∗

c0 = sup

inf

(Ti ¯ T −1

T ∗

vj − T −1

vi. Assumption: Tiv2 ≤ ω(Tiv, v), ω ∈ (0, 2) = ⇒ c0 > 0. If Ti = Pi, then c0 = sup

inf

Piwi2

vj. 4.3. Outline of the proof. The proof is based on a sequence of identities (1) Let E0 = I and Ei = (I − Ti)Ei−1 for i = 1 : J, then v2 − Ev2 =

( ¯ T iEi−1v, Ei−1v) = (I

T

where ¯ T

T1, ¯ T2, . . . , ¯ TJ), I

     I I . . . I      , E

     I E1 . . . EJ−1      , L

       I . . . T1 I . . . T1 T2 I . . . . . . . . . . . . ... . . . T1 T2 T3 . . . I        (2) The following identity holds v2 − Ev2 = (T

= ⇒ E2 = 1 − inf

SLIDE 11

where S

T

T

T

T

     T ∗

T ∗

. . . T ∗

     : V → V . In fact S

(T

T

T

= T

T

T

T

T

hence I

T

I

T

I

T

as desired (3) One key observation: for any w ∈ V, ˜ w = (S

v = T

w: (S

w, ˜ φ) = 0, ∀˜ φ ∈ N(T

⇒ (S

w, ˜ w) = inf

v, ˜ v). Indeed, (S

w, ˜ φ) = ((S

w, ˜ φ) = (T

φ) − (T

w, ˜ φ) = 0 Hence for arbitrary ˜ v, we have (S

v, ˜ v) = (S

w, ˜ w) + (S

w − ˜ v, ˜ w − ˜ v), because ˜ w − ˜ v ∈∈ N(T

SLIDE 12

Thus

−1 = sup

(T

(T

= sup

inf

(S

v, ˜ v) v2 + 1 = sup

inf

(S

v, ˜ v) v2 + 1 ≡ c0 + 1.

(7) (T ∗

T −1

T ∗

− I) = ¯ T −1

− I. PROOF. Let Si = T ∗

− I. It follows that (T ∗

T −1T ∗

− I) = (T ∗

− I)−1(T −1

− I) = (Si − T −1

)S−1

(Si − (T ∗

= Si − T −1

− (T ∗

S−1

(T ∗

= −I + ¯ T −1

. Q.E.D. 5.1. A special decomposition. In some important applications, such as multigrid methods, a special decomposition v =

a sequence of linear operators. If we assume that there exist operators: Πi : V → Vi such that R(Πi − Πi−1) ⊂ Vi, ∀v ∈ V with Π0 = 0 and ΠJ = I, we have a telescopic decomposition v =

vi with vi = (Πi − Πi−1)v.

SLIDE 13

Then, in view of the definition of c0, we have c0 ≤ sup

(Ti ¯ T −1

with (8) wi =

vj − T −1

vi = (I − Πi)v + (I − T −1

)(Πi − Πi−1)v. Furthermore c0 ≤ sup

(Ti ¯ T −1

≤ 2 sup

T −1

+(( ¯ T −1

− I)(Πi − Πi−1)v, (Πi − Πi−1)v)

An important special case is when Πi = Pi. In this case, the expression (8) is reduced to wi =

(I − T −1

)(Pi − Pi−1)v and (9) c0 ≤ sup

(( ¯ T −1

− I)(Pi − Pi−1)v, (Pi − Pi−1)v)

5.2. A multigrid method convergence proof. The estimate below estimate is very well known before and basic, but the proof here is different from those in the literature (Braess and Hackbusch and Bramble and Pasciak). We assume that Ω has been triangulated with a nested sequence of quasi-uniform triangulations Tk = {τ i

quasi-uniformity constants are independent of k. These triangulations should be nested in the sense that any triangle τ l

hk ∼ γ−k. Associated with each Tk, a finite element space Vk ⊂ H1

(10) V0 ⊂ V1 ⊂ . . . ⊂ Vk ⊂ . . . ⊂ VJ.

SLIDE 14

We are interested in multigrid method for solving the following finite element equation: Find uh ∈ V ≡ VJ satisfying (11) a(uh, v) = f(v) ∀v ∈ V where a(u, v) =

∇u · ∇v. Define Ai : Vi → Vi by: (Aiui, vi)0 = a(ui, vi), ∀ui, vi ∈ Vi. The subspace solver Ti in a multigrid method is often given by applying a number of, say m, smoothings using, for example, local relaxation method such as Gauss-Seidel iteration (with perhaps exception of i = 0 in which exact solver may be used, namely T0 = P0). Let Ri ≈ A−1

denote one such a smoothing. With m smoothings, one strategy is to apply Ri and Rt

Ri represents forward Gauss-Seidel iteration, then Rt

Ti =

if m is even Pi − Ki(K∗

if m is odd and ¯ T i = (I − Km

where Ki = I − RiAi, and K∗

and Ki,m =

if m is even KiK∗

if m is odd Note that K∗

In this example, we shall make two assumptions. The first assumption is on the approximation of the finite element subspaces: (12) (I − Pi−1)vi2 ≤ c1 λi a(vi, vi), ∀v ∈ Vi, where λi = ρ(Ai). The second assumption is on the smoother Ri: (13) c2 λi (v, v) ≤ ( ¯ Riv, v) ∀v ∈ Vi, where ¯ Ri is the corresponding symmetrization of Ri. We note that the first assumption is satisfied, for example, when aij are smooth and Ω is smooth or a convex Lipschitz domain and the second assumption is satisfied for Gauss-Seidel iteration. To estimate the constant c0, we consider the following decomposition v =

(14) vi = (Pi − Pi−1)v.

SLIDE 15

Since, thanks to (10), Pi−1 = Pi−1Pi = PiPi−1, (12) implies that (15) λi(vi, vi) ≤ c1a(vi, vi). We note that a( ¯ T −1

T i)vi, vi) = a((I − Km

= ( ¯ R−1

¯ RiAi(I − Km

≤ λi c2 ((I − Ki,m)(I − Km

≤ λi c2 max

≤ λi c2m(vi, vi). Hence, by (15), we have

a( ¯ T −1

T i)vi, vi) ≤

λi c2m(vi, vi) ≤

c1 c2ma(vi, vi) = c1 c2ma(v, v). Thus, by c0 ≤ sup

(( ¯ T −1

− I)(Pi − Pi−1)v, (Pi − Pi−1)v)

we have c0 ≤ c1 c2m. Consequently, the method of successive subspace corrections, based on multilevel subspaces (10) with a smoother satisfying (13), has the following convergence estimate: EJ2

c1 c1 + c2m. 5.3. More general elliptic PDEs.

Rd is Lipschitz domain), satisfying: a(u, v) = f(v), for all v ∈ H1

where a(u, v) =

∂v ∂xi αij ∂u ∂xj dx, f(v) =

fv dx, A = (αij) is symmetric and strictly positive definite for all x ∈ Ω, A and f are smooth enough so that the variational problem is well posed.

SLIDE 16

5.4. Discretization and subspace splitting.

a(uh, v) = f(v), for all v ∈ Vh, where Vh is a finite dimensional subspace.

Au = f, A S.P.D.

5.5. Basic multigrid iteration.

uk+1 = uk + B(f − Auk)

0. If k = 1 then B1g = A−1

in the 1D example 1. Pre-smoothing: x1 = ST

(I − Pnh) . . . (I − P1)(u − uk) 2. Coarse grid correction: (I − PH) 1. q0 = (P k

2. q1 = Bk−1q0; 3. x2 = x1 + P k

3. Post-smoothing: Bkg = x2 + Sk(g − Akx2). – 5.6. Basic facts from Finite element theory.

∀wk ∈ Vk

∀wk ∈ Vk

∈ Vk): ([Qτv], wk)τ = (v, wk)τ, ∀wk ∈ Vk|τ

SLIDE 17

5.7. Error estimates.

⋆ (I − Pk)va hkv2,Ω, ⋆ (I − Qk)v0 hkva ⋆ (I − Qτ)v0,τ hkv1,τ, ⋆ (I − Pk)v0 h2

⋆ (I − Pk)v0 hkv1,Ω, (if Aubin-Nitsche argument works). ⋆ (I − Qk)va va (Qk is energy stable) ⋆ λk = λmax(Ak) = h−2, rescaling (i.e. using ℓ2 instead of L2, gives λmin(Ak) = h2

5.8. On the addtive methods. nH subdomains ∪Ωk = Ω, H1

B := T = nH

Tk (B−1v, v)a ≈ (v, v)a. Theorem (Vassilevski, Wang?). Let and T = nH

(16) vk = ¯ Tk ¯ T −1v,

and (17) min

( ¯ T −1

wi, wi)a = (T −1v, v)a. A sketch of proof for ¯ Tk = Pk. inf

|vi|2

= inf

|wi + ξi|2

= inf

(|wi|2

(Awi, ξi) = inf

(|wi|2

(AiA−1

= inf

(|wi|2

ξi) = inf

(|wi|2

Finally we observe that

(Awi, wi) = (Aiwi, wi) = (AiA−1

wi) = (T −1v, v).

SLIDE 18

5.9. Multiplicative overlapping DD-method. (1) Partition of unity: J

suppθi ⊂ Ωi ∪ ∂Ω, 0 ≤ θi ≤ 1, maxx∈¯

and L2 projection Q0 : V → V0: h−1

Pk

vi2

(

θi)(v − Q0v))2

≤ Λ1

max

|

∇θi(x)|v − Q0v2

Λ1C0|v|2

Λ0 C0v2

Note: convergence rate does not depend on the possible osscilations in aij! 5.10. Convergence for finite element equations.

f ∈ V ∗

results in Gauss-Seidel method.

ei = M−1

Πi(f − Au), Mi ≈ Ai. 5.11. Convergence via subspace correction framework.

EJ = (I − TJ) . . . (I − T1), with Tk = M−1

Theorem.[XU AND Z, J. of AMS (2003)] The energy norm of the error transfer operator EJA satisfies EJ2

1 c0 + 1, where c0 = sup

inf

( ¯ T −1

vj − T −1

Tk ≡ T ∗

⇒ c0 < ∞.

SLIDE 19

K(v, ¯ Tk) = v2

( ¯ T −1

Tk).

K(v, ¯ Tk) =

( ¯ T −1

where for xk are xk =

T =

Tk. 5.12. Further simplifications: Richardson smoother.

K(v, Tk) 2 ≤ K(v, ¯ Tk) ≤ K(v, Tk)

K sup

inf

(T −1

5.13. Some known results.

K(v) ≤

λk(Pk − Pk−1)v2

K(v) ≤

λk(Qk − Qk−1)v2

(Pk − Qk)v2

SLIDE 20

5.14. Stability and convergence.

expansion: K(v) ≤

λk(Qk − Qk−1)v2

(Pk − Qk)v2

is true for nested FE spaces, and hence K(v) ≤

λk(Qk − Qk−1)v2

(Pk − Qk)v2

Such estimate without any regularity assumptions gives that the convergence rate depends only logarithmically on h (J ≈ | log2 h|). BASIC REFERENCES

MR 714691 (85h:65233)

(95b:65002)

Philadelphia, PA, 2000. MR 1774296 (2001h:65002)

(87e:65082)

Mathematics (SIAM), Philadelphia, PA, 2011.

Verlag, Berlin, 2005. MR 2104179 (2005g:65006)

uller, Multigrid, Academic Press Inc., San Diego, CA, 2001, With contributions by A. Brandt, P. Oswald and K. St¨

element equations. MR 2427040 (2010b:65007)