IDR as a Deflation Method Gerard Sleijpen Department of - - PowerPoint PPT Presentation

idr as a deflation method
SMART_READER_LITE
LIVE PREVIEW

IDR as a Deflation Method Gerard Sleijpen Department of - - PowerPoint PPT Presentation

Netherlands-Japan, NLA2012, Delft, April 10, 2012 IDR as a Deflation Method Gerard Sleijpen Department of Mathematics http://www.staff.science.uu.nl/ sleij101/ Gerard Sleijpen Joint work with Martin van Gijzen Delft University of


slide-1
SLIDE 1

Netherlands-Japan, NLA2012, Delft, April 10, 2012

IDR as a Deflation Method

Gerard Sleijpen Department of Mathematics

http://www.staff.science.uu.nl/∼sleij101/

slide-2
SLIDE 2

Gerard Sleijpen

Joint work with Martin van Gijzen Tijmen Collignon Delft University of Technology,

Delft Institue of Applied Mathematics

Delft, the Netherlands

slide-3
SLIDE 3

Program

  • Introduction
  • Sonneveld spaces
  • Inducing dimension reduction
  • IDR and deflation
slide-4
SLIDE 4

Program

  • Introduction
  • Sonneveld spaces
  • Inducing dimension reduction
  • IDR and deflation
slide-5
SLIDE 5

Solve Ax = b for x. A is n × n non-singular (complex).

slide-6
SLIDE 6

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b)

slide-7
SLIDE 7

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on angle b and AKk(A, b)

slide-8
SLIDE 8

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis

slide-9
SLIDE 9

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Extraction quality depends on specific method practical issues as efficiency, stability

slide-10
SLIDE 10

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Two ways to improve spectral properties:

  • Preconditioning. To cluster eigenvalues, . . .

Apply iterative method to K−1Ax = K−1b

  • Deflation. To replace small eigenvalues by 0

Remove known components

slide-11
SLIDE 11

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Two ways to improve spectral properties:

  • Preconditioning. To cluster eigenvalues, . . .

Apply iterative method to K−1Ax = K−1b

  • Deflation. To replace small eigenvalues by 0

Remove known components

slide-12
SLIDE 12

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Two ways to improve spectral properties:

  • Preconditioning. To cluster eigenvalues, . . .

Apply iterative method to K−1Ax = K−1b

  • Deflation. To replace small eigenvalues by 0

Remove known components

slide-13
SLIDE 13

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Two ways to improve spectral properties:

  • Preconditioning. To cluster eigenvalues, . . .

Apply iterative method to K−1Ax = K−1b

  • Deflation. To replace small eigenvalues by 0

Remove known components

slide-14
SLIDE 14

Krylov subspace methods

To solve Ax = b for x

  • Expansion.

Built a Krylov subspace Kk(A, b)

  • Extraction.

Extract an approximate solution xk from Kk(A, b) Convergence depends on spectral properties of A: – how well eigenvalues are clustered – on the “conditioning” of the eigenvector basis Two ways to improve spectral properties:

  • Preconditioning. To cluster eigenvalues, . . .

Apply iterative method to K−1Ax = K−1b

  • Deflation. To replace small eigenvalues by 0

Remove known components

slide-15
SLIDE 15

Program

  • Introduction
  • Sonneveld spaces
  • Inducing dimension reduction
  • IDR and deflation
slide-16
SLIDE 16

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix. Terminology. R is called the initial shadow residual or the IDR test matrix

slide-17
SLIDE 17

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix. Definition.

  • Block Krylov subspace of order k generated by A∗ and

R Kk(A∗, R) ≡

  

  • j<k

(A∗)j R γj | γj ∈ Cs

  

slide-18
SLIDE 18

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix. Definition.

  • Block Krylov subspace of order k generated by A∗ and

R Kk(A∗, R) ≡

  

  • j<k

(A∗)j R γj | γj ∈ Cs

  

Example. Bi-CG generates residuals rBi-CG

k

⊥ Kk(A∗, r0), here R = [ r0]: rBi-CG

k+1 = rBi-CG k

− Aukαk ⊥ Kk+1(A∗, r0) with uk such that Auk ⊥ Kk(A∗, r0).

slide-19
SLIDE 19

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

slide-20
SLIDE 20

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

  • Example.

