A rational Krylov method for non-linear matrix problems Karl - - PowerPoint PPT Presentation

A rational Krylov method for non-linear matrix problems

Karl Meerbergen (Joint work with Roel Van Beeumen and Wim Michiels)

KU Leuven

SIMAX – Valencia

Motivation

Solution of linear eigenvalue problem Ax = λBx is pretty well understood

◮ Appealing method: shift-and-invert Arnoldi ◮ Reliable for computing several eigenvalues near a shift point. ◮ Rational Krylov: more than one interpolation points = poles.

Polynomial eigenvalue problems

◮ We can use Krylov methods (when we ‘linearize’)

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 2 / 24

Motivation

Solution of linear eigenvalue problem Ax = λBx is pretty well understood

◮ Appealing method: shift-and-invert Arnoldi ◮ Reliable for computing several eigenvalues near a shift point. ◮ Rational Krylov: more than one interpolation points = poles.

Polynomial eigenvalue problems

◮ We can use Krylov methods (when we ‘linearize’)

Solution of non-linear eigenvalue problem A(λ)x = 0

◮ Residual inverse iteration and Newton iterations are locally convergent

methods

◮ Globally convergent methods by using polynomial approximation of

A(λ)

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 2 / 24

Outline

1

Motivation

2

Polynomial eigenvalue problem

3

Nonlinear eigenvalue problem

4

Rational Krylov

5

A statement on moment matching and projection

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 3 / 24

Polynomial eigenvalue problem

Polynomial eigenvalue problem: A0x + λA1x + · · · + λNANx = 0 Transformed to

(A − λB)x =           A0 A1 · · · Ap−1 I ... I      − λ      · · · · · · −AN I ... . . . I                x λx . . . λp−1x      =

Is a linear problem, so we can use Arnoldi’s method Is a well known technique.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 4 / 24

Linearization for polynomial eigenvalue problems

(Shift and) invert transform A−1B: λ−1      x λx . . . λp−1x      =      S1 · · · · · · Sp I ... . . . I           x λx . . . λp−1x      Krylov sequence:    v . . .   

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 5 / 24

Linearization for polynomial eigenvalue problems

(Shift and) invert transform A−1B: λ−1      x λx . . . λp−1x      =      S1 · · · · · · Sp I ... . . . I           x λx . . . λp−1x      Krylov sequence:    v . . .         w = S1v v . . .     

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 5 / 24

Linearization for polynomial eigenvalue problems

(Shift and) invert transform A−1B: λ−1      x λx . . . λp−1x      =      S1 · · · · · · Sp I ... . . . I           x λx . . . λp−1x      Krylov sequence:    v . . .         w v . . .             t = S1w + S2v w v . . .       

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 5 / 24

Linearization for polynomial eigenvalue problems

(Shift and) invert transform A−1B: λ−1      x λx . . . λp−1x      =      S1 · · · · · · Sp I ... . . . I           x λx . . . λp−1x      Krylov sequence:    v . . .         w v . . .             t w v . . .                 u = S1t + S2w + S3v t w v . . .         

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 5 / 24

Linearization for polynomial eigenvalue problems

Shift-and-invert transform: λ−1      x λx . . . λp−1x      =      S1 · · · · · · Sp I ... . . . I           x λx . . . λp−1x      Krylov sequence:               v1 v2 v3 v4 · · · vk v1 v2 v3 vk−1 . . . v1 v2 vk−2 . . . v1 . . . . . . . . . . . . v1              

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 6 / 24

Nonlinear eigenvalue problem

Straightforward idea: choose an approximation of a fixed degree and solve the corresponding polynomial eigenvalue problem:

◮ Polynomial expansion, e.g. Taylor, Chebyshev, . . .

A(λ) ≈ A0p0(λ) + · · · + ANpN(λ) [Amiraslani, Corless, Lancaster, 2009], [Effenberger, Kressner, 2011]

◮ Rational expansion, e.g. (Chebyshev) Pad´

e A(λ) ≈ A0 α0 λ − σ0 + · · · + AN αN λ − σN [Su & Bai, 2011]

◮ Spectral discretization

