Rational Krylov methods for f ( A ) b Michael Eiermann, Oliver G. - - PowerPoint PPT Presentation

▶
rational krylov methods for f a b
SMART_READER_LITE
LIVE PREVIEW

Rational Krylov methods for f ( A ) b Michael Eiermann, Oliver G. - - PowerPoint PPT Presentation

Dec 27, 2023 298 likes •731 views

Rational Krylov methods for f ( A ) b Michael Eiermann, Oliver G. Ernst, and Stefan Gttel DWCAA09 September 8th, 2009 Pr sr t tr f ( A ) b r A s r

slide-1
SLIDE 1

Rational Krylov methods for f(A)b

Michael Eiermann, Oliver G. Ernst, and Stefan Güttel DWCAA09 September 8th, 2009

slide-2
SLIDE 2

Pr♦❜❧❡♠

❲❡ ❝♦♥s✐❞❡r t❤❡ ✈❡❝t♦r f(A)b✱ ✇❤❡r❡

◮ A ✐s ❛ ❧❛r❣❡ N✲❜②✲N ♠❛tr✐①✱ ◮ b ✐s ❛ ✈❡❝t♦r ♦❢ ❧❡♥❣t❤ N✱ ◮ f ✐s ❛ s✉✐t❛❜❧❡ ❢✉♥❝t✐♦♥✳

❈♦♠♣✉t❡ ❛♣♣r♦①✐♠❛t✐♦♥ fm ≈ f(A)b ❢r♦♠ ❛ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡✳

slide-3
SLIDE 3

❲❤❛t ✐s ❛ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡❄

▲❡t {ξ1, ξ2, . . .} ⊆  ⊂ C ❜❡ ❛ ❣✐✈❡♥ s❡q✉❡♥❝❡ ♦❢ ♣♦❧❡s✳ ❉❡✜♥❡ t❤❡ ♣♦❧②♥♦♠✐❛❧s qm(z) =

m

  • j=1

ξj=∞

(z − ξj) ∈ Pm. ❆ss✉♠❡ t❤❛t qm(A)−1 ❡①✐sts✳ ❚❤❡♥ Qm+1(A, b) = Km+1(A, qm(A)−1b) ✐s t❤❡ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡ ❛ss♦❝✐❛t❡❞ ✇✐t❤ (A, b, qm)✳

slide-4
SLIDE 4

❙♣❡❝✐❛❧ ❝❛s❡s

◮  = {∞} ⇒ ♣♦❧②♥♦♠✐❛❧ ❑r②❧♦✈ Qm+1 = Km+1

❬◆❛✉ts ✫ ❲②❛tt ✽✸❪ ❬✈❛♥ ❞❡r ❱♦rst ✽✼❪ ❬❉r✉s❦✐♥ ✫ ❑♥✐③❤♥❡r♠❛♥ ✽✽❪ ❬●❛❧❧♦♣♦✉❧♦s ✫ ❙❛❛❞ ✾✷❪ ❬❍♦❝❤❜r✉❝❦ ✫ ▲✉❜✐❝❤ ✾✼❪ ❬❊✐❡r♠❛♥♥ ✫ ❊r♥st ✵✻❪

◮  = {ξ} ⇒ s❤✐❢t✲✐♥✈❡rt ❑r②❧♦✈

❬▼♦r❡t ✫ ◆♦✈❛t✐ ✵✹❪ ❬✈❛♥ ❞❡♥ ❊s❤♦❢ ✫ ❍♦❝❤❜r✉❝❦ ✵✻❪

◮  = {0, ∞} ⇒ ❡①t❡♥❞❡❞ ❑r②❧♦✈

❬❉r✉s❦✐♥ ✫ ❑♥✐③❤♥❡r♠❛♥ ✾✽❪ ❬❑♥✐③❤♥❡r♠❛♥ ✫ ❙✐♠♦♥❝✐♥✐ ✵✽❪

◮  ❛r❜✐tr❛r② ⇒ r❛t✐♦♥❛❧ ❑r②❧♦✈

❬❘✉❤❡ ✽✹❪ ❬❇❡❛tt✐❡ ✵✹❪ ❬❇❡❝❦❡r♠❛♥♥ ✫ ❘❡✐❝❤❡❧ ✵✽❪ ❬❑♥✐③❤♥❡r♠❛♥ ❡t ❛❧ ✵✽❪

slide-5
SLIDE 5

❘❛t✐♦♥❛❧ ❆r♥♦❧❞✐ ❛❧❣♦r✐t❤♠ ❬❘✉❤❡ ✽✹✴✾✹❪

Input✿ A, b, {ξ1, ξ2, . . . , ξm} v1 := b/b for j = 1, 2, . . . , m do x := (I − A/ξj)−1Avj H(1 : j, j) := [v1, . . . , vj]∗x x := x − [v1, . . . , vj]H(1 : j, j) H(j + 1, j) := x vj+1 := x/H(j + 1, j) end ❨✐❡❧❞s ❞❡❝♦♠♣♦s✐t✐♦♥ AVm+1(Im + HmX−1

m ) = Vm+1Hm✳

