Krylov subspace methods for eigenvalue problems David S. Watkins - - PowerPoint PPT Presentation

Krylov subspace methods for eigenvalue problems

David S. Watkins

watkins@math.wsu.edu

Department of Mathematics Washington State University

October 16, 2008 – p.

Problem: Linear Elasticity

October 16, 2008 – p.

Problem: Linear Elasticity

Elastic Deformation (3D, anisotropic, composite materials)

October 16, 2008 – p.

Problem: Linear Elasticity

Elastic Deformation (3D, anisotropic, composite materials) Singularities at cracks, interfaces

October 16, 2008 – p.

Problem: Linear Elasticity

Elastic Deformation (3D, anisotropic, composite materials) Singularities at cracks, interfaces Lamé Equations (spherical coordinates)

October 16, 2008 – p.

Problem: Linear Elasticity

Elastic Deformation (3D, anisotropic, composite materials) Singularities at cracks, interfaces Lamé Equations (spherical coordinates) Separate radial variable.

October 16, 2008 – p.

Problem: Linear Elasticity

Elastic Deformation (3D, anisotropic, composite materials) Singularities at cracks, interfaces Lamé Equations (spherical coordinates) Separate radial variable. Get quadratic eigenvalue problem.

(λ2M + λG + K)v = 0 M∗ = M > 0 G∗ = −G K∗ = K < 0

October 16, 2008 – p.

Linear Elasticity, Continued

Discretize θ and ϕ variables.

October 16, 2008 – p.

Linear Elasticity, Continued

Discretize θ and ϕ variables. (finite element method)

October 16, 2008 – p.

Linear Elasticity, Continued

Discretize θ and ϕ variables. (finite element method)

(λ2M + λG + K)v = 0 MT = M > 0 GT = −G KT = K < 0

matrix quadratic eigenvalue problem (large, sparse)

October 16, 2008 – p.

Linear Elasticity, Continued

Discretize θ and ϕ variables. (finite element method)

(λ2M + λG + K)v = 0 MT = M > 0 GT = −G KT = K < 0

matrix quadratic eigenvalue problem (large, sparse) Find few smallest eigenvalues (and corresponding eigenvectors).

October 16, 2008 – p.

Linear Elasticity, Continued

Discretize θ and ϕ variables. (finite element method)

(λ2M + λG + K)v = 0 MT = M > 0 GT = −G KT = K < 0

matrix quadratic eigenvalue problem (large, sparse) Find few smallest eigenvalues (and corresponding eigenvectors). Respect the structure. (symmetric/skew-symmetric)

October 16, 2008 – p.

Hamiltonian Structure

October 16, 2008 – p.

Hamiltonian Structure

October 16, 2008 – p.

Reduction to First Order

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0 w = λv,

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0 w = λv, Mw = λMv

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0 w = λv, Mw = λMv

• −K

−M v w

• − λ
• G

M −M v w

• = 0

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0 w = λv, Mw = λMv

• −K

−M v w

• − λ
• G

M −M v w

• = 0

Ax − λBx = 0

October 16, 2008 – p.

Reduction to First Order

λ2Mv + λGv + Kv = 0 w = λv, Mw = λMv

• −K

−M v w

• − λ
• G

M −M v w

• = 0

Ax − λBx = 0

symmetric/skew-symmetric

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

B = RTJR

• J =
• I

−I

• sometimes easy, always possible

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

B = RTJR

• J =
• I

−I

• sometimes easy, always possible

A − λRTJR

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

B = RTJR

• J =
• I

−I

• sometimes easy, always possible

A − λRTJR R−T AR−1 − λJ

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

B = RTJR

• J =
• I

−I

• sometimes easy, always possible

A − λRTJR R−T AR−1 − λJ JT R−TAR−1 − λI

October 16, 2008 – p.

Reduction to Hamiltonian Matrix

A − λB

(symmetric/skew-symmetric)

B = RTJR

• J =
• I

−I

• sometimes easy, always possible

A − λRTJR R−T AR−1 − λJ JT R−TAR−1 − λI H = JT R−TAR−1 is Hamiltonian.

October 16, 2008 – p.

in our case . . .

October 16, 2008 – p.

in our case . . .

B =

• G

M −M

• October 16, 2008 – p.

in our case . . .

B =

• G

M −M

• B = RTJR =
• I

−1

2G

M I −I I

1 2G

M

• October 16, 2008 – p.

in our case . . .

B =

• G

M −M

• B = RTJR =
• I

−1

2G

M I −I I

1 2G

M

• H =
• I

−1

2G

I M−1 −K I −1

2G

I

• October 16, 2008 – p.