A(λ) ≈ A0 + λA0p0(λ) + · · · + ANpN(λ) α0p0(λ) + · · · + αNpN(λ) [Trefethen 2000], [Michiels, Niculescu 2007], [Breda, Maset, Vermiglio 2005]

In this talk, we explore a ‘flexible’ choice of N and pj.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 7 / 24

Taylor-Arnoldi

[Jarlebring, Meerbergen, Michiels, 2010] Power sequence on Taylor series A(λ) ≈ A0 + λA1 + λ2A2 + · · · For k iterations of Arnoldi, use Companion linearization with N > k Arnoldi:    v . . .    Breakthrough for solving the delay (equation) eigenvalue problem

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 8 / 24

Taylor-Arnoldi

[Jarlebring, Meerbergen, Michiels, 2010] Power sequence on Taylor series A(λ) ≈ A0 + λA1 + λ2A2 + · · · For k iterations of Arnoldi, use Companion linearization with N > k Arnoldi:    v . . .         w = −A−1

0 A1v

v . . .      Breakthrough for solving the delay (equation) eigenvalue problem

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 8 / 24

Taylor-Arnoldi

[Jarlebring, Meerbergen, Michiels, 2010] Power sequence on Taylor series A(λ) ≈ A0 + λA1 + λ2A2 + · · · For k iterations of Arnoldi, use Companion linearization with N > k Arnoldi:    v . . .         w v . . .             t = −A−1

0 (A1w + A2)v

w v . . .        Breakthrough for solving the delay (equation) eigenvalue problem

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 8 / 24

Taylor-Arnoldi

[Jarlebring, Meerbergen, Michiels, 2010] Power sequence on Taylor series A(λ) ≈ A0 + λA1 + λ2A2 + · · · For k iterations of Arnoldi, use Companion linearization with N > k Arnoldi:    v . . .         w v . . .             t w v . . .                 u = −A−1

0 (A1t + A2w + A3)v

t w v . . .          Breakthrough for solving the delay (equation) eigenvalue problem

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 8 / 24

Other polynomials

Use polynomials in order to

◮ avoid derivatives or ◮ interpolation in more than one point.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 9 / 24

Other polynomials

Use polynomials in order to

◮ avoid derivatives or ◮ interpolation in more than one point.

Newton polynomial Lagrange polynomial Chebyshev polynomial

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 9 / 24

Other polynomials

Use polynomials in order to

◮ avoid derivatives or ◮ interpolation in more than one point.

Newton polynomial Lagrange polynomial Chebyshev polynomial A (partly) successful generalization of Taylor-Arnoldi to rational Krylov

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 9 / 24

Newton polynomial basis

Newton polynomials: N0 := 1, N1 := (λ − σ0), N2 := (λ − σ0)(λ − σ1), . . . A(λ) = A0N0 + A1N1 + · · · with

◮ A0 = [A0] = A(σj) and ◮ Aj = [A0, . . . , Aj] = [A1, . . . , Aj] − [A0, . . . , Aj−1]

σj − σ0 (divided difference).

σ0, σ1, . . . are interpolation points. Linearization:

              A0 A1 · · · AN σ0I I σ1I I ... ... σN−2I I        − λ       · · · · · · I . . . ... . . . I                      x N1(λ)x N2(λ)x . . . . . .         = 0

A−1B is not a block upper Hessenberg matrix

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 10 / 24

Rational Krylov

Shift-and-invert step: (A − σjB)w = Bv First step: multiply by B: shift vector downwards B · v = w       · · · · · · I . . . ... . . . I       ·        v1 v2 . . . vN−1 vN        =        v1 v2 . . . vN−1        Then apply (A − σjB)−1: A − σjB =        A0 A1 · · · AN (σ0 − σj)I I (σ1 − σj)I I ... ... (σN−1 − σj)I I       

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 11 / 24

Rational Krylov

µj

i = σi − σj.

Shift-and-invert-step:

        A0 A1 A2 A3 A4 A5 σ0I I σ1I I σ2I I σ3I I σ4I I        

−1 

       ⋆         =         ⋆ ⋆ ⋆ ⋆ ⋆ ⋆        

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors This is rational Krylov with poles σ1, σ2, . . . Linear recurrence relation of the form: AVk+1Hk = BVk+1Lk

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 12 / 24