slide-6
SLIDE 6

❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❞❡❝♦♠♣♦s✐t✐♦♥s

❚❤❡♦r❡♠ ✭●✳✱ ✷✵✵✾✮✿ ▲❡t ❛ ❣❡♥❡r❛❧ ❞❡❝♦♠♣♦s✐t✐♦♥ AVm+1Km = Vm+1Hm ❜❡ ❣✐✈❡♥✱ ✇❤❡r❡ Vm+1 ❤❛s m + 1 ❧✐♥❡❛r❧② ✐♥❞❡♣❡♥❞❡♥t ❝♦❧✉♠♥s✱ Km ∈ C(m+1)×m✱ Hm ∈ C(m+1)×m✱ ❛♥❞ Hm ✐s ♦❢ r❛♥❦ m✳ ❚❤❡♥ ✐s ♦❢ r❛♥❦ ✳ ❢♦r ❛ ✈❡❝t♦r ✳ ❋♦r ❡✈❡r② ✈❡❝t♦r t❤❡r❡ ❡①✐sts ❛ ✉♥✐q✉❡ ♣♦❧②♥♦♠✐❛❧ ✱ ✱ s✉❝❤ t❤❛t ✳ ❍❡♥❝❡✱ ✐❢ ✐s ✐♥✈❡rt✐❜❧❡✱ ✳

slide-7
SLIDE 7

❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❞❡❝♦♠♣♦s✐t✐♦♥s

❚❤❡♦r❡♠ ✭●✳✱ ✷✵✵✾✮✿ ▲❡t ❛ ❣❡♥❡r❛❧ ❞❡❝♦♠♣♦s✐t✐♦♥ AVm+1Km = Vm+1Hm ❜❡ ❣✐✈❡♥✱ ✇❤❡r❡ Vm+1 ❤❛s m + 1 ❧✐♥❡❛r❧② ✐♥❞❡♣❡♥❞❡♥t ❝♦❧✉♠♥s✱ Km ∈ C(m+1)×m✱ Hm ∈ C(m+1)×m✱ ❛♥❞ Hm ✐s ♦❢ r❛♥❦ m✳ ❚❤❡♥

  • 1. Km ✐s ♦❢ r❛♥❦ m✳
  • 2. colspan(Vm+1) = Km+1(A, q) ❢♦r ❛ ✈❡❝t♦r q✳
  • 3. ❋♦r ❡✈❡r② ✈❡❝t♦r b ∈ colspan(Vm+1) t❤❡r❡ ❡①✐sts ❛ ✉♥✐q✉❡

♣♦❧②♥♦♠✐❛❧ qm✱ deg(qm) ≤ m✱ s✉❝❤ t❤❛t b = qm(A)q✳ ❍❡♥❝❡✱ ✐❢ qm(A) ✐s ✐♥✈❡rt✐❜❧❡✱ colspan(Vm+1) = Qm+1(A, b)✳

slide-8
SLIDE 8

❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❛♣♣r♦①✐♠❛t✐♦♥s

❆ s♣❡❝✐❛❧ ❝❛s❡ ✐s t❤❡ ✭r❡❞✉❝❡❞✮ ❞❡❝♦♠♣♦s✐t✐♦♥ AVmKm = Vm+1Hm. ❆s ❛♥ ❛♣♣r♦①✐♠❛t✐♦♥ t♦ f(A)b ✇❡ ❝♦♥s✐❞❡r fm := Vmf(HmK−1

m )V† mb.

❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮✿ ❚❤❡r❡ ❤♦❧❞s fm = rm(A)b = pm−1 qm−1 (A)b, ✇❤❡r❡ rm ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ(HmK−1

m )✳

slide-9
SLIDE 9

❊①❛♠♣❧❡

❚❤❡ ✐t❡r❛t✐♦♥ v1 = b, βjvj+1 = (I − A/ξj)−1(A − αjI)vj, j = 1, . . . , m, ②✐❡❧❞s ❛ ❞❡❝♦♠♣♦s✐t✐♦♥ AVm+1Km = Vm+1Hm ✇✐t❤ Vm+1 = [v1, . . . , vm+1]✱

Km =             1

β1/ξ1

1

β2/ξ2

✳✳✳ ✳✳✳ 1

βm/ξm

            ❛♥❞ Hm =             α1 β1 α2 β2 ✳✳✳ ✳✳✳ αm βm             .

slide-10
SLIDE 10

❊①❛♠♣❧❡

❚❤❡ ✐t❡r❛t✐♦♥ v1 = b, βjvj+1 = (I − A/ξj)−1(A − αjI)vj, j = 1, . . . , m, ❝❛♥ ❜❡ ✉s❡❞ ❢♦r ❡①♣❧✐❝✐t r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥✿ ❇② ❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮ ✇❡ ❦♥♦✇ t❤❛t fm = Vmf(HmK−1

m )e1 = rm(A)b =