in our case . . .

B =

• G

M −M

• B = RTJR =
• I

−1

2G

M I −I I

1 2G

M

• H =
• I

−1

2G

I M−1 −K I −1

2G

I

• Do not compute H explicitly.

October 16, 2008 – p.

in our case . . .

B =

• G

M −M

• B = RTJR =
• I

−1

2G

M I −I I

1 2G

M

• H =
• I

−1

2G

I M−1 −K I −1

2G

I

• Do not compute H explicitly. (nor M−1)

October 16, 2008 – p.

Working with H

October 16, 2008 – p.

Working with H

Krylov subspace methods

October 16, 2008 – p.

Working with H

Krylov subspace methods just need to apply the operator:

October 16, 2008 – p.

Working with H

Krylov subspace methods just need to apply the operator: x → Hx

October 16, 2008 – p.

Working with H

Krylov subspace methods just need to apply the operator: x → Hx

H =

• I

−1

2G

I M−1 −K I −1

2G

I

• October 16, 2008 – p.

Working with H

Krylov subspace methods just need to apply the operator: x → Hx

H =

• I

−1

2G

I M−1 −K I −1

2G

I

• October 16, 2008 – p.

Working with H

Krylov subspace methods just need to apply the operator: x → Hx

H =

• I

−1

2G

I M−1 −K I −1

2G

I

• H−1

=

• I

1 2G

I (−K)−1 M I

1 2G

I

• October 16, 2008 – p.

Working with H−1

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

Do not compute (−K)−1

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

Do not compute (−K)−1 Cholesky decomposition:

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

Do not compute (−K)−1 Cholesky decomposition: (−K) = RTR To compute w = −K−1v,

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

Do not compute (−K)−1 Cholesky decomposition: (−K) = RTR To compute w = −K−1v, Solve (−K)w = v.

October 16, 2008 – p.

Working with H−1

H−1 =

• I

1 2G

I (−K)−1 M I

1 2G

I

• x → H−1x

Do not compute (−K)−1 Cholesky decomposition: (−K) = RTR To compute w = −K−1v, Solve (−K)w = v.

RT Rw = v

Backsolve!

October 16, 2008 – p.

K =

October 16, 2008 – p. 1

K =

20 40 60 80 100 120 140 20 40 60 80 100 120 140 nz = 670

October 16, 2008 – p. 1

R =

October 16, 2008 – p. 1

R =

20 40 60 80 100 120 140 20 40 60 80 100 120 140 nz = 896

October 16, 2008 – p. 1

K−1 =

October 16, 2008 – p. 1

K−1 =

20 40 60 80 100 120 140 20 40 60 80 100 120 140 nz = 21316

October 16, 2008 – p. 1

Second Application

October 16, 2008 – p. 1

Second Application

Nonlinear Optics

October 16, 2008 – p. 1

Second Application

Nonlinear Optics Schrödinger eigenvalue problem

October 16, 2008 – p. 1

Second Application

Nonlinear Optics Schrödinger eigenvalue problem

−2 2m ∇2ψ + V ψ = λψ

October 16, 2008 – p. 1

Second Application

Nonlinear Optics Schrödinger eigenvalue problem

−2 2m ∇2ψ + V ψ = λψ

Solve numerically (finite elements)

October 16, 2008 – p. 1

Second Application

Nonlinear Optics Schrödinger eigenvalue problem

−2 2m ∇2ψ + V ψ = λψ

Solve numerically (finite elements)

Kv = λMv K = KT > 0, M = MT > 0

October 16, 2008 – p. 1

Second Application

Nonlinear Optics Schrödinger eigenvalue problem

−2 2m ∇2ψ + V ψ = λψ

Solve numerically (finite elements)

Kv = λMv K = KT > 0, M = MT > 0

Matrices are large and sparse.

October 16, 2008 – p. 1

Kv = λMv

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors.

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv R−T MR−1(Rv) = λ−1(Rv)

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv R−T MR−1(Rv) = λ−1(Rv) A = R−TMR−1, AT = A > 0

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv R−T MR−1(Rv) = λ−1(Rv) A = R−TMR−1, AT = A > 0 x → Ax

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv R−T MR−1(Rv) = λ−1(Rv) A = R−TMR−1, AT = A > 0 x → Ax

backsolve

October 16, 2008 – p. 1

Kv = λMv

Want few smallest eigenvalues and associated eigenvectors. Invert the problem.

K = RT R RT Rv = λMv R−T MR−1(Rv) = λ−1(Rv) A = R−TMR−1, AT = A > 0 x → Ax

backsolve Do not form A explicitly.

October 16, 2008 – p. 1