R = [ r0 ], rBi-CGSTAB

k

= Pk(A)rBi-CG

k

Note that rBi-CG

k

⊥ Kk(A∗, r0). In Bi-CGSTAB: Pk+1(λ) = (1 − ωkλ)Pk(λ), where, with r′

k ≡ Pk(A)rBi-CG k+1 ,

ωk ≡ argminω∈Cr′

k − ωAr′ k2

Theorem. rBi-CGSTAB

k

∈ S(Pk, A, r0).

slide-21
SLIDE 21

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

  • [Sonneveld, van Gijzen, 2008, S. Sonneveld, van Gijzen 2010]
  • Property. Bi-CGSTAB ∼ IDR(s) for s = 1 (i.e. R = [

r0])

slide-22
SLIDE 22

100 200 300 400 500 600 −16 −14 −12 −10 −8 −6 −4 −2 2 4 # MV log10||rk||2 Convergence history linear solvers. Matrix "meier01" is 1000 x 1000 bicgstab BiCGstab(2) bi−idrs(8)

uxx + uyy + uzz + 1000 ux = f, f is defined by the solution u(x, y, z) = exp(xyz) sin(πx) sin(πy) sin(πz). Discretized with (10 × 10 × 10) volumes. No preconditioning.

slide-23
SLIDE 23

100 200 300 400 500 600 −16 −14 −12 −10 −8 −6 −4 −2 2 4 # MV log10||rk||2 Convergence history linear solvers. Matrix is 1000 x 1000 bicgstab BiCGstab(2) idr(2)stab(2) bi−idrs(8)

uxx + uyy + uzz + 1000 ux = f, f is defined by the solution u(x, y, z) = exp(xyz) sin(πx) sin(πy) sin(πz). Discretized with (10 × 10 × 10) volumes. No preconditioning.

slide-24
SLIDE 24

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

  • IDR Theorem. With

Pk+1(λ) ≡ (α − ωλ)Pk(λ), ω = 0

  • (αI − ωA)
  • S(Pk, A,

R) ∩ R⊥ = S(Pk+1, A, R)

  • S(Pk+1, A,

R) ⊂ S(Pk, A, R)

  • if A has no eigenvector in

R⊥ and S(Pk, A, R) = {0}.

slide-25
SLIDE 25

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

  • IDR Theorem. With

Pk+1(λ) ≡ (α − ωλ)Pk(λ), ω = 0

  • (αI − ωA)
  • S(Pk, A,

R) ∩ R⊥ = S(Pk+1, A, R)

  • S(Pk+1, A,

R) ⊂ S(Pk, A, R)

  • if A has no eigenvector in

R⊥ and S(Pk, A, R) = {0}. If zeros Pk+1 = eigenvalues A, then, increase k leads to dimension reduction S(Pk, A, R) = dimension increase Kk(A∗, R)

slide-26
SLIDE 26

s ∈ N, typically s = 1, 2, 4, 8. R is a full rank n × s matrix.

  • Definition. Pk polynomial of exact degree k.
  • Pk Sonneveld subspace order k generated by A and

R S(Pk, A, R) ≡

  • Pk(A) v | v ⊥ Kk(A∗,

R)

  • IDR Theorem. With

Pk+1(λ) ≡ (α − ωλ)Pk(λ), ω = 0

  • (αI − ωA)
  • S(Pk, A,

R) ∩ R⊥ = S(Pk+1, A, R)

  • S(Pk+1, A,

R) ⊂ S(Pk, A, R)

  • if A has no eigenvector in

R⊥ and S(Pk, A, R) = {0}. Corollary Gk ≡ S(Pk, A, R), G′

k ≡ Gk ∩ R⊥.

With Pk+1(λ) ≡ (1 − ωkλ)Pk(λ), we have Gk+1 ≡ (I − ωkA)G′

k ⊂ Gk

and AG′

k ⊂ Gk

slide-27
SLIDE 27

Program

  • Introduction
  • Sonneveld spaces
  • Inducing dimension reduction
  • IDR and deflation
slide-28
SLIDE 28

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk in Gk iteratively by increase k

slide-29
SLIDE 29

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk in Gk The first residual to be in Gk by construction is called a primary residual. The implementation may rely on other residuals: secondary residuals