pm−1 qm−1 (A)b, ✇❤❡r❡ rm ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ(HmK−1

m ) = {α1, . . . , αm}✳

P A I N m e t h

  • d

Preassigned Poles and Interpolation Nodes

slide-11
SLIDE 11

❊①❛♠♣❧❡

❚❤❡ ✐t❡r❛t✐♦♥ v1 = b, βjvj+1 = (I − A/ξj)−1(A − αjI)vj, j = 1, . . . , m, ❝❛♥ ❜❡ ✉s❡❞ ❢♦r ❡①♣❧✐❝✐t r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥✿ ❇② ❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮ ✇❡ ❦♥♦✇ t❤❛t fm = Vmf(HmK−1

m )e1 = rm(A)b =

pm−1 qm−1 (A)b, ✇❤❡r❡ rm ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ(HmK−1

m ) = {α1, . . . , αm}✳

P A I N m e t h

  • d

Preassigned Poles and Interpolation Nodes

slide-12
SLIDE 12

❘❡♠❛r❦s

◮ 2 ✈❡❝t♦rs st♦r❛❣❡ ♥❡❡❞✱ 0 ✐♥♥❡r✲♣r♦❞✉❝ts ◮ ■❢ ❛❧❧ ξj = ∞ ⇒ ♣♦❧②♥♦♠✐❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ❛t {α1, . . . , αm} ◮ P♦❧②♥♦♠✐❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ♠❡t❤♦❞s ❤❛✈❡ ❜❡❡♥ ❝♦♥s✐❞❡r❡❞

❜❡❢♦r❡ ❬❍✉✐s✐♥❣❛ ❡t ❛❧ ✾✾❪ ❬❇❡r❣❛♠❛s❝❤✐✱ ❈❛❧✐❛r✐ ✫ ❱✐❛♥❡❧❧♦ ✵✹❪

◮ ❋♦r {α1, . . . , αm} ✉s❡ ▲❡❥❛ ♣♦✐♥ts✱ s❝❛❧❡❞ t♦ ❛ s❡t ♦❢ ✉♥✐t

❝❛♣❛❝✐t② ❢♦r st❛❜✐❧✐t② ❬❘❡✐❝❤❡❧ ✾✵❪✳

◮ ◆♦ s✉❝❤ s❝❛❧✐♥❣ ✐s ♥❡❝❡ss❛r② ✇✐t❤ t❤❡ P❆■◆ ♠❡t❤♦❞✿

s✐♠♣❧② ❝❤♦♦s❡ βj s✉❝❤ t❤❛t vj+1 = 1✱ j = 1, . . . , m✳

◮ ❋♦r r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ✉s❡ ▲❡❥❛✲❇❛❣❜② ♣♦✐♥ts✳

slide-13
SLIDE 13

❈♦♠♣✉t❡✿ f(A)b =

  • Ab✱ A = diag(1, . . . , 1000)✱ b = [1, . . . , 1]T✳
slide-14
SLIDE 14

❈♦♠♣✉t❡✿ f(A)b =

  • Ab✱ A = diag(1, . . . , 1000)✱ b = [1, . . . , 1]T✳
slide-15
SLIDE 15

20 40 60 80 100 120 140 160 10

−15

10

−10

10

−5

10

  • rder m

error

R ≈ 0.36 R =

κ−1 κ+1 ≈ 0.94

slide-16
SLIDE 16

20 40 60 80 100 120 140 160 10

−15

10

−10

10

−5

10

  • rder m

error

❩♦❧♦t❛r❡✈✬s ♣♦❧❡s

slide-17
SLIDE 17

■❢ ❣♦♦❞ ✭♦r ❜❡st✮ r❛t✐♦♥❛❧ ❛♣♣r♦①✐♠❛t✐♦♥ r∗

m t♦ f ✐s ❦♥♦✇♥

❡①♣❧✐❝✐t❧②✱ ♦♥❡ ❝❛♥ ❞✐r❡❝t❧② ❡✈❛❧✉❛t❡ r∗

m(A)b ≈ f(A)b✳ ❬❚r❡❢❡t❤❡♥ ❡t ❛❧ ✵✻❪ ❬❋r♦♠♠❡r ❡t ❛❧ ✵✻❪ ❬❙❝❤♠❡❧③❡r ❡t ❛❧ ✵✼❪ ❬❍❛❧❡ ❡t ❛❧ ✵✽❪

❍♦✇❡✈❡r✱ ✉s✐♥❣ t❤❡ ♣♦❧❡s ξj ♦❢ r∗

m ❛♥❞ s✉✐t❛❜❧❡ ✐♥t❡r♣♦❧❛t✐♦♥