Rational Krylov

µj

i = σi − σj.

Shift-and-invert-step: A − σ1B

        A0 A1 A2 A3 A4 A5 µ1

0I

I I µ1

2I

I µ1

3I

I µ1

4I

I        

−1 

       ⋆         =         ⋆ ⋆        

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors This is rational Krylov with poles σ1, σ2, . . . Linear recurrence relation of the form: AVk+1Hk = BVk+1Lk

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 12 / 24

Rational Krylov

µj

i = σi − σj.

Shift-and-invert-step: A − σ2B

        A0 A1 A2 A3 A4 A5 µ2

0I

I µ2

1I

I I µ2

3I

I µ2

4I

I        

−1 

       ⋆ ⋆         =         ⋆ ⋆ ⋆        

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors This is rational Krylov with poles σ1, σ2, . . . Linear recurrence relation of the form: AVk+1Hk = BVk+1Lk

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 12 / 24

Rational Krylov

µj

i = σi − σj.

Shift-and-invert-step: A − σ3B

        A0 A1 A2 A3 A4 A5 µ3

0I

I µ3

1I

I µ3

2I

I I µ3

4I

I        

−1 

       ⋆ ⋆ ⋆         =         ⋆ ⋆ ⋆ ⋆        

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors This is rational Krylov with poles σ1, σ2, . . . Linear recurrence relation of the form: AVk+1Hk = BVk+1Lk

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 12 / 24

Rational Krylov

µj

i = σi − σj.

Shift-and-invert-step: A − σ4B

        A0 A1 A2 A3 A4 A5 µ4

0I

I µ4

1I

I µ4

2I

I µ4

3I

I I        

−1 

       ⋆ ⋆ ⋆ ⋆         =         ⋆ ⋆ ⋆ ⋆ ⋆        

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors This is rational Krylov with poles σ1, σ2, . . . Linear recurrence relation of the form: AVk+1Hk = BVk+1Lk

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 12 / 24

Rational Krylov

With constant continuation vector v1, (A − σ1B)−1Bv1, (A − σ2B)−1Bv1, . . . On each iteration, we must solve the system: A(σj)wj = [A(σ0), A(σj)]v1 So, only two matrices have to be stored: A(σ0) and A(σj). Note: [A(σ0), A(σj)] ≈ A′(σj) Recurrence relation has small residual norm (rounding errors). Caveat: choice of continuation vector may be unstable. When the σj become Ritz values or converged Ritz value show up again, be worried.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 13 / 24

Rational Krylov

Usual choice of continuation vector v1, (A − σ1B)−1Bv1, (A − σ2B)−1Bv2, . . . where vj is the last subspace vector. On each iteration, we must solve the system: A(σj)wj =

j

• i=1

Aiwj−1,i So, k divided differences have to be stored or computed Note: j

i=1 Aiwj−1,i ≈ A′(σj)w0,i

Recurrence relation has small residual norm (rounding errors). Is known to be stable in practice.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 14 / 24

Newton basis for scalar equation

Compute roots of α0 + α1N1(z) + · · · + αdNd(z) = 0 Linearization in iteration k (pole σk):                          

α0 · · · αd · · · σ0 1 . . . ... 1 σd 1 ... ... . . . σk−1 1

         − λ        

1 ... 1 ... ... 1

                                 

x0 . . . xd . . . . . . xk

         = 0 The d roots λ are eigenvalues ∞ is a defective eigenvalue (multiplicity k + 1 − d)

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 15 / 24

Newton basis for scalar equation

RKS Recurrence relation: AVkHk − BVkKk = fkeT

k

with fk = (A − σkB)vk+1 = 0 In this case, Vk = Ik+1,k Last RKS iteration: (A − σkB)−1B = C dT

• with

C ∈ Ck×k From recurrence relation: (A − σkB)−1BVk − VkHk(Kk − σkHk)−1 = vk+1gT

k

C − Hk(Kk − σkHk)−1 = Thus the Ritz values computed from Kkz = λHkz are roots of the polynomial RKS residual is not zero!

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 16 / 24

Lagrange polynomial basis

Lagrange polynomials: Lj := ωj

