Residual Arnoldi Method for solving large eigenvalue problems

(1)

Residual Arnoldi Method

for solving large eigenvalue problems

Che-Rung Lee

roger@umd.edu

(2)

University of Maryland, College Park

Outline

• _Introduction • _Theory • _Experiments • _Conclusion

(3)

Introduction

• _Eigenproblem

• _{Subspace methods}

• _{Krylov subspace and the Arnoldi process} • _{Problems of the Arnoldi process}

• _{Residual Arnoldi method}

(4)

Eigenproblem

• _Let A be a matrix of order n. If a scalar λ and a

nonzero vector x satisfy

Ax = λx,

• λ is an eigenvalue of A, and

• x is the corresponding (right) eigenvector. • (λ, x) is called an eigenpair.

• _{Eigenproblem: find all or some eigenpairs of} A.

(5)

Eigenproblem

• _Let A be a matrix of order n. If a scalar λ and a

nonzero vector x satisfy

Ax = λx,

• λ is an eigenvalue of A, and

• x is the corresponding (right) eigenvector. • (λ, x) is called an eigenpair.

• _{Eigenproblem: find all or some eigenpairs of} A. • _{In this presentation, we assume that} kxk = 1, and

(6)

Subspace methods

• _When A is large (and maybe sparse), subspace

methods are usually used.

• _Steps

1. Generate a subspace.

2. Extract approximations from that subspace. 3. Convergence test.

(7)

Krylov subspace

• _{Given an unit vector} u₁_{. The Krylov subspace}

Kk(A, u1) is a subspace spanned by

u1, Au1, . . . , Ak−1u1.

• _{A Krylov subspace usually contains good}

approximations to the eigenvectors corresponding to the eigenvalues on the edge of spectrum.

(8)

The Arnoldi process

• _{Algorithm that generates orthonormal bases of a}

series of Krylov subspaces.

(9)

The Arnoldi process

• _Steps

1. U1 = u1

2. For i=1, 2, . . .

3. Compute v = Au_i

4. Orthogonalization: u_i+1 = (I − U_iU_i∗)v 5. Normalization : u_i+1 = u_i+1/ku_i+1k

(10)

The Arnoldi process

• _Steps

1. U1 = u1

2. For i=1, 2, . . .

3. Compute v = Au_i

4. Orthogonalization: u_i+1 = (I − U_iU_i∗)v 5. Normalization : u_i+1 = u_i+1/ku_i+1k

6. Expand U by u : U_i+1 = (U_i u_i+1)

• _{Arnoldi relation:} AU_k = U_kH_k + β_ku_k+1e∗

k.

(11)

Problem of the Arnoldi process

• _{When there are errors in the computation of} Au,

convergence will stagnate.

(12)

Problem of the Arnoldi process

• _{Example: Let} A be a 100 × 100 nonsymmetric

matrix with eigenvalues 1, 0.95, . . . , 0.9599.

40 40

(13)

Problem of the Arnoldi process

0 10 20 30 40 10−15 10−10 10−5 100 40

(14)

Problem of the Arnoldi process

0 10 20 30 40 10−15 10−10 10−5 100 10 20 30 40 10−15 10−10 10−5 100

x-axis: dimension of subspace. y-axis: error. Add relative error ǫ = 10−3 to Au.

(15)

Residual Arnoldi method

• _{Subspace is expanded by residuals.}

• _Let (µ, z) be an eigenpair approximation. Its

residual is r = Az − µz.

• _{The pair} (µ, z) is called a candidate.

(16)

Residual Arnoldi method

0 10 20 30 40 10−15 10−10 10−5 100 40

(17)

Residual Arnoldi method

0 10 20 30 40 10−15 10−10 10−5 100 0 10 20 30 40 10−15 10−10 10−5 100 3

(18)

Change candidate

• _{What happens if we select a different candidate}

during the computation?

40

(19)

Change candidate