♥♦❞❡s αj ❛s ✐♥♣✉ts ❢♦r P❆■◆✱ ✇❡ ❝❛♥ ❛❝❤✐❡✈❡ ❡ss❡♥t✐❛❧❧② t❤❡ s❛♠❡ ❛❝❝✉r❛❝② ❛t t❤❡ s❛♠❡ ❝♦♠♣✉t❛t✐♦♥❛❧ ❝♦st✳ ▼♦r❡♦✈❡r✱ t❤❡ P❆■◆ ♠❡t❤♦❞ ✐s ✐♠♣❧✐❝✐t❧② ❜❛s❡❞ ♦♥ ❡①❛❝t ✐♥t❡r♣♦❧❛t✐♦♥ ♦❢ f ❛♥❞ ❤❡♥❝❡ r♦❜✉st t♦ ♣❡rt✉r❜❛t✐♦♥s ✐♥ r∗

m✿

limsupm→∞ f(A)b − fm 1/m ≤ R < 1 ✐❢ αj✱ ξj ❛r❡ ❡q✉✐❧✐❜r✐✉♠✲❞✐str✐❜✉t❡❞ ♦♥ ✱ ✳

slide-18
SLIDE 18

❘❛②❧❡✐❣❤✲❘✐t③ ❡①tr❛❝t✐♦♥

❚❤❡r❡ ✐s ❛ ✇❛② t♦ ❛✉t♦♠❛t✐❝❛❧❧② ❝❤♦♦s❡ ♥❡❛r✲♦♣t✐♠❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ♣♦✐♥ts {α1, . . . , αm} ❛t ✐t❡r❛t✐♦♥ m✿

  • 1. ❈♦♠♣✉t❡ ♦rt❤♦♥♦r♠❛❧ ❜❛s✐s Vm ♦❢ Qm = q−1

m−1Km✳

  • 2. ✏❉❡t❡r♠✐♥❡✑ ❘❛②❧❡✐❣❤ q✉♦t✐❡♥t Am = V∗

mAVm✳

  • 3. ❈♦♠♣✉t❡ fm = Vmf(Am)V∗

mb✳

❚❤❡♦r❡♠✿ {α1, . . . , αm} = Λ(Am)✳ ❚❤❡♦r❡♠✿ f(A)b − fm ≤ Cminp∈Pm−1 f − p/qm−1F(A)✳ Pr✐❝❡✿ m ✈❡❝t♦rs st♦r❛❣❡ ♥❡❡❞✱ m2/2 ✐♥♥❡r✲♣r♦❞✉❝ts✳

slide-19
SLIDE 19

❈♦♠♣✉t❡✿ f(A)b = log(A)b✱ A ♥♦r♠❛❧ ✇✐t❤ ✶✵✵✵ ❡✐❣❡♥✈❛❧✉❡s

slide-20
SLIDE 20

❈♦♠♣✉t❡✿ f(A)b = log(A)b✱ A ♥♦r♠❛❧ ✇✐t❤ ✶✵✵✵ ❡✐❣❡♥✈❛❧✉❡s

slide-21
SLIDE 21

5 10 15 20 25 10

−15

10

−10

10

−5

10

  • rder m

error Leja−Bagby Rayleigh−Ritz

slide-22
SLIDE 22

5 10 15 20 25 10

−15

10

−10

10

−5

10

  • rder m

error Leja−Bagby Rayleigh−Ritz

♣♦❧②♥♦♠✐❛❧ ♠❡t❤♦❞

slide-23
SLIDE 23

P❛r❛♠❡t❡r✲❞❡♣❡♥❞❡♥t ♣r♦❜❧❡♠s

■♥ ♣r❛❝t✐❝❡✱ ♦♥❡ ♦❢t❡♥ ✐s ♥♦t ✐♥t❡r❡st❡❞ ✐♥ f(A)b ❜✉t ✐♥ f τ(A)b✱ ❢♦r τ ∈ T ❢r♦♠ s♦♠❡ ♣❛r❛♠❡t❡r s❡t T✳

  • ✐✈❡♥ ❛ s✐♥❣❧❡ r❛t✐♦♥❛❧ ❑r②❧♦✈ ❞❡❝♦♠♣♦s✐t✐♦♥ ✭❛s ❜❡❢♦r❡✮

AVmKm = Vm+1Hm, R(Vm) = Qm = q−1

m−1Km,

✇❡ ❝♦♠♣✉t❡ s❡✈❡r❛❧ ❛♣♣r♦①✐♠❛t✐♦♥s fτ

m = Vmf τ(HmK−1 m )V† mb = rτ m(A)b =

pτ

m−1

qm−1 (A)b, ✇❤❡r❡ rτ

m ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f τ ❛t Λ(HmK−1 m )✳

slide-24
SLIDE 24

❊①❛♠♣❧❡✿ ❚r❛♥s❢❡r ❢✉♥❝t✐♦♥

f τ(z) = (z − τ)−1✱ s♣❡❝tr✉♠  = [0, +∞)✱ ♣❛r❛♠❡t❡rs T = i[1, c]✳ ▲❡t ωm(z) = (z − α1) · · · (z − αm)✱ t❤❡♥ rτ

m(z) =

1 − qm−1(τ)

qm−1(z) ωm(z) ωm(τ)

z − τ = pτ