slide-30
SLIDE 30

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk in Gk All our updates for residuals are of the form r ← r − αc with c = Au and u available. Hence, x ← x + αu. Updates r and x are consistent. Approximate solution x gets a(n almost) free ride:

  • nly vector updates, no MVs, no inner products.

Focuss on r

slide-31
SLIDE 31

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: rk ∈ Gk, then construct residual rk+1 in Gk+1 Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 is a skew projection that projects onto R⊥

slide-32
SLIDE 32

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: rk ∈ Gk, then construct residual rk+1 in Gk+1 Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 is a skew projection that projects onto R⊥

slide-33
SLIDE 33

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: rk ∈ Gk, then construct residual rk+1 in Gk+1 Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 is a skew projection that projects onto R⊥ Π1 = I − Vσ−1 R∗ with σ ≡ R∗V s × s non-singular, and V ≡ V

k

is n × s matrix with span(V) ⊂ Gk [rk, V

k] full rank

slide-34
SLIDE 34

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: rk ∈ Gk, then construct residual rk+1 in Gk+1 Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 is a skew projection that projects onto R⊥ Π1 = I − Vσ−1 R∗ with σ ≡ R∗V s × s non-singular, and V ≡ V

k

is n × s matrix with span(V) ⊂ Gk [rk, V

k] full rank

Π1(vj) = 0 for the columns of V = [v1, . . . , vs]

slide-35
SLIDE 35

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k =????

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

Π1 is a skew projection that projects onto R⊥

slide-36
SLIDE 36

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k =????

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

Π1 is a skew projection that projects onto R⊥ w ∈ G′

k

⇒ Aw ∈ Gk ⇒ Π1Aw ∈ G′

k

slide-37
SLIDE 37

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k =????

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

Π1 is a skew projection that projects onto R⊥ w ∈ G′

k

⇒ Aw ∈ Gk ⇒ Π1Aw ∈ G′

k

⇒ Ks+1(Π1A, r′

k) ⊂ G′ k

slide-38
SLIDE 38

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k =????

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

Π1 is a skew projection that projects onto R⊥ w ∈ G′

k

⇒ Aw ∈ Gk ⇒ Π1Aw ∈ G′

k

⇒ Ks+1(Π1A, r′

k) ⊂ G′ k

  • V′

k with span(r′ k, V′ k) = Ks+1(Π1A, r′ k)

slide-39
SLIDE 39

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k =????

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

Π1 is a skew projection that projects onto R⊥ w ∈ G′

k

⇒ Aw ∈ Gk ⇒ Π1Aw ∈ G′

k

⇒ Ks+1(Π1A, r′

k) ⊂ G′ k

  • V′

k with span(r′ k, V′ k) = Ks+1(Π1A, r′ k);

AV′

k is side product

slide-40
SLIDE 40

Induced Dimension Reduction

Select (ωk),

  • R. Let Pk+1(λ) = (1 − ωkλ)Pk(λ).

G0 = Cn, Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩

R⊥. IDR: construct residuals rk+1 in Gk+1, V

k+1 span in Gk+1.

Gk V

k

IDR step → Π1 G′

k ≡ Gk ∩

R V′

k = !!

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

V

k+1 = V′ k − ωkAV′ k

  • Theorem. Assume

span([rk, V

k]) ⊂ Gk.

If [r′

k, V′ k] spans Ks+1(Π1A, r′ k) and no break-down occurs,

then span([rk+1, V

k+1]) ⊂ Gk+1.

This is essentially the only way to move to the next Gk

slide-41
SLIDE 41

[r′

k, V′ k] spans Ks+1(Π1A, r′ k). How to select the basis V k?

(In exact arithmetic) any basis leads to the same projection.

slide-42
SLIDE 42

Program

  • Introduction
  • Sonneveld spaces
  • Inducing dimension reduction
  • IDR and deflation
slide-43
SLIDE 43

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗

slide-44
SLIDE 44

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Deflated system: Solve Π1Ax′ = r′

k

for x′ ⊥ R. (∗) Then x = xk + Qrk + (I − QA)x′ solves Ax = b. Note. Π1A : R⊥ → R⊥ and Π1rk ∈ R⊥. Kk(Π1A, Π1rk) leads to approximate solutions of (∗) in R⊥.