• _{What happens if we select a different candidate}

during the computation?

• _{Example: candidate is changed at iteration 30.}

0 10 20 30 40

10−15 10−10 10−5 100

(20)

Shift-invert enhancement

• _{Given a shift} σ, the subspace is generated by the

following steps:

1. Select a candidate and compute its residual r. 2. Solve the linear system (A − σI)v = r.

3. Use v is subspace expansion

40

(21)

Shift-invert enhancement

• _{Given a shift} σ, the subspace is generated by the

following steps:

1. Select a candidate and compute its residual r. 2. Solve the linear system (A − σI)v = r.

3. Use v is subspace expansion

• _Example:

Use σ = 1.3 and

solve linear systems to precision 10−3. 0 10 20 30 40 10−15 10−10 10−5 100

(22)

Theory

• _{Perturbation theory of eigenproblem} • _{Algorithm of residual Arnoldi method} • _{Residual Arnoldi relation}

• _{The backward error} • _{Convergence theory}

• _{Shift-invert enhancement}

(23)

Perturbation theory

• _{Let ˜}A = A + E. For an eigenpair (λ, x) of A,

there exists an eigenpair (˜λ, ˜x) of ˜A, such that when kEk is small enough,

˜

λ ≃ λ + (y∗Ex) ˜

x ≃ x + X(λI − L)−1Y ∗Ex,

where (x X) is a nonsingular matrix, (y Y )∗ is its inverse, and L = Y ∗AX.

(24)

Condition number

˜ λ ≃ λ + (y∗Ex) ˜ x ≃ x + X(λI − L)−1Y ∗Ex,

(25)

Condition number

˜

λ ≃ λ + (y∗Ex) ˜

x ≃ x + X(λI − L)−1Y ∗Ex,

• _cond(λ) = kykkxk, and

(26)

Condition number

˜ λ ≃ λ + (y∗Ex) ˜ x ≃ x + X(λI − L)−1Y ∗Ex,

cond(x) = sep−1(λ, L) = k(λI − L)−1k.

• _Therefore, |˜λ − λ| ≤ cond(λ)kEk

k˜x − xk ≤ C1cond(x)kEk.

(27)

Condition number

˜

λ ≃ λ + (y∗Ex) ˜

x ≃ x + X(λI − L)−1Y ∗Ex,

cond(x) = sep−1(λ, L) = k(λI − L)−1k.

• _Therefore, |˜λ − λ| ≤ cond(λ)kEk

k˜x − xk ≤ C1cond(x)kEk.

• _If E has some special structure, kExk ≪ kEk,

|˜λ − λ| ≤ cond(λ)kExk

(28)

Start from the algorithm

1. Compute an eigenpair approximation

2. Compute its residual (inexactly).

3. Orthogonalize residual against U_k

(29)

Start from the algorithm

• _{Using Rayleigh–Ritz method.}

• _{Rayleigh quotient:} H_k = U∗

k AUk.

• _{Ritz pair:} (µ_k, U_ky_k), where H_ky_k = µ_ky_k_.

(30)

Start from the algorithm

k AUk.

• r˜_k = r_k + f_k = AU_ky_k − µ_kU_ky_k + f_k_.

• _{Relative error condition:} kf_kk ≤ ǫkr_kk.

(31)

Start from the algorithm

k AUk.

• r˜_k = r_k + f_k = AU_ky_k − µ_kU_ky_k + f_k_.

• _{Relative error condition:} kf_kk ≤ ǫkr_kk.

• r_k + f⊥

k = Ukgk + βkuk+1e∗k

• f⊥

(32)

Residual Arnoldi relation

For i = 1, · · · , k

• _{Put all} y_i _{into an upper triangular matrix} Y_k • _{Put all} µ_i _{into a diagonal matrix} M_k

• _{Put all} g_i _and β_i _(except β_k_{) into an upper}

Hessenberg matrix G_k