What our applications have in common

October 16, 2008 – p. 1

What our applications have in common

large, sparse matrices

October 16, 2008 – p. 1

What our applications have in common

large, sparse matrices use of matrix factorization (Cholesky decomposition)

October 16, 2008 – p. 1

What our applications have in common

large, sparse matrices use of matrix factorization (Cholesky decomposition) some kind of structure

October 16, 2008 – p. 1

What our applications have in common

large, sparse matrices use of matrix factorization (Cholesky decomposition) some kind of structure . . . not enough time to discuss this

October 16, 2008 – p. 1

Classification

• f Eigenvalue Problems

October 16, 2008 – p. 1

Classification

• f Eigenvalue Problems

small

October 16, 2008 – p. 1

Classification

• f Eigenvalue Problems

small medium

October 16, 2008 – p. 1

Classification

• f Eigenvalue Problems

small medium large

October 16, 2008 – p. 1

Small Matrices

October 16, 2008 – p. 1

Small Matrices

store conventionally

October 16, 2008 – p. 1

Small Matrices

store conventionally similarity transformations

October 16, 2008 – p. 1

Small Matrices

store conventionally similarity transformations

QR algorithm

October 16, 2008 – p. 1

Small Matrices

store conventionally similarity transformations

QR algorithm

get all eigenvalues/vectors

October 16, 2008 – p. 1

Small Matrices

store conventionally similarity transformations

QR algorithm

get all eigenvalues/vectors

n ≈ 103

October 16, 2008 – p. 1

Medium Matrices

October 16, 2008 – p. 1

Medium Matrices

store as sparse matrix

October 16, 2008 – p. 1

Medium Matrices

store as sparse matrix no similarity transformations

October 16, 2008 – p. 1

Medium Matrices

store as sparse matrix no similarity transformations matrix factorization okay

October 16, 2008 – p. 1

Medium Matrices

store as sparse matrix no similarity transformations matrix factorization okay shift and invert

October 16, 2008 – p. 1

20 40 60 80 100 120 140 20 40 60 80 100 120 140 nz = 670

October 16, 2008 – p. 1

20 40 60 80 100 120 140 20 40 60 80 100 120 140 nz = 1815

October 16, 2008 – p. 2

Medium Matrices, Continued

October 16, 2008 – p. 2

Medium Matrices, Continued

store matrix factor as sparse matrix

October 16, 2008 – p. 2

Medium Matrices, Continued

store matrix factor as sparse matrix

n ≈ 105

October 16, 2008 – p. 2

Medium Matrices, Continued

store matrix factor as sparse matrix

n ≈ 105

get selected eigenvalues/vectors

October 16, 2008 – p. 2

Medium Matrices, Continued

store matrix factor as sparse matrix

n ≈ 105

get selected eigenvalues/vectors Krylov subspace methods

October 16, 2008 – p. 2

Medium Matrices, Continued

store matrix factor as sparse matrix

n ≈ 105

get selected eigenvalues/vectors Krylov subspace methods Jacobi-Davidson methods

October 16, 2008 – p. 2

Large Matrices

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations no shift-and-invert

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations no shift-and-invert

n ≈ 107

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations no shift-and-invert

n ≈ 107

get selected eigenvalues/vectors

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations no shift-and-invert

n ≈ 107

get selected eigenvalues/vectors Krylov subspace methods

October 16, 2008 – p. 2

Large Matrices

store as sparse matrix no similarity transformations no shift-and-invert

n ≈ 107

get selected eigenvalues/vectors Krylov subspace methods Jacobi-Davidson methods

October 16, 2008 – p. 2

Krylov Subspace Methods

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1 Pick a vector v

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1 Pick a vector v

v, Av, A2v, A3v, . . .

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1 Pick a vector v

v, Av, A2v, A3v, . . .

Krylov subspace:

Kj(A, v) = span

• v, Av, A2v, . . . , Aj−1v
• October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1 Pick a vector v

v, Av, A2v, A3v, . . .

Krylov subspace:

Kj(A, v) = span

• v, Av, A2v, . . . , Aj−1v
• Look in here for approximate eigenvectors.

October 16, 2008 – p. 2

Krylov Subspace Methods

x → Ax

Example: A = R−TMR−1 Pick a vector v

v, Av, A2v, A3v, . . .

Krylov subspace:

Kj(A, v) = span

• v, Av, A2v, . . . , Aj−1v
• Look in here for approximate eigenvectors.

. . . but need better basis

October 16, 2008 – p. 2

Arnoldi Process

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

hij = Avj, vi

(Gram-Schmidt)

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

hij = Avj, vi