slide-45
SLIDE 45

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Deflated system: Solve Π1Ax′ = r′

k

for x′ ⊥ R. (∗) Then x = xk + Qrk + (I − QA)x′ solves Ax = b. Note. Π1A : R⊥ → R⊥ and Π1rk ∈ R⊥. Kk(Π1A, Π1rk) leads to approximate solutions of (∗) in R⊥. Solve (∗) with s steps of some Krylov method:

  • Advantage. At the same time

a basis V′ of Ks+1(Π1A, Π1rk) and a smaller residual.

slide-46
SLIDE 46

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Deflated system: Solve Π1Ax′ = r′

k

for x′ ⊥ R. (∗) Then x = xk + Qrk + (I − QA)x′ solves Ax = b. Note. Π1A : R⊥ → R⊥ and Π1rk ∈ R⊥. Kk(Π1A, Π1rk) leads to approximate solutions of (∗) in R⊥. However, if ss is an s-step Krylov residual, then Π1(ss − ωkAss) = rk+1.

slide-47
SLIDE 47

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗

  • Interpretation. Solve Au = r

Search for an approximate solution in span(U) R∗AUα = R∗r Then u ≈ Uα = U(R∗AU)−1R∗r = Qr

slide-48
SLIDE 48

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗

  • Interpretation. Solve Au = r

Search for an approximate solution in span(U) R∗AUα = R∗r Then u ≈ Uα = U(R∗AU)−1R∗r = Qr In multigrid: R∗ represents the restriction operator U represents the prolongation Q is the coarse grid correction

slide-49
SLIDE 49

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗

  • Interpretation. Solve Au = r

Search for an approximate solution in span(U) R∗AUα = R∗r Then u ≈ Uα = U(R∗AU)−1R∗r = Qr Then the update of the residual r is r′

k = rk − AQrk

with update of the approximate solution x′

k = xk + Qrk

slide-50
SLIDE 50

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Interpretation. Coarse grid correction: r′

k = rk − AQrk and x′ k = xk + Qrk

Smoothing: rk+1 = r′

k − ωAr′ k and xk+1 = x′ k + ωr′ k

slide-51
SLIDE 51

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Interpretation. Coarse grid correction: r′

k = rk − AQrk and x′ k = xk + Qrk

Smoothing: rk+1 = r′

k − ωAr′ k and xk+1 = x′ k + ωr′ k

Combine: rk+1 = (I − ωA)Π1rk = rk − A(Q + ωΠ1)r xk+1 = xk + (Q + ωΠ1)rk With P ≡ Q + ωΠ, I − AP is the residual reduction operator

slide-52
SLIDE 52

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Interpretation. Coarse grid correction: r′

k = rk − AQrk and x′ k = xk + Qrk

Smoothing: rk+1 = r′

k − ωAr′ k and xk+1 = x′ k + ωr′ k

Combine: rk+1 = (I − ωA)Π1rk = rk − A(Q + ωΠ1)r xk+1 = xk + (Q + ωΠ1)rk With P ≡ Q + ωΠ, I − PA is the error reduction operator

slide-53
SLIDE 53

Gk rk IDR step → Π1 G′

k ≡ Gk ∩

R⊥ r′

k = Π1rk

  • pol. step

→ I − ωkA Gk+1 ≡ (I − ωkA)G′

k

rk+1 = r′

k − ωkAr′ k

Π1 = I − AUσ−1 R∗ with σ ≡ R∗AU and V = AU Π1 = I − AQ with Q ≡ Uσ−1 R∗ Interpretation. Coarse grid correction: r′

k = rk − AQrk and x′ k = xk + Qrk

Smoothing: rk+1 = r′

k − ωAr′ k and xk+1 = x′ k + ωr′ k

Combine: rk+1 = (I − ωA)Π1rk = rk − A(Q + ωΠ1)r xk+1 = xk + (Q + ωΠ1)rk With P ≡ Q + ωΠ, I − PA is the error reduction operator IDR can be viewed as Richardson with a flexible preconditioner

slide-54
SLIDE 54

Spectrum of IDR’s error reduction operator

Λ(I − PA) = Λ(I − AP) = Λ(Π1(I − ωA)) The eigenvalues of the IDR error reduction operator I − PA are related to the eigenvalues of the IDR deflated matrix Π1A.