• _{Put all} f⊥

k into Fk.

(33)

Residual Arnoldi relation

For i = 1, · · · , k

Hessenberg matrix G_k • _{Put all} f⊥ k into Fk. AU_k + F_kY_k−1 = U_k(G_k + Y_kM_k)Y_k−1 + βk ηk u_k+1eT_k

(34)

Residual Arnoldi relation

For i = 1, · · · , k

Hessenberg matrix G_k • _{Put all} f⊥ k into Fk. AU_k + F_kY_k−1 = U_k(G_k + Y_kM_k)Y_k−1 + βk ηk u_k+1eT_k • (G_k + Y_kM_k)Y −1 k is upper Hessenberg.

(35)

Backward error

(36)

Backward error

AUk + FkY_k−1 = Uk(Gk + YkMk)Y_k−1 + β_η_kk uk+1eT_k • _Let E_k = F_kY −1 k Uk∗. (A + E_k)U_k = U_k(G_k + Y_kM_k)Y_k−1 + βk ηk uk+1e T k .

(37)

Backward error

• _{By the uniqueness of Arnoldi relation,} E_k _{is the}

backward error of the residual Arnoldi method, and U_k spans a Krylov subspace of A + E_k.

(38)

Backward error

• _{By the uniqueness of Arnoldi relation,} E_k _{is the}

backward error of the residual Arnoldi method, and U_k spans a Krylov subspace of A + E_k.

• _Empirically, kE_kk is around the level of ǫ. With

that, we can prove kE_kxk ≤ C2ǫkr_kk.

(39)

Some notation

Let (˜λ_k, ˜x_k) be an eigenpair of A + E_k corresponding to (λ, x) of A, and let (˜µ_k, ˜z_k) be the candidate

(40)

Convergence theory

Some assumptions for our proof.

(41)

Convergence theory

(42)

Convergence theory

1. The target eigenpair (λ, x) is simple.

2. There exists a constant C3 > 0 such that

sep−1(µ_k, L) < C3.

• _Matrix L = X∗AX, where X is an

orthonormal basis of span{I − xx∗}.

(43)

Convergence theory

sep−1(µ_k, L) < C3.

3. There exists a positive constant _C4 such that if

kEk ≤ C4, then there are descending constants

˜

κ1, ˜κ2, . . . independent of E, with lim_kκ˜_k = 0

(44)

Convergence theory

sep−1(µ_k, L) < C3.

3. There exists a positive constant _C4 such that if

kEk ≤ C4, then there are descending constants

˜

κ1, ˜κ2, . . . independent of E, with lim_kκ˜_k = 0

such that k˜z_k − ˜x_kk ≤ ˜κ_k. 4. kE_kk ≤ ǫC5.

(45)

Put everything together

• _{Perturbation theory:} k˜x_k − xk ≤ C₆kE_kxk. • _{Backward error:} kE_kxk ≤ C₂ǫkr_kk. • _{Assumption 3:} k˜z_k − ˜x_kk ≤ ˜κ_k_. • _{Residual bound:} kr_kk ≤ 2kz_k − xk. • _{Invariant property} z˜_k = z_k_.

(46)

Put everything together

• _{Perturbation theory:} k˜x_k − xk ≤ C₆kE_kxk. • _{Backward error:} kE_kxk ≤ C₂ǫkr_kk. • _{Assumption 3:} k˜z_k − ˜x_kk ≤ ˜κ_k_. • _{Residual bound:} kr_kk ≤ 2kz_k − xk. • _{Invariant property} z˜_k = z_k_. If ǫ ≤ 1/2C2C6, kr_kk ≤ 2˜κk 1 − 2C2C6ǫ .

(47)

Shift-invert enhancement

• _Let S = (A − σI)−1 _and T_k = (σI − M_k)−1Y −1

k .

The SIRA relation is

SUk+FkTk = Uk(Hk−Yk)Tk+

βk