m−1

qm−1 ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f τ ❛t {α1, . . . , αm}✳ ❍❡♥❝❡✱ fτ

m = rτ m(A)b✳

❚❤❡ r❡❧❛t✐✈❡ ❡rr♦r ✐s ❛♥❞ ✐❢ ✱ ✐ts ♠✐♥✐♠✐③❛t✐♦♥ ✐s r❡❧❛t❡❞ t♦ t❤❡ ❆❉■ ♣r♦❜❧❡♠✳

❬❑♥✐③❤♥❡r♠❛♥✱ ❉r✉s❦✐♥ ✫ ❩❛s❧❛✈s❦② ✵✽❪

slide-25
SLIDE 25

❊①❛♠♣❧❡✿ ❚r❛♥s❢❡r ❢✉♥❝t✐♦♥

f τ(z) = (z − τ)−1✱ s♣❡❝tr✉♠  = [0, +∞)✱ ♣❛r❛♠❡t❡rs T = i[1, c]✳ ▲❡t ωm(z) = (z − α1) · · · (z − αm)✱ t❤❡♥ rτ

m(z) =

1 − qm−1(τ)

qm−1(z) ωm(z) ωm(τ)

z − τ = pτ

m−1

qm−1 ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f τ ❛t {α1, . . . , αm}✳ ❍❡♥❝❡✱ fτ

m = rτ m(A)b✳

❚❤❡ r❡❧❛t✐✈❡ ❡rr♦r ✐s [f τ(z) − rτ

m(z)]/f τ(z) =

qm−1(τ) qm−1(z) ωm(z) ωm(τ) , z ∈ , τ ∈ T, ❛♥❞ ✐❢  = T✱ ✐ts ♠✐♥✐♠✐③❛t✐♦♥ ✐s r❡❧❛t❡❞ t♦ t❤❡ ❆❉■ ♣r♦❜❧❡♠✳

❬❑♥✐③❤♥❡r♠❛♥✱ ❉r✉s❦✐♥ ✫ ❩❛s❧❛✈s❦② ✵✽❪

slide-26
SLIDE 26

−20 −10 10 20 30 40 50 5 10

slide-27
SLIDE 27

−20 −10 10 20 30 40 50 5 10 10 10

1

10

2

10

3

0.5 1 parameter spread c

convrate R

slide-28
SLIDE 28

■♥ ♣r❛❝t✐❝❛❧ ❝♦♠♣✉t❛t✐♦♥s ✐t ✐s ♠♦r❡ ❝♦♥✈❡♥✐❡♥t t♦ ❤❛✈❡ r❡❛❧ ♣♦❧❡s✱ ✐✳❡✳✱  = [−∞, 0)✳ ❍♦✇ t♦ s❡❧❡❝t {ξ1, ξ2, . . .} ⊂ ❄ −20 −10 10 20 30 40 50 5 10

slide-29
SLIDE 29

■♥ ♣r❛❝t✐❝❛❧ ❝♦♠♣✉t❛t✐♦♥s ✐t ✐s ♠♦r❡ ❝♦♥✈❡♥✐❡♥t t♦ ❤❛✈❡ r❡❛❧ ♣♦❧❡s✱ ✐✳❡✳✱  = [−∞, 0)✳ ❍♦✇ t♦ s❡❧❡❝t {ξ1, ξ2, . . .} ⊂ ❄ ❯s❡ st❛♥❞❛r❞ t♦♦❧ t♦ s♦❧✈❡ ♥♦♥st❛♥❞❛r❞ ❛♣♣r♦①✐♠❛t✐♦♥ ♣r♦❜❧❡♠✿ ❆ss✉♠❡ t♦ ❤❛✈❡ ❛ s✐♥❣❧❡ r❡♣❡❛t❡❞ ♣♦❧❡ ξ✳ ❚❤❡♥ rτ

m(z) =

pτ

m−1(z)

qm−1(z) = pτ

m−1(z)

(z − ξ)m−1 = pτ

m−1(

z),

  • z = (z − ξ)−1,

✐✳❡✳✱ ✇❡ ❤❛✈❡ ❛ ♣♦❧②♥♦♠✐❛❧ ♣r♦❜❧❡♠✿ ❛♠♦♥❣ p ∈ Pm−1 ♠✐♥✐♠✐③❡

  • f τ(

z−1 + ξ) − p( z)

  •  ,
  •  =

z : z ∈  .

slide-30
SLIDE 30

❆♣♣❧② ❲❛❧s❤✬s t❤❡♦r② ♦♥ ♣♦❧②♥♦♠✐❛❧ ❛♣♣r♦①✐♠❛t✐♦♥ t♦ ♦❜t❛✐♥ t❤❡ ❛s②♠♣t♦t✐❝ ❝♦♥✈❡r❣❡♥❝❡ r❛t❡ R1(ξ, τ)✳ ❋♦r t❤❡ tr❛♥s❢❡r ❢✉♥❝t✐♦♥✿ R1(ξ, τ) =

  • 1 +
  • 8d3/4 + 4d1/2 +
  • 8d1/4