(Gram-Schmidt) Normalization: hj+1,j = ˆ

vj+1 , vj+1 = ˆ vj+1/hj+1,j

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

hij = Avj, vi

(Gram-Schmidt) Normalization: hj+1,j = ˆ

vj+1 , vj+1 = ˆ vj+1/hj+1,j

Collect coefficients hij

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

hij = Avj, vi

(Gram-Schmidt) Normalization: hj+1,j = ˆ

vj+1 , vj+1 = ˆ vj+1/hj+1,j

Collect coefficients hij

H4 =      h11 h12 h13 h14 h21 h22 h23 h24 h32 h33 h34 h43 h44     

October 16, 2008 – p. 2

Arnoldi Process

build orthonormal basis v1, v2, v3, . . .

jth step: ˆ vj+1 = Avj − j

i=1 vihij

hij = Avj, vi

(Gram-Schmidt) Normalization: hj+1,j = ˆ

vj+1 , vj+1 = ˆ vj+1/hj+1,j

Collect coefficients hij

H4 =      h11 h12 h13 h14 h21 h22 h23 h24 h32 h33 h34 h43 h44     

eigenvalues are Ritz values

October 16, 2008 – p. 2

Example: 479 × 479 matrix

−150 −100 −50 50 100 150 −2000 −1500 −1000 −500 500 1000 1500 2000

October 16, 2008 – p. 2

Example: 479 × 479 matrix

−150 −100 −50 50 100 150 −2000 −1500 −1000 −500 500 1000 1500 2000

October 16, 2008 – p. 2

Example: 479 × 479 matrix

−150 −100 −50 50 100 150 −2000 −1500 −1000 −500 500 1000 1500 2000

October 16, 2008 – p. 2

Example: 479 × 479 matrix

−150 −100 −50 50 100 150 −2000 −1500 −1000 −500 500 1000 1500 2000

October 16, 2008 – p. 2

For better accuracy . . .

October 16, 2008 – p. 2

For better accuracy . . . . . . take more steps.

October 16, 2008 – p. 2

For better accuracy . . . . . . take more steps. but,

October 16, 2008 – p. 2

For better accuracy . . . . . . take more steps. but, vectors take up space.

October 16, 2008 – p. 2

For better accuracy . . . . . . take more steps. but, vectors take up space. Alternate plan:

October 16, 2008 – p. 2

For better accuracy . . . . . . take more steps. but, vectors take up space. Alternate plan: Start over with a better vector.

October 16, 2008 – p. 2

Implicit Restarts

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps.

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values.

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values. Keep the best ones (e.g. 10) . . .

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values. Keep the best ones (e.g. 10) . . . . . . and associated invariant subspace.

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values. Keep the best ones (e.g. 10) . . . . . . and associated invariant subspace. Restart at step 11.

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values. Keep the best ones (e.g. 10) . . . . . . and associated invariant subspace. Restart at step 11. (neat details omitted)

October 16, 2008 – p. 3

Implicit Restarts

Take, say, 30 steps. Get 30 Ritz values. Keep the best ones (e.g. 10) . . . . . . and associated invariant subspace. Restart at step 11. (neat details omitted) Build back up to 30 and then restart again.

October 16, 2008 – p. 3

Example

October 16, 2008 – p. 3

Example

460 × 460 Hamiltonian matrix (toy problem)

October 16, 2008 – p. 3

Example

460 × 460 Hamiltonian matrix (toy problem)

asking for 24 eigenpairs

October 16, 2008 – p. 3

Example

460 × 460 Hamiltonian matrix (toy problem)

asking for 24 eigenpairs building to dimension 84

October 16, 2008 – p. 3

Example

460 × 460 Hamiltonian matrix (toy problem)

asking for 24 eigenpairs building to dimension 84 restarting with 28

October 16, 2008 – p. 3

Example

460 × 460 Hamiltonian matrix (toy problem)

asking for 24 eigenpairs building to dimension 84 restarting with 28 Hamiltonian Lanczos process

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 3

−3 −2 −1 1 2 3 −3 −2 −1 1 2 3

October 16, 2008 – p. 4

Concluding Remarks

October 16, 2008 – p. 4

Concluding Remarks

I used these methods to solve my “medium” sized eigenvalue problems.

October 16, 2008 – p. 4

Concluding Remarks

I used these methods to solve my “medium” sized eigenvalue problems. Simple ideas lead to powerful methods.

October 16, 2008 – p. 4

Concluding Remarks

I used these methods to solve my “medium” sized eigenvalue problems. Simple ideas lead to powerful methods. Thank you for your attention.

October 16, 2008 – p. 4