(σ − µk)ηk

(48)

Shift-invert enhancement

k .

SUk+FkTk = Uk(Hk−Yk)Tk+ βk (σ − µk)ηk uk+1eT_k . • _{Backward error} E_k = F_kT_kU∗ k .

(49)

Shift-invert enhancement

k .

(50)

Shift-invert enhancement

k .

• _{We can prove that} kE_kxk ≤ C₈ǫkr_kk. • _If ǫ < 1/C₉ _{for some constant} C₉_{, then}

krkk ≤

2˜κ_k 1 − C9ǫ

.

(51)

Experiments

• _RAPACK

• _{Compare with ARPACK} • _{Compare with SRRIT} • _{Inexact Krylov method}

(52)

RAPACK

• _{A numerical package implementing the residual}

Arnoldi method.

(53)

RAPACK

Arnoldi method.

(54)

RAPACK

Arnoldi method.

• _{Two computational modes: RA and SIRA.} • _{Uses reverse communication to get matrix}

operation results.

(55)

RAPACK

Arnoldi method.

operation results.

• _{Implements Krylov Schur restarting method for}

(56)

RAPACK

Arnoldi method.

operation results.

• _{Implements Krylov Schur restarting method for}

memory management.

• _{Allows an arbitrary initial subspace.}

(57)

ARPACK

(58)

ARPACK

• _{Implements implicitly restarted Arnoldi method.} • _{Can solve standard eigenproblems, generalized}

eigenproblems, and singular value

decompositions for symmetric, nonsymmetric and complex matrices.

(59)

ARPACK

• _{Four computational modes.}

• _{Mode 1: Standard eigenproblem.}

• _{Mode 2: Generalized eigenproblem.}

(60)

ARPACK

• _{Four computational modes.}

• _{Mode 1: Standard eigenproblem.}

• _{Mode 2: Generalized eigenproblem.}

• _{Mode 3, 4: With shift-invert enhancement.} • _{We only compare mode 1 with the RA mode and}

mode 3 with the SIRA mode.

(61)

Test problem

• _{A real nonsymmetric eigenmat} A of order 10000. • _{First 100 eigenvalues are} 1, 0.95, · · · , 0.9599_. • _{Other eigenvalues are in} (0.25, 0.75).

(62)

Test problem

• _{The condition number is around} 105_. • _Tasks

• _{Compute 6 largest eigenvalues using mode 1}

of ARPACK and the RA mode of RAPACK.

• _{Compute 6 smallest eigenvalues using mode 3}

of ARPACK and the SIRA mode of RAPACK.

(63)

Test problem

• _{The condition number is around} 105_. • _Tasks

• _{Compute 6 largest eigenvalues using mode 1}

of ARPACK and the RA mode of RAPACK.

• _{Compute 6 smallest eigenvalues using mode 3}

of ARPACK and the SIRA mode of RAPACK.

• _Setting

• _{Maximum dimension of subspace 20.} • _{Convergence precision} 10−13_.

(64)

Mode 1 and the RA mode

Mode 1 RA mode ETime (second) 4.6860 8.4242

MVM 113 138

Let x_i be the ith eigen-vector of A. The error is measured by kxi − U U∗xik. 0 50 100 150 10−15 10−10 10−5 100

(65)

Mode 3 and the SIRA mode

Use GMRES to solve linear systems with shift=0.

Mode 3 SIRA mode

Etime (second) 378 168

MVM 11842 4606

Outer iterations 68 144 Precision for solving ₁₀−13 ₁₀−3

(66)

Mode 3 and the SIRA mode

Mode 3 SIRA mode

MVM 11842 4606

20 40 60 80 100 120 0 50 100 150 200 250 150

(67)

Mode 3 and the SIRA mode

Mode 3 SIRA mode

MVM 11842 4606

50 100 150 200 250 10−10 10−5 100

(68)

SRRIT