1 + d −1/2 , d = −τ2/ξ2.

10 10

1

10

2

0.2 0.4 0.6 0.8 1 imag(τ) R1(ξ,τ)

ξ = −1

slide-31
SLIDE 31

❆♣♣❧② ❲❛❧s❤✬s t❤❡♦r② ♦♥ ♣♦❧②♥♦♠✐❛❧ ❛♣♣r♦①✐♠❛t✐♦♥ t♦ ♦❜t❛✐♥ t❤❡ ❛s②♠♣t♦t✐❝ ❝♦♥✈❡r❣❡♥❝❡ r❛t❡ R1(ξ, τ)✳ ❋♦r t❤❡ tr❛♥s❢❡r ❢✉♥❝t✐♦♥✿ R1(ξ, τ) =

  • 1 +
  • 8d3/4 + 4d1/2 +
  • 8d1/4

1 + d −1/2 , d = −τ2/ξ2.

10 10

1

10

2

0.2 0.4 0.6 0.8 1 imag(τ) R1(ξ,τ)

ξ = −10

slide-32
SLIDE 32

❆♣♣❧② ❲❛❧s❤✬s t❤❡♦r② ♦♥ ♣♦❧②♥♦♠✐❛❧ ❛♣♣r♦①✐♠❛t✐♦♥ t♦ ♦❜t❛✐♥ t❤❡ ❛s②♠♣t♦t✐❝ ❝♦♥✈❡r❣❡♥❝❡ r❛t❡ R1(ξ, τ)✳ ❋♦r t❤❡ tr❛♥s❢❡r ❢✉♥❝t✐♦♥✿ R1(ξ, τ) =

  • 1 +
  • 8d3/4 + 4d1/2 +
  • 8d1/4

1 + d −1/2 , d = −τ2/ξ2.

10 10

1

10

2

0.2 0.4 0.6 0.8 1 imag(τ) R1(ξ,τ)

ξ = −100

slide-33
SLIDE 33

❈♦♥s✐❞❡r p ♣♦❧❡s {ξ1, . . . , ξp} r❡♣❡❛t❡❞ ❝②❝❧✐❝❛❧❧②✳ ❚❤❡ ♣r♦❞✉❝t ❢♦r♠ ♦❢ t❤❡ ❡rr♦r f τ(z) − rτ

m(z) = qm−1(τ) qm−1(z) ωm(z) ωm(τ)

z − τ ❛❧❧♦✇s t♦ ❝♦♥❝❧✉❞❡ t❤❛t R({ξ1, . . . , ξp}, τ) =

p

  • j=1

R1(ξj, τ)1/p ✐s t❤❡ ❛s②♠♣t♦t✐❝ ❝♦♥✈❡r❣❡♥❝❡ r❛t❡ ❢♦r t❤✐s ♣♦❧❡ s❡q✉❡♥❝❡✳ ⇒ ❋✐♥❞ {ξ∗

1 , . . . , ξ∗ m} ♠✐♥✐♠✐③✐♥❣ t❤❡ ✇♦rst✲❝❛s❡ r❛t❡

maxτ∈T R({ξ1, . . . , ξm}, τ).

slide-34
SLIDE 34

❋✐♥❞ ♦♣t✐♠✳ ♣♦❧❡s ❜② ♥♦♥♥❡❣✳ ♠✐♥✐♠✐③❛t✐♦♥ e − Mx∞✳ ❍❡r❡ ✐s t❤❡ ♦♣t✐♠❛❧ ♦✈❡r❛❧❧✲❝♦♥✈❡r❣❡♥❝❡ r❛t❡ ♦♥ T = i[1, c] ❞❡♣❡♥❞✐♥❣ ♦♥ c✳

10 10

1

10

2

10

3

10

4

10

5

0.5 1 1.5 2 parameter spread c convrate R(c) real poles

slide-35
SLIDE 35

❋✐♥❞ ♦♣t✐♠✳ ♣♦❧❡s ❜② ♥♦♥♥❡❣✳ ♠✐♥✐♠✐③❛t✐♦♥ e − Mx∞✳ ❍❡r❡ ✐s t❤❡ ♦♣t✐♠❛❧ ♦✈❡r❛❧❧✲❝♦♥✈❡r❣❡♥❝❡ r❛t❡ ♦♥ T = i[1, c] ❞❡♣❡♥❞✐♥❣ ♦♥ c✳

10 10

1

10

2

10

3

10

4

10

5

0.5 1 1.5 2 parameter spread c convrate R(c) real poles imag poles

slide-36
SLIDE 36

❋✐♥❞ ♦♣t✐♠✳ ♣♦❧❡s ❜② ♥♦♥♥❡❣✳ ♠✐♥✐♠✐③❛t✐♦♥ e − Mx∞✳ ❍❡r❡ ✐s t❤❡ ♦♣t✐♠❛❧ ♦✈❡r❛❧❧✲❝♦♥✈❡r❣❡♥❝❡ r❛t❡ ♦♥ T = i[1, c] ❞❡♣❡♥❞✐♥❣ ♦♥ c✳