[Erlangga, Nabben, 2008]

  • Theorem. For a λ ∈ C, λ = 0 we have that

if Π1Av = λv then Π1(I − ωA)v = (1 − ωλ)v If Π1Av = 0 then (I − PA)v = 0.

slide-55
SLIDE 55

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk as in IDR

slide-56
SLIDE 56

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk Theorem

  • 0 is an eigenvalue of Π1A with geometric multiplicity ≥ s.
  • If Pk(µ) = 0 (i.e., µ = 1/ωj for some j < k), then

µ is an eigenvalue of Π1A with geometric multiplicity ≥ s

  • If Pk(µ) = Pk(µ) = . . . = P (ℓ−1)

k

(µ) = 0 then µ is an eigenvalue of Π1A with algebraic multiplicity ≥ ℓs If span(V) ⊂ Gk, then V = Pk(A)V for some n × s matrix V ⊥ Kk(A∗, R).

  • If

R∗(A − λI)−1V is singular, then λ is an eigenvalue of Π1A.

slide-57
SLIDE 57

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk If span(V) ⊂ Gk, then V = Pk(A)V for some n × s matrix V ⊥ Kk(A∗, R).

  • Theorem. V is independent of Pk in IDR.
slide-58
SLIDE 58

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk

  • Comment. V

k only spans an s-dimensional subspace.

(with expected deflation as indicated by the s-fold eig. 0)

slide-59
SLIDE 59

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk

  • Comment. V

k only spans an s-dimensional subspace.

(with expected deflation as indicated by the s-fold eig. 0) Though V

k has a “memory”

  • f all ks vectors in the preceeding V

j

(j < k): it allows control (clustering) of ks eigenvalues.

slide-60
SLIDE 60

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk Theorem

  • If Pk(µ) = 0 (i.e., µ = 1/ωj for some j < k), then

1 − ωkµ is an eigenvalue of (I − ωkA)Π1 with multiplicity s, 0 is an eigenvalue of (I − ωkA)Π1 with multiplicity s.

slide-61
SLIDE 61

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk Theorem

  • If Pk(µ) = 0 (i.e., µ = 1/ωj for some j < k), then

1 − ωkµ is an eigenvalue of (I − ωkA)Π1 with multiplicity s, 0 is an eigenvalue of (I − ωkA)Π1 with multiplicity s.

  • Comment. When selecting ωj = 1 (j ≤ k),

0 is an eigenvalue of (I−ωkA)Π1 with multiplicity (k +1)s.

slide-62
SLIDE 62

The spectrum of IDR’s deflated operator

With (ωk), Pk+1(λ) = (1 − ωkλ)Pk(λ), Gk+1 ≡ (I − ωkA)G′

k

with G′

k ≡ Gk ∩ R⊥,

Let V be such that span(V) ⊂ Gk Theorem

  • If Pk(µ) = 0 (i.e., µ = 1/ωj for some j < k), then

1 − ωkµ is an eigenvalue of (I − ωkA)Π1 with multiplicity s, 0 is an eigenvalue of (I − ωkA)Π1 with multiplicity s.

  • If V

k is constructed with the IDR scheme, then

the remaining n − (k + 1)s eigenvalues do not depend on ωj, j < k.

slide-63
SLIDE 63

Conclusions

  • IDR is a [class of] very effective methods
  • IDR can be viewed as a deflation method with a

flexible preconditioner

  • The deflated matrices as produced in IDR have

remarkable spectral properties

  • It is not clear how to exploit these elegant relations

and properties to explain the effectiveness of IDR

slide-64
SLIDE 64

References

  • 1. Sleijpen, Sonneveld, and van Gijzen,

“Bi-CGSTAB as an induced dimension reduction method”, in APNUM 60, 1100–1114, 2010

  • 2. Sleijpen and van Gijzen,

“BiCGstab(ℓ) strategies to induce dimension reduction”, SISC, 32 (5), 2687–2709, 2010

  • 3. Sleijpen and van Gijzen,

“Subspaces for inducing dimension reduction”, Proceeding of the Krylov Forum, Kyoto, March, 2010.

  • 4. Collignon, Sleijpen, and van Gijzen,

“Interpreting IDR(s) as a deflation method”, Report, Delft University, October 2010.