• _{Implements Schur–Rayleigh–Ritz iteration}

method.

• _{Can compute the dominant invariant subspace.}

• _{Can use an arbitrary subspace to start the process.}

(69)

SRRIT

• _{Implements Schur–Rayleigh–Ritz iteration}

method.

• _{Can compute the dominant invariant subspace.}

• _{Can use an arbitrary subspace to start the process.} • _{Compare it with the RA mode for using an}

existing subspace as initialization.

• _{Use matrix} S = A−1_{, where} A is previously

defined, to compute 6 smallest eigenvalues.

(70)

Successive inner-outer process

Properties of Krylov subspaces 1. Stagnate around error level 2. Invariant convergent curves 3. Superlinear convergence 0 10 20 30 40 10−15 10−10 10−5 100 1.e−3 1.e−6 1.e−9 1.e−12

(71)

Successive inner-outer process

Properties of Krylov subspaces 1. Stagnate around error level 2. Invariant convergent curves 3. Superlinear convergence 0 10 20 30 40 10−15 10−10 10−5 100 1.e−3 1.e−6 1.e−9 1.e−12 Algorithm:

1. Divide the process into 4 stages with increasing precision requirement 10−3, 10−6, 10−9, 10−12. 2. Stage i computes matrix-vector multiplication

(solves linear system) to precision ǫ_i.

3. Each stage uses previously generated subspace as an initial subspace.

(72)

SRRIT and the RA mode

Stage 10−3 10−6 10−9 10−12 Total SRRIT 305 120 152 213 783 RA mode 39 40 40 34 153 50 100 150 10−10 10−5 100

(73)

Inexact Krylov method

• _{Allows increasing errors in matrix-vector}

(74)

Inexact Krylov method

multiplication.

• _{Implemented in the RA mode with} S = A−1_{, and}

tolerable error size max ǫ, _mkrǫτ

i −1k

ǫ: relative error

τ : convergence precision

m: maximum number of subspace ri: residual in the ith iteration

(75)

Inexact Krylov method

multiplication.

i −1k

ǫ: relative error (=10−3) τ : convergence precision (=10−12) m: maximum number of subspace (=50)

(76)

Inexact Krylov method

multiplication.

i −1k

ǫ: relative error (=10−3) τ : convergence precision (=10−12) m: maximum number of subspace (=50)

ri: residual in the ith iteration

• _{Compare it with mode 3 and the SIRA mode to}

compute 6 smallest eigenvalues of A.

(77)

Result of inexact Krylov method

Inexact Mode 3 SIRA

Etime 80 106 48

TMVM 5240 7083 2829 Iteration 43 50 89

(78)

Result of inexact Krylov method

Inexact Mode 3 SIRA

Etime 80 106 48 TMVM 5240 7083 2829 Iteration 43 50 89 20 40 60 80 0 50 100 150 Inexact Krylov Mode 3 SIRA Mode 0 10 20 30 40 50 10−15 10−10 10−5 100

(79)

Conclusion

• _{Residual Arnoldi method for eigenproblems}

allows errors in the computation, and can work on an appropriate initial subspace.

(80)

Conclusion

• _{With shift-invert enhancement, residual Arnoldi}

method can reduce a lot of computational cost.

(81)

Conclusion

• _{RAPACK can compute few selected eigenpairs}

for real matrices efficiently, and only requires moderate memory.

(82)

Conclusion

• _{RAPACK can compute few selected eigenpairs}

for real matrices efficiently, and only requires moderate memory.

• _{Many other algorithms can be implemented in}

RAPACK and get better performance.

(83)

Future work

• _{Block residual Arnoldi method.}

• _{Using other eigenvector approximations, such as}

refine Ritz vector or harmonic Ritz vector.

• _{More inexact Krylov subspace methods.} • _{Extension of RAPACK to solve other}

eigenproblems, such as generalized eigenproblem.

• _{Parallelization of RAPACK.}