• i=j

(λ − σi) , ωj =  

i=j

(σj − σi)  

−1

, j = 0, . . . , N A(λ) = A0L0 + A1L1 + · · · with

◮ Aj = A(σj), j = 0, . . . , N and ◮ σ0, σ1, . . . the interpolation points.

Linearization:

           A0 · · · AN σ0I −ω0 ... . . . σN−1I −ωNI      − λ       · · · · · · I . . . ... . . . I                  L0x L1x . . . LNx      = 0

A−1B is not a block upper Hessenberg matrix

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 17 / 24

Rational Krylov

Apply shift: (A − σjB)u = λBu has the following form:

                A0 A1 · · · AN ⋆ −ω0I . . . ... ⋆ −ωNI         − λ        · · · · · · I I ... . . . I                       x L1x L2x . . . LNx        = 0

(A − σjB)−1B applied to vj respects the desired zero structure in the vectors The method does not see the last column of A (ωj): does not work for finding roots of a polynomial.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 18 / 24

Chebyshev series

A(λ) = A0T0 + A1T1 + · · · Linearization:

                    −A0 −A1 · · · ⋆ I I I ... ... ... ... ... I I           − λ        · · · · · · · · · −AN I 2I ... . . . 2I                  x = 0

No block upper Hessenberg structure Arnoldi iterations require computation of A(0), A′(0) expressed in Chebyshev basis. i.e. moment matching requires the full Chebyshev expansion (use a functional framework [Jarlebring, Michiels, M. 2011]) Shifting as for Newton basis does not help.

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 19 / 24

Numerical example: Hopf bifurcation

Nonlinear equation: ˙ x = f (x, λ) Hopf bifurcation of dynamical system: A(λ)u = iµu with µ real λ is an eigenvalue of (A(λ) ⊗ I + I ⊗ A(λ)x = 0 Usually: continuation method for λ Modified Olmstead model (n = 40) of the form: (A + λEu + (λ + 1)−1Eu + (λ + 2)−1/2Ev)u = 0 Pick ten poles in a Leja fashion in [6, 10], then use Ritz value. For related work, see [Guckerheimer et al. 1993, 1997], [Govaerts, 2000] [M. & Spence 2010], [Elman, M., Spence, Wu, 2012], [Elman, Wu, 2012]

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 20 / 24

Hopf bifurcations

20 iterations rational Krylov λ B−1

N ANx − λx2

A(λ)x2 8.127 1.0 · 10−12 6.0 · 10−12 8.127 1.1 · 10−11 7.4 · 10−11 7.405 8.6 · 10−11 4.3 · 10−10 Finding more than three eigenvalues is not possible for this problem with this method

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 21 / 24

Gun problem

Problem from the NLEVP collection (K − λM + i

• λ − σ2

1W1 + i

• λ − σ2

2W2)x = 0

where K has order 9956 and rank(W1) + rank(W2) = 84. The goal is to compute all eigenvalues in a halve circle in the λ-plane We chose five interpolation points equally distributed in the λ-plane

100 150 200 250 300 350 20 40 60 80 100 120 Re √ λ Im √ λ 10 20 30 40 50 60 10−16 10−14 10−12 10−10 10−8 10−6 10−4 10−2 100 iteration relative residual norm E(λ, x)

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 22 / 24

Moment matching

Linear system: A(s)x = f Krylov methods match moments of the transfer function AN(s)x = f :

◮ Arnoldi on Companion form: subspace spanned by x(j)(0),

j = 0, . . . , N.

◮ Rational Krylov: subspace spanned by x(σj) = A−1

N (σj)f .

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 23 / 24

Conclusions

Also attend talk by Roel Van Beeumen (CP38, Thursday 17:50) Advantages of rational Krylov:

◮ Low degree polynomials ◮ With or without derivatives ◮ Exploit structure, use symbolic tools (See talk Roel Van Beeumen) ◮ Not slower than Newton

Disadvantages of rational Krylov:

◮ Derivatives required for computing many eigenvalues (Hermite

interpolation)

◮ Newton polynomials require divided differences

• K. Meerbergen (KU Leuven)

Nonlinear Rational Krylov SIMAX – Valencia 24 / 24