10 10

1

10

2

10

3

10

4

10

5

0.5 1 1.5 2 parameter spread c convrate R(c) real poles imag poles log(imag)/log(real)

❊①✿ ■❢ s♦❧✈❡ ♦❢ ❝♦♠♣❧❡① s②st❡♠ ✐s ✶✳✺ ❛s ❡①♣❡♥s✐✈❡ ❛s r❡❛❧ s♦❧✈❡✱ ✉s❡ ✐♠❛❣✐♥❛r② ♣♦❧❡s ♦♥❧② ✐❢ c < 10✦

slide-37
SLIDE 37

❈♦♠♣✉t❡✿ f τ(A)b = (A − τI)−1b✱ A = diag(0, . . . , 1e4)✱ τ ∈ T = i[10, 1000]✱ b = randn✳ 10 20 30 40 50 60 70 10

−15

10

−10

10

−5

  • rder m

error τ = 10i τ = 1000i P❆■◆ ♠❡t❤♦❞

slide-38
SLIDE 38

❈♦♠♣✉t❡✿ f τ(A)b = (A − τI)−1b✱ A = diag(0, . . . , 1e4)✱ τ ∈ T = i[10, 1000]✱ b = randn✳ 10 20 30 40 50 60 70 10

−15

10

−10

10

−5

  • rder m

error τ = 10i τ = 1000i ❘❛②❧❡✐❣❤✲❘✐t③

slide-39
SLIDE 39

❈♦♠♣✉t❡✿ f τ(A)b = (A − τI)−1b✱ A = diag(0, . . . , 1e4)✱ τ ∈ T = i[10, 1000]✱ b = randn✳ 10 20 30 40 50 60 70 10

−15

10

−10

10

−5

  • rder m

error τ = 10i τ = 1000i ❘❛②❧❡✐❣❤✲❘✐t③ ⇐= s✉♣❡r❧✐♥❡❛r✦

slide-40
SLIDE 40

❈♦♠♣✉t❡✿ f τ(A)b = (A − τI)−1b✱ A = diag(0, . . . , 1e4)✱ τ ∈ T = i[10, 1000]✱ b = randn✳ 10 20 30 40 50 60 70 10

−15

10

−10

10

−5

  • rder m

error τ = 10i τ = 1000i ❘❛②❧❡✐❣❤✲❘✐t③

slide-41
SLIDE 41

❙✉♠♠❛r②

◮ ❍❛✈❡ ❝❤❛r❛❝t❡r✐③❡❞ t❤❡ ❣❡♥❡r❛❧ ❢♦r♠ ♦❢ ❛ r❛t✐♦♥❛❧ ❑r②❧♦✈

❞❡❝♦♠♣♦s✐t✐♦♥✳

◮ ❆❧❧ ❡①✐st✐♥❣ r❛t✐♦♥❛❧ ❑r②❧♦✈ ♠❡t❤♦❞s ✜t ✐♥t♦ t❤✐s ❢r❛♠❡✇♦r❦✳ ◮ Pr♦♣♦s❡ ✏P❆■◆✑ ❛s ❛♥ ❡❢✜❝✐❡♥t ❛♥❞ r♦❜✉st r❛t✐♦♥❛❧ ❑r②❧♦✈

♠❡t❤♦❞ ❢♦r ♣r♦❜❧❡♠s ✇✐t❤ ❦♥♦✇♥ s♣❡❝tr❛❧ ♣r♦♣❡rt✐❡s✳

◮ ❍❛✈❡ ♣r❡s❡♥t❡❞ s✐♠♣❧❡ ♠❡t❤♦❞ ❢♦r ✜♥❞✐♥❣ ❝♦♥str❛✐♥❡❞ ♣♦❧❡

s❡q✉❡♥❝❡s ②✐❡❧❞✐♥❣ ❛s②♠♣t♦t✐❝❛❧❧② ♦♣t✐♠❛❧ ❝♦♥✈❡r❣❡♥❝❡✳

◮ ❚❤✐s ♠❡t❤♦❞ ♠❛② ❜❡ ❛♣♣❧✐❡❞ ❢♦r ❣❡♥❡r❛❧ f ❜② ✉s✐♥❣ ❈❛✉❝❤②

✐♥t❡❣r❛❧ r❡♣r❡s❡♥t❛t✐♦♥✳

◮ ❈❛♥ ❡①♣❧❛✐♥ s✉♣❡r❧✐♥❡❛r ❝♦♥✈❡r❣❡♥❝❡ ♦❜s❡r✈❡❞ ✇✐t❤

❘❛②❧❡✐❣❤✲❘✐t③ ❡①tr❛❝t✐♦♥ ❢♦r ❍❡r♠✐t✐❛♥ ♣r♦❜❧❡♠s ✉s✐♥❣ ✇❡✐❣❤t❡❞ ♣♦t❡♥t✐❛❧ t❤❡♦r② ❬❇❡❝❦❡r♠❛♥♥✱ ●✳ ✫ ❱❛♥❞❡❜r✐❧ ✵✾❪✳

slide-42
SLIDE 42

❙❡❡ ♠② ♣♦st❡r ❢♦r ✏r❛t✐♦♥❛❧ ❘✐t③ ✈❛❧✉❡s✑ ❛♥❞ ✏✐♥❡①❛❝t s♦❧✈❡s✑✿

Rational Krylov methods and approximation of f (A)b

  • B. Beckermann‡, M. Eiermann†, O. G. Ernst†, Stefan Güttel†, and R. Vandebril‡

†TU Bergakademie Freiberg, Germany ‡Université Lille 1, France

Institut für Numerische Mathematik und Optimierung Equipe d’Analyse Numérique et d’Optimisation

Matrix Functions

Given a square matrix A of size N × N, a vector b of length N and a scalar function f (z), f (A)b := p(A)b, where p ∈ N−1 is a polynomial of degree ≤ N − 1 that Hermite- interpolates f at the eigenvalues of A. In typical applications the matrix A is large and sparse. Some Applications

  • f (z) = (z − iω)−1: model reduction in the frequency domain,
  • f (z) = exp(−tz): time-integration of linear ODE’s, exponential

integrators, e.g., in geophysics or chemistry,

  • f (z) = tz: simulation of Brownian motion of molecules or

sampling from Gaussian Markov random fields,

  • f (z) = sign(z): simulations in quantum chromodynamics.

Rational Krylov Spaces

Definition: Given a sequence of polynomials qm−1(z) =

m−1

  • j=1

ξj=∞

(z − ξj), m = 1,2,..., where ξj ∈ C \ Λ(A). Then the associated rational Krylov spaces

  • f order m are defined as

m(A,b) := qm−1(A)−1m(A,b), where m(A,b) = span{A0b,A1b,...,Am−1b }. Properties: Let M be the invariance index of (A, ). Then

Rational Ritz values

... are the eigenvalues of the Rayleigh quotient Am = Q∗

mAQm,

denoted by Θ = {θ1,...,θm}. Let A be Hermitian. Then the θk’s lie in the spectral interval of A and interlace the eigenvalues Λ(A) = {λ1,...,λN}: (*) In any interval (θκ,θκ+1) there is at least one eigenvalue λk of A. Moreover, the rational Ritz values are zeros of orthogonal rational functions and may be characterized as (see, e.g., [2, 3]) (**) The θk’s are the zeros of the minimizer of p(A)qm−1(A)−1b among all monic p ∈ ∞

m .

Logarithmic potential theory can explain the asymptotic distribu- tion of the rational Ritz values. Therefore we consider

  • a sequence of Hermitian matrices {AN}, each of size N × N,

whose eigenvalue counting measures converge to a Borel prob- ability measure σ in the weak-* sense,

  • a sequence of vectors {bN}, each of length N,
  • a ray sequence of integers {mN} such that

mN/N → t ∈ (0,1) as N → +∞,

  • a sequence of polynomials {qN}, each of degree mN − 1, whose

zero counting measures converge to a Borel measure ν, ν = t,

  • the sequence {ΘN} of rational Ritz values of order mN.

Tools from Potential Theory Associated with a (signed) Borel measure µ1 is the logarithmic potential Uµ1(z) :=

  • 1

log|x − z| dµ1(x).

−1 −0.5 0.5 1 −0.2 0.2 0.4 0.6 0.8

ν Σt* Uµ−ν J µ σ U−ν F

t = 0.7

Inexact solves & error estimators

In each iteration of the rational Arnoldi method a linear system of the form (A− ξjI)xj = qj is solved. If the residuals are collected in a matrix Rm, then (1) becomes AQm+1Km = Qm+1Hm + Rm. (2) Setting Em := −RmKm†Q∗

m+1, we observe that we have computed

an exact Arnoldi decomposition (A+ Em)Qm+1Km = Qm+1Hm for the matrix A+ Em. The Rayleigh quotient Am computed from the data Km and Hm satisfies

  • Am = Q∗

m(A+ Em)Qm

= Q∗

mAQm + Q∗ m(−RmKm†Q∗ m+1)Qm

= Am −Q∗

mRmKm†Im.

Here, Am := Q∗

mAQm is referred to as the corrected Rayleigh quo-

tient, because it is a compression of A instead of A+ Em. It can be computed from Am without explicit projection, only by additional inner-products Q∗

mRm.

We now decompose the error f (A)b − fm ≤ f (A)b − f (A+ Em)b

  • sensitivity error

+ f (A+ Em)b − fm

  • approximation error

, and estimate sensitivity error ≈ f ( Am)Q∗

mb − f (

Am)Q∗

mb .

It is advisable to terminate the rational Arnoldi method if the ap- proximation error falls below the sensitivity error, because after