Krylov sequence methods

(1)

師大

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

March 22, 2009

(2)

師大

Outline

1 Krylov decompositions

Householder transformation Arnoldi decompositions Lanczos decompositions Krylov decompositions

Computation of refined and Harmonic Ritz vectors

2 Restarted Arnoldi method

The implicitly restarted Arnoldi method Krylov-Schur restarting

3 The Lanczos Algorithm

(3)

師大 Householder transformation

Definition

A Householder transformation or elementary reflector is a matrix of

H = I − uu^∗ where kuk2 =√

2.

Note that H is Hermitian and unitary.

Theorem

Let x be a vector such that kxk₂ = 1and x₁ is real and nonnegative. Let

u = (x + e1)/√ 1 + x1. Then

Hx = (I − uu^∗)x = −e1.

(4)

Proof:

I − uu^∗x = x − (u^∗x)u = x − x^∗x + x₁

√1 + x₁ · x + e₁

√1 + x₁

= x − (x + e1) = −e1

Theorem

Let x be a vector with x₁6= 0. Let u = ρ_kxk^x

2 + e₁ q1 + ρ_kxk^x¹

2

,

where ρ = ¯x₁/|x₁|. Then

Hx = − ¯ρkxk2e1.

(5)

Proof: Since

[ ¯ρx^∗/kxk₂+ e^T₁][ρx/kxk₂+ e₁]

= ρρ + ρx¯ ₁/kxk₂+ ¯ρ¯x₁/kxk₂+ 1

= 2[1 + ρx1/kxk2], it follows that

u^∗u = 2 ⇒ kuk₂=√ 2 and

u^∗x = ρkxk¯ 2+ x1

q1 + ρ_kxk^x¹

2

.

(6)

Hence,

Hx = x − (u^∗x)u = x − ρkxk¯ 2+ x1

q1 + ρ_kxk^x¹

2

ρ_kxk^x

2 + e1

q1 + ρ_kxk^x¹

2

=

"

1 −( ¯ρkxk2+ x1)_kxk^ρ

2

1 + ρ_kxk^x¹

2

#

x −ρkxk¯ 2+ x1

1 + ρ_kxk^x¹

2

e1

= −ρkxk¯ 2+ x1

1 + ρ_kxk^x¹

2

e1

= −¯ρkxk2e1.

(7)

Definition

A complex m × n-matrix R = [rij]is called an upper (lower) triangular matrix, if r_ij = 0for i > j (i < j).

Definition

Given A ∈ C^m×n, Q ∈ C^m×m unitary and R ∈ C^m×nupper triangular such that A = QR. Then the product is called a QR-factorization of A.

Theorem

Any complex m × n matrix A can be factorized by the product A = QR, where Q is m × m-unitary. R is m × n upper triangular.

(8)

Proof: Let A⁽⁰⁾ = A = [a⁽⁰⁾₁ |a⁽⁰⁾₂ | · · · |a⁽⁰⁾_n ]. Find Q₁ = (I − 2w₁w^∗₁)such that Q₁a⁽⁰⁾₁ = ce₁. Then

A⁽¹⁾ = Q₁A⁽⁰⁾= [Q₁a⁽⁰⁾₁ , Q₁a⁽⁰⁾₂ , · · · , Q₁a⁽⁰⁾_n ]

=







c₁ ∗ · · · ∗ 0

... a⁽¹⁾₂ · · · a⁽¹⁾n

0







. (1)

Find Q₂=

1 0

0 I − w₂w^∗₂

such that (I − 2w₂w^∗₂)a⁽¹⁾₂ = c₂e₁. Then

A⁽²⁾= Q2A⁽¹⁾=







c₁ ∗ ∗ · · · ∗ 0 c2 ∗ · · · ∗ 0 0

... ... a⁽²⁾₃ · · · a⁽²⁾_n 0 0





 .

(9)

We continue this process. Then after l = min(m, n) steps A^(l) is an upper triangular matrix satisfying

A^(l−1)= R = Ql−1· · · Q₁A.

Then A = QR, where Q = Q^∗₁· · · Q^∗_l−1. Theorem

Let A be a nonsingular n × n matrix. Then the QR- factorization is essentially unique. That is, if A = Q₁R1 = Q2R2, then there is a unitary diagonal matrix D = diag(di)with |di| = 1 such that Q₁ = Q₂Dand DR₁ = R₂.

Proof: Let A = Q1R1 = Q2R2. Then Q^∗₂Q1= R2R⁻¹₁ = Dmust be a diagonal unitary matrix.

(10)

師大 Arnoldi decompositions

Suppose that the columns of K_k+1 are linearly independent and let

K_k+1= U_k+1R_k+1

be the QR factorization of K_k+1. Then the columns of U_k+1 are results of successively orthogonalizing the columns of K_k+1. Theorem

Let ku₁k₂= 1and the columns of K_k+1(A, u₁)be linearly independent. Let U_k+1 = [ u₁ · · · u_k+1 ]be the Q-factor of Kk+1. Then there is a (k + 1) × k unreduced upper Hessenberg matrix ˆH_k such that

AU_k= U_k+1Hˆ_k. (2)

Conversely, if U_k+1is orthonormal and satisfies (2), where ˆHk

is a (k + 1) × k unreduced upper Hessenberg matrix, then U_k+1 is the Q-factor of K_k+1(A, u₁).

(11)

Proof: (“⇒”) Let K_k= U_kR_k be the QR factorization and Sk= R_k⁻¹. Then

AU_k= AK_kS_k = K_k+1

0 S_k

= U_k+1R_k+1

0 S_k

= U_k+1Hˆ_k, where

Hˆ_k= R_k+1

0 Sk

.

It implies that ˆH_kis a (k + 1) × k Hessenberg matrix and h_i+1,i = r_i+1,i+1s_ii= ri+1,i+1

r_ii . Thus by the nonsingularity of R_k, ˆH_kis unreduced.

(“⇐”) If k = 1, then

Au1 = h11u1+ h21u2 ⇒ u2 = −h₁₁

h₂₁ u1+ 1 h₂₁Au1.

(12)

Since [ u₁ u₂ ]is orthonormal and u₂ is a linear combination of u₁ and Au₁, [ u₁ u₂ ]is the Q-factor of K₂.

Assume U_k is the Q-factor of K_k. If we partition Hˆ_k =

Hˆk−1 hk

0 h_k+1,k

, then from (2)

Auk = Ukhk+ hk+1,kuk+1.

Thus u_k+1is a linear combination of Au_k and the columns of Uk. Hence U_k+1is the Q-factor of K_k+1.

(13)

Definition

Let U_k+1∈ C^n×(k+1)be orthonormal. If there is a (k + 1) × k unreduced upper Hessenberg matrix ˆHk such that

AU_k = U_k+1Hˆ_k, (3)

then (3) is called an Arnoldi decomposition of order k. If ˆH_kis reduced, we say the Arnoldi decomposition is reduced.

Partition

Hˆk=

Hk

h_k+1,ke^T_k

, and set

βk= hk+1,k. Then (3) is equivalent to

AU_k= U_kH_k+ β_ku_k+1e^T_k.

(14)

Theorem

Suppose the Krylov sequence K_k+1(A, u1)does not terminate at k + 1. Then up to scaling of the columns of U_k+1, the Arnoldi decomposition of K_k+1is unique.

Proof: Since the Krylov sequence K_k+1(A, u₁)does not terminate at k + 1, the columns of K_k+1(A, u₁)are linearly independent. By Theorem 8, there is an unreduced matrix H_k and β_k6= 0 such that

AU_k= U_kH_k+ β_ku_k+1e^T_k, (4) where U_k+1= [Ukuk+1]is an orthonormal basis for

K_k+1(A, u1). Suppose there is another orthonormal basis U˜_k+1= [ ˜U_ku˜_k+1]for K_k+1(A, u₁), unreduced matrix ˜H_kand β˜k6= 0 such that

A ˜U_k= ˜U_kH˜_k+ ˜β_ku˜_k+1e^T_k.

(15)

Then we claim that

U˜_k^Hu_k+1 = 0.

For otherwise there is a column ˜u_j of ˜U_ksuch that

˜

u_j = αu_k+1+ U_ka, α 6= 0.

Hence

A˜u_j = αAu_k+1+ AU_ka

which implies that A˜u_j contains a component along A^k+1u₁. Since the Krylov sequence K_k+1(A, u1)does not terminate at k + 1, we have

K_k+2(A, u1) 6= K_k+1(A, u1).

Therefore, A˜uj lies in K_k+2(A, u1)but not in K_k+1(A, u1)which is a contradiction.

Since U_k+1 and ˜Uk+1are orthonormal bases for K_k+1(A, u1) and ˜U_k^Hu_k+1= 0, it follows that

(16)

R(U_k) = R( ˜Uk) and U_k^Hu˜k+1 = 0, that is

U_k= ˜U_kQ for some unitary matrix Q. Hence

A( ˜UkQ) = ( ˜UkQ)(Q^HH˜kQ) + ˜βku˜k+1(e^T_kQ), or

AU_k= U_k(Q^HH˜_kQ) + ˜β_ku˜_k+1e^T_kQ. (5) On premultiplying (4) and (5) by U_k^H, we obtain

H_k = U_k^HAU_k = Q^HH˜_kQ.

Similarly, premultiplying by u^H_k+1, we obtain βke^T_k = u^H_k+1AUk = ˜βk(u^H_k+1u˜k+1)e^T_kQ.

(17)

It follows that the last row of Q is ω_ke^T_k, where |ω_k| = 1. Since the norm of the last column of Q is one, the last column of Q is ω_ke_k. Since H_kis unreduced, it follows from the implicit Q theorem that

Q =diag(ω₁, · · · , ωk), |ω_j| = 1, j = 1, . . . , k.

Thus up to column scaling U_k = ˜U_kQis the same as ˜U_k. Subtracting (5) from (4), we find that

β_ku_k+1 = ω_kβ˜_ku˜_k+1

so that up to scaling u_k+1and ˜u_k+1are the same.

(18)

Theorem

Let the orthonormal matrix U_k+1satisfy AUk= Uk+1Hˆk,

where ˆH_kis Hessenberg. Then ˆH_kis reduced if and only if R(U_k)contains an eigenspace of A.

Proof: (“⇒”) Suppose that ˆH_kis reduced, say that h_j+1,j= 0.

Partition Hˆ_k=

H₁₁ H₁₂ 0 H₂₂

and U_k= [ U₁₁ U₁₂ ], where H₁₁is an j × j matrix and U₁₁is consisted the first j columns of U_k+1. Then

A[ U₁₁ U₁₂ ] = [ U₁₁ U₁₂ u_k+1 ]

H11 H12

0 H₂₂

.

(19)

It implies that

AU11= U11H11

so that U₁₁is an eigenbasis of A.

(“⇐”) Suppose that A has an eigenspace that is a subset of R(U_k)and ˆH_kis unreduced. Let (λ, U_kw)for some w be an eigenpair of A. Then

0 = (A − λI)Ukw = (Uk+1Hˆk− λU_k)w

=

U_k+1Hˆ_k− λU_k+1

I 0

w = U_k+1Hˆ_λw, where

Hˆ_λ=

Hk− λI h_k+1,ke^T_k

.

Since ˆH_λis unreduced, the matrix U_k+1Hˆ_λ is of full column rank. It follows that w = 0 which is a contradiction.

(20)

Write the k-th column of the Arnoldi decomposition AU_k= U_kH_k+ β_ku_k+1e^T_k,

in the form

Auk= Ukhk+ βkuk+1. Then from the orthonormality of U_k+1, we have

hk= U_k^HAuk. Since

βkuk+1= Auk− U_khk

and ku_k+1k₂ = 1, we must have

βk= kAuk− U_khkk₂ and

u_k+1 = β_k⁻¹(Au_k− U_kh_k).

(21)

Algorithm (Arnoldi process) 1. for k = 1, 2, . . .

2. h_k = U_k^HAu_k 3. v = Auk− U_khk

4. β_k= h_k+1,k = kvk2

5. u_k+1 = v/β_k 6. Hˆ_k=

Hˆ_k−1 h_k 0 h_k+1,k

7. end for k

The computation of u_k+1 is actually a form of the well-known Gram-Schmidt algorithm.

In the presence of inexact arithmetic cancelation in statement 3 can cause it to fail to produce orthogonal vectors.

The cure is process called reorthogonalization.

(22)

Algorithm (Reorthogonalized Arnoldi process) for k = 1, 2, . . .

hk= U_k^HAuk

v = Au_k− U_kh_k w = U_k^Hv hk= hk+ w v = v − U_kw β_k= h_k+1,k= kvk₂ uk+1= v/βk

Hˆ_k=

Hˆ_k−1 h_k 0 hk+1,k

end for k

(23)

師大 Lanczos decompositions

Let A be Hermitian and let

AU_k= U_kT_k+ β_ku_k+1e^T_k (6) be an Arnoldi decomposition. Since T_kis upper Hessenberg and T_k= U_k^HAU_kis Hermitian, it follows that T_kis tridiagonal and can be written in the form

Tk =







α1 β¯1

β1 α2 β¯2

β₂ α₃ β¯₃ . .. ... . ..

β_k−2 α_k−1 β¯_k−1 βk−1 αk





 .

Equation (6) is called a Lanczos decomposition. The first column of (6) is

Au₁= α₁u₁+ β₁u₂,

(24)

or

u2 = Au1− α₁u1

β₁ .

From the orthonormality of u₁and u₂, it follows that α₁= u^H₁ Au₁

and

β1 = kAu1− α₁u1k₂.

More generality, from the j-th column of (6) we get the relation u_j+1= Auj− α_juj− ¯βj−1uj−1

β_j where

αj = u^H_j Auj and βj = kAuj− α_juj− ¯βj−1uj−1k₂. This is the Lanczos three-term recurrence.

(25)

Algorithm (Lanczos recurrence)

Let u₁ be given. This algorithm generates the Lanczos decomposition

AU_k= U_kT_k+ β_ku_k+1e^T_k where T_kis Hermitian tridiagonal.

1. u0 = 0; β0 = 0;

2. for j = 1 to k 3. u_j+1 = Au_j 4. αj = u^H_j uj+1

5. v = uj+1− α_juj− β_j−1uj−1

6. β_j = kvk₂ 7. uj+1 = v/βj

8. end for j

(26)

師大 Krylov decompositions

Definition

Let u₁, u₂, . . . , u_k+1be linearly independent and let U_k= [u₁ · · · u_k].

AU_k= U_kB_k+ u_k+1b^H_k+1

is called a Krylov decomposition of order k. R(U_k+1)is called the space spanned by the decomposition. Two Krylov

decompositions spanning the same spaces are said to be equivalent.

Let [V v]^H be any left inverse for U_k+1. Then it follows that Bk = V^HAUk and b^H_k+1 = v^HAUk.

In particular, B_k is a Rayleigh quotient of A.

(27)

Let

AU_k= U_kB_k+ u_k+1b^H_k+1

be a Krylov decomposition and Q be nonsingular. That is AU_k= U_k+1Bˆ_k with Bˆ_k=

Bk

b^H_k+1

. (7)

Then we get an equivalent Krylov decomposition of (7) in the form

A(U_kQ) =

U_k+1

Q 0 0 1

Q 0

0 1

−1

Bˆ_kQ

!

=

U_kQ u_k+1

Q⁻¹BkQ b^H_k+1Q

= (U_kQ)(Q⁻¹BQ) + u_k+1(b^H_k+1Q). (8) The two Krylov decompositions (7) and (8) are said to be similar.

(28)

Let

γ ˜u_k+1= u_k+1− U_ka.

Since u₁, . . . , u_k, u_k+1are linearly independent, we have γ 6= 0.

Then it follows that

AU_k= U_k(B_k+ ab^H_k+1) + ˜u_k+1(γb^H_k+1).

Since R([U_ku_k+1]) = R([U_ku˜_k+1]), this Krylov decomposition is equivalent to (7).

Theorem

Every Krylov decomposition is equivalent to a (possibly reduced) Arnoldi decomposition.

Proof: Let

AU = U B + ub^H be a Krylov decomposition and let

(29)

U = ˜U R be the QR factorization of U . Then

A ˜U = A(U R⁻¹) = (U R⁻¹)(RBR⁻¹) + u(b^HR⁻¹) ≡ ˜U ˜B + u˜b^H is an equivalent decomposition. Let

˜

u = γ⁻¹(u − U a)

be a vector with k˜uk₂ = 1such that U^Hu = 0. Then˜ A ˜U = ˜U ( ˜B + a˜b^H) + ˜u(γ˜b^H) ≡ ˜U ˆB + ˜uˆb^H

is an equivalent orthonormal Krylov decomposition. Let Q be a unitary matrix such that

ˆb^HQ = kˆbk2e^T_k

and Q^HBQˆ is upper Hessenberg. Then the equivalent decomposition

(30)

A ˆU ≡ A( ˜U Q) = ( ˜U Q)(Q^HBQ) + ˜ˆ u(ˆb^HQ) ≡ ˆU ¯B + kˆbk₂ueˆ ^T_k is a possibly reduced Arnoldi decomposition where

Uˆ^Hu = Qˆ ^HU˜^Hu = Q˜ ^HR^−HU^Hu = 0.˜

Reduction to Arnoldi form Let

AU = U B + ub^H

be the Krylov decomposition with B ∈ C^k×k. Let H1 be a Householder transformation such that

b^HH1 = βe_k.

(31)

Reduce H₁^HBH₁ to Hessenberg form as the following illustration:

B :=







× × × ×







⇒B := BH₂=







⊗ ⊗ ⊗ ×

⊗ ⊗ ⊗ × 0 0 ⊗ ×







⇒ B := H₂^HB =







+ + + + + + + + + + + + 0 0 ⊗ ×







⇒B := BH3=







⊕ ⊕ + +

⊕ ⊕ + + 0 ⊕ + + 0 0 ⊗ ×







⇒ B := H₃^HB =







∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ 0 ⊕ + + 0 0 ⊗ ×







(32)

Let

Q = H1H2· · · H_k−1. Then Q^HBQis upper Hessenberg and

b^HQ = (b^HH₁)(H₂· · · H_k−1) = βe^T_k(H₂· · · H_k−1) = βe^T_k. Therefore, the Krylov decomposition

A(U Q) = (U Q)(Q^HBQ) + βue^T_k (9) is an Arnoldi decomposition.

(33)

師大 Computation of refined and Harmonic Ritz vectors

Assume that

AU = U B + ub^H is a n orthonormal Krylov decomposition.

Refined Ritz vectors

If µ is a Ritz value, then the refined Ritz vector associated with µis the right singular vector of (A − µI)U whose singular value is smallest. From (9), we have

(A − µI)U = U (B − µI) + ub^H =

U u

B − µI b^H

≡

U u Bˆµ.

Since [U u] is orthonormal, the right singular vectors of (A − µI)Uare the same as the right singular vectors of ˆB_µ. Thus the computation of a refined Ritz vector can be reduced to computing the singular value decomposition of ˆB_µ.

(34)

師大 Computation of refined and Harmonic Ritz vectors

Harmonic Ritz vectors

Recall: (κ + δ, U w) is a harmonic Ritz pair if

U^H(A − κI)^H(A − κI)U w = δU^H(A − κI)^HU w.

Since

(A − κI)U = U (B − κI) + ub^H, we have

U^H(A − κI)^H(A − κI)U = (B − κI)^H(B − κI) + bb^H and

U^H(A − κI)^HU = (B − κI)^H. It follows that

(B − κI)^H(B − κI) + bb^H w = δ(B − κI)^Hw which is a small generalized eigenvalue problem.

(35)

師大

Let

AU_k= U_kH_k+ β_ku_k+1e^T_k be an Arnoldi decomposition.

1 In principle, we can keep expanding the Arnoldi decomposition until the Ritz pairs have converged.

2 Unfortunately, it is limited by the amount of memory to storage of U_k.

3 Restarted the Arnoldi process once k becomes so large that we cannot store U_k.

Implicitly restarting method Krylov-Schur decomposition

(36)

師大 The implicitly restarted Arnoldi method

Choose a new starting vector for the underlying Krylov sequence

A natural choice would be a linear combination of Ritz vectors that we are interested in.

Filter polynomials

Assume A has a complete system of eigenpairs (λi, xi)and we are interested in the first k of these eigenpairs. Expand u₁in the form

u1 =

k

X

i=1

γixi+

n

X

i=k+1

γixi. If p is any polynomial, we have

p(A)u1 =

k

X

i=1

γip(λi)xi+

n

X

i=k+1

γip(λi)xi.

(37)

Choose p so that the values p(λ_i) (i = k + 1, . . . , n)are small compared to the values p(λ_i) (i = 1, . . . , k).

Then p(A)u1is rich in the components of the xi that we want and deficient in the ones that we do not want.

pis called a filter polynomial.

Suppose we have Ritz values µ₁, . . . , µ_m and µ_k+1, . . . , µ_m are not interesting. Then take

p(t) = (t − µ_k+1) · · · (t − µ_m).

Implicitly restarted Arnoldi: Let

AU_m = U_mH_m+ β_mu_m+1e^T_m (10) be an Arnoldi decomposition with order m. Choose a filter polynomial p of degree m − k and use the implicit restarting process to reduce the decomposition to a decomposition

A ˜U_k= ˜U_kH˜_k+ ˜β_ku˜_k+1e^T_k of order k with starting vector p(A)u₁.

(38)

Let κ₁, . . . , κ_m be eigenvalues of H_m and suppose that

κ1, . . . , κm−kcorrespond to the part of the spectrum we are not interested in. Then take

p(t) = (t − κ₁)(t − κ₂) · · · (t − κ_m−k).

The starting vector p(A)u₁ is equal to

p(A)u₁ = (A − κ_m−kI) · · · (A − κ₂I)(A − κ₁I)u₁

= (A − κ_m−kI) [· · · [(A − κ2I) [(A − κ1I)u1]]] . In the first, we construct an Arnoldi decomposition with starting vector (A − κ₁I)u₁. From (10), we have

(A − κ₁I)U_m = U_m(H_m− κ₁I) + β_mu_m+1e^T_m (11)

= U_mQ₁R₁+ β_mu_m+1e^T_m, where

H_m− κ₁I = Q₁R₁

is the QR factorization of H_m− κ₁I. Postmultiplying by Q₁,

(39)

we get

(A − κ₁I)(U_mQ₁) = (U_mQ₁)(R₁Q₁) + β_mu_m+1(e^T_mQ₁).

It implies that

AU_m⁽¹⁾ = U_m⁽¹⁾H_m⁽¹⁾+ βmum+1b^(1)H_m+1, where

U_m⁽¹⁾ = U_mQ₁, H_m⁽¹⁾= R₁Q₁+ κ₁I, b^(1)H_m+1 = e^T_mQ₁. (H_m⁽¹⁾: one step of single shifted QR algorithm)

(40)

Theorem

Let H_m be an unreduced Hessenberg matrix. Then H_m⁽¹⁾has the form

H_m⁽¹⁾=

"

Hˆm⁽¹⁾ hˆ₁₂ 0 κ₁

# ,

where ˆH_m⁽¹⁾is unreduced.

Proof: Let

H_m− κ₁I = Q₁R₁ be the QR factorization of H_m− κ₁I with

Q1= G(1, 2, θ1) · · · G(m − 1, m, θm−1)

where G(i, i + 1, θ_i)for i = 1, . . . , m − 1 are Given rotations.

(41)

Since H_mis unreduced upper Hessenberg, i.e., the subdiagonal elements of H_m are nonzero, we get

θ_i 6= 0 for i = 1, . . . , m − 1 (12) and

(R1)ii6= 0 for i = 1, . . . , m − 1. (13) Since κ₁ is an eigenvalue of H_m, we have that H_m− κ₁I is singular and then

(R1)mm = 0. (14)

Using the results of (12), (13) and (14), we get

H_m⁽¹⁾ = R₁Q₁+ κ₁I = R₁G(1, 2, θ₁) · · · G(m − 1, m, θ_m−1) + κ₁I

=

"

Hˆm⁽¹⁾ ˆh₁₂ 0 κ1

# , where ˆHm⁽¹⁾is unreduced.

(42)

Remark

Um⁽¹⁾ is orthonormal.

Since H_m is upper Hessenberg and Q₁ is the Q-factor of the QR factorization of H_m− κ₁I, it implies that Q₁ and Hm⁽¹⁾ are also upper Hessenberg.

The vector b^(1)H_m+1 = e^T_mQ₁ has the form b^(1)H_m+1 =

h

0 · · · 0 q⁽¹⁾_m−1,m q⁽¹⁾m,m

i

;

i.e., only the last two components of b⁽¹⁾_m+1are nonzero.

(43)

For on postmultiplying (11) by e₁, we get

(A − κ₁I)u₁= (A − κ₁I)(U_me₁) = U_m⁽¹⁾R₁e₁ = r₁₁⁽¹⁾u⁽¹⁾₁ . Since H_m is unreduced, r⁽¹⁾₁₁ is nonzero. Therefore, the first column of Um⁽¹⁾is a multiple of (A − κ₁I)u₁.

By the definition of Hm⁽¹⁾, we get

Q1H_m⁽¹⁾Q^H₁ = Q1(R1Q1+ κ1I)Q^H₁ = Q1R1+ κ1I = Hm. Therefore, κ₁, κ₂, . . . , κ_mare also eigenvalues of Hm⁽¹⁾.

(44)

Similarly,

(A − κ₂I)U_m⁽¹⁾ = U_m⁽¹⁾(H_m⁽¹⁾− κ₂I) + β_mu_m+1b^(1)H_m+1 (15)

= U_m⁽¹⁾Q2R2+ βmum+1b^(1)H_m+1, where

H_m⁽¹⁾− κ₂I = Q₂R₂

is the QR factorization of Hm⁽¹⁾− κ₂I with upper Hessenberg matrix Q₂. Postmultiplying by Q₂, we get

(A − κ2I)(U_m⁽¹⁾Q2) = (U_m⁽¹⁾Q2)(R2Q2) + βmum+1(b^(1)H_m+1Q2).

It implies that

AU_m⁽²⁾ = U_m⁽²⁾H_m⁽²⁾+ β_mu_m+1b^(2)H_m+1, where

U_m⁽²⁾ ≡ U_m⁽¹⁾Q2

is orthonormal,

(45)

H_m⁽²⁾ ≡ R₂Q2+ κ2I =





H_m−2⁽²⁾ ∗ ∗ κ2 ∗ κ₁





is upper Hessenberg with unreduced matrix H_m−2⁽²⁾ and b^(2)H_m+1 ≡ b^(1)H_m+1Q2 = q_m−1,m⁽¹⁾ e^H_m−1Q2+ q_m,m⁽¹⁾ e^T_mQ2

=

0 · · · 0 × × × . For on postmultiplying (15) by e1, we get

(A − κ₂I)u⁽¹⁾₁ = (A − κ₂I)(U_m⁽¹⁾e₁) = U_m⁽²⁾R₂e₁ = r₁₁⁽²⁾u⁽²⁾₁ . Since Hm⁽¹⁾is unreduced, r⁽²⁾₁₁ is nonzero. Therefore, the first column of Um⁽²⁾is a multiple of

(A − κ2I)u⁽¹⁾₁ = 1/r⁽¹⁾₁₁(A − κ2I)(A − κ1I)u1.

(46)

Repeating this process with κ₃, . . . , κ_m−k, the result will be a Krylov decomposition

AU_m^(m−k)= U_m^(m−k)H_m^(m−k)+ βmum+1b^(m−k)H_m+1 with the following properties

1 Um^(m−k)is orthonormal.

2 H_m^(m−k)is upper Hessenberg.

3 The first k − 1 components of b^(m−k)H_m+1 are zero.

4 The first column of Um^(m−k) is a multiple of (A − κ₁I) · · · (A − κ_m−kI)u₁.

(47)

Corollary

Let κ₁, . . . , κ_m be eigenvalues of H_m. If the implicitly restarted QRstep is performed with shifts κ₁, . . . , κ_m−k, then the matrix Hm^(m−k) has the form

H_m^(m−k)=

"

H_kk^(m−k) H_k,m−k^(m−k) 0 T^(m−k)

# ,

where T^(m−k) is an upper triangular matrix with Ritz value κ₁, . . . , κ_m−kon its diagonal.

(48)

For k = 3 and m = 6, A

u u u u u u

=

u u u u u u







× × × × × ×

0 × × × × ×

0 0 × × × ×

0 0 0 × × ×

0 0 0 0 × ×





 +u

0 0 q q q q .

Therefore, the first k columns of the decomposition can be written in the form

AU_k^(m−k)= U_k^(m−k)H_kk^(m−k)+ h_k+1,ku^(m−k)_k+1 e^T_k + β_kq_mkum+1e^T_k, where U_k^(m−k)consists of the first k columns of Um^(m−k), H_kk^(m−k) is the leading principal submatrix of order k of H_m^(m−k), and q_km is from the matrix Q = Q₁· · · Q_m−k.

(49)

Hence if we set

U˜_k = U_k^(m−k), H˜_k = H_kk^(m−k),

β˜_k = kh_k+1,ku^(m−k)_k+1 + β_kq_mku_m+1k₂,

˜

u_k+1 = β˜_k⁻¹(h_k+1,ku^(m−k)_k+1 + β_kq_mku_m+1), then

A ˜Uk= ˜UkH˜k+ ˜βku˜k+1e^T_k

is an Arnoldi decomposition whose starting vector is proportional to (A − κ₁I) · · · (A − κ_m−kI)u₁.

Avoid any matrix-vector multiplications in forming the new starting vector.

Get its Arnoldi decomposition of order k for free.

For large n the major cost will be in computing U Q.

(50)

師大 Krylov-Schur restarting

If a Krylov decomposition can be partitioned in the form A

U₁ U₂ = U₁ U₂

B₁₁ B₁₂ 0 B₂₂

+ u

b^H₁ b^H₂ , then

AU₁ = U₁B₁₁+ ub^H₁ is also a Krylov decomposition.

The process of Krylov-Schur restarting:

Compute the Schur decomposition of the Rayleigh quotient Move the desired eigenvalues to the beginning

Throw away the rest of the decomposition

(51)

Exchanging eigenvalues and eigenblocks

• Move an eigenvalue from one place to another.

Let a triangular matrix be partitioned in the form

R ≡





A B C

0 S D

0 0 E



,

where

S =

s₁₁ s₁₂ 0 s₂₂

. Suppose that Q is a unitary matrix such that

Q^HSQ =

s₂₂ sˆ₁₂ 0 s₁₁

,

(52)

then the eigenvalues s₁₁and s₂₂in the matrix

diag I Q^H I R diag I Q I =





A BQ C

0 Q^HSQ Q^HD

0 0 E



 will have traded places.

• How to find such unitary matrix Q?

Let

S =

S₁₁ S₁₂ 0 S22

,

where Siiis of order ni(i = 1, 2). Therefore are four cases to consider.

1 n1 = 1, n2= 1.

2 n₁ = 2, n₂= 1.

3 n1 = 1, n2= 2.

4 n₁ = 2, n₂= 2.

(53)

For the first two cases (n₁ = 1, n₂ = 1or n₁= 2, n₂= 1):

Let

S =

S₁₁ s₁₂ 0 s₂₂

,

where S₁₁is of order one or two. Let x be a normalized eigenvector corresponding to s₂₂and let Q = [x Y ] be orthogonal. Then

Q^TSQ =

x^T Y^T

S

x Y =

x^TSx x^TSY Y^TSx Y^TSY

=

s22 sˆ^T₁₂ 0 Sˆ11

.

Note that ˆS₁₁and S₁₁have the same eigenvalues.

(54)

For the third case (n₁ = 1, n2= 2):

Let

S =

s₁₁ s^T₁₂ 0 S22

,

where S₂₂is of order two. Let y be a normalized left eigenvector corresponding to s11and let Q = [X y] be orthogonal. Then Q^TSQ =

X^T y^T

S

X y =

X^TSX X^TSy y^TSX y^TSy

=

Sˆ22 ˆs12

0 s11

.

(55)

For the last case (n₁= 2, n₂ = 2):

Let

S =

S11 S12

0 S₂₂

. Let (S₂₂, X)be an orthonormal eigenpair, i.e.,

SX = X U S₂₂U⁻¹

for some nonsingular U , and let Q = [X Y ] be orthogonal. Then Q^TSQ =

X^TSX X^TSY Y^TSX Y^TSY

=

X^TXU S₂₂U⁻¹ X^TSY Y^TXU S22U⁻¹ Y^TSY

=

U S₂₂U⁻¹ Sˆ₁₂ 0 Sˆ11

.

(56)

Question

How to compute the orthonormal eigenbasis X?

Let the eigenbasis be

P I

, where P is to be determined.

Then

S₁₁ S₁₂ 0 S₂₂

P I

=

P I

S₂₂. Hence P can be solved from the Sylvester equation

S11P − P S22= −S12.

The orthonormal eigenbasis X can be computed by the QR factorization

P I

=

X Y

R 0

.

(57)

The Krylov-Schur cycle Assume A ∈ C^n×n.

1 Write the corresponding Krylov decomposition in the form AUm= UmTm+ βmum+1e^T_m.

2 Compute the Schur decomposition of T_m, Sm= Q^HTmQ

where S_m is upper triangular.

3 Transform the decomposition to the form A ˆU_m = ˆU_mS_m+ u_m+1b^H_m+1.

4 Select m − k Ritz values and move them to the end of S_m, accumulating the transformations in Q₁.

5 Truncate the decomposition, i.e.,

S_k := S_m[1 : k, 1 : k], b^H_k := b^H_m+1Q₁[:, 1 : k], U_k:= ˆU_mQ₁[:, 1 : k].

(58)

Deflation

We say a Krylov decomposition has been deflated if it can be partitioned in the form

A

U₁ U₂ = U₁ U₂

B₁₁ B₁₂ 0 B₂₂

+ u

0 b^H₂ . It implies that

AU₁ = U₁B₁₁, so that U11spans an eigenspace of A.

(59)

Criterion of Deflation:

Theorem

Let

AU = U B + ub^H

be an orthonormal Krylov decomposition, and let

[M, ˜U ] = [M, U W ]be an orthonormal pair. Let [W, W⊥]be unitary, and set

B =˜

W^H W_⊥^H

B

W W⊥ ≡

B˜₁₁ B˜₁₂ B˜₂₁ B˜₂₂

and

˜b^H = b^H

W W⊥ = ˜b^H₁ ˜b^H₂ .

(60)

Then

kA ˜U − ˜U M k²_F = k ˜B₂₁k²_F + k˜b₁k²_F + k ˜B₁₁− M k²_F. Proof: Let

U˜ U˜⊥ = U W W⊥ . Then

A ˜U − ˜U M = U BW + ub^HW − U W M

= U

W W⊥

W^H W_⊥^H

B

W W⊥

I 0

−

I 0

M

+ub^H

W W⊥

I 0

= U˜ U˜⊥

B˜11− M B˜21

+ u˜b^H₁ = U˜ U˜⊥ u





B˜11− M B˜₂₁

˜b^H₁



.

(61)

Since u^HU = 0, we have

u^H U˜ U˜⊥ = u^HU

W W⊥ = 0.

It implies that [ ˜U , ˜U⊥, u]is an orthonormal matrix. Therefore, kA ˜U − ˜U M k²_F = k





B˜₁₁− M B˜21

˜b^H₁



k²_F

= k ˜B₂₁k²_F + k˜b₁k²_F + k ˜B₁₁− M k²_F. Suppose that A ˜U − ˜U M is small. Transform the Krylov decomposition to the form

A U˜ U˜⊥ = ˜U U˜⊥

B˜₁₁ B˜₁₂ B˜₂₁ B˜₂₂

+ u ˜b^H₁ ˜b^H₂

= U˜ U˜⊥

B˜11 B˜12

0 B˜22

+ u

0 ˜b^H₂ + ˜U⊥B˜21+ u˜b^H₁

.

(62)

From Theorem 16, we have k

B˜21

˜b^H₁

k_F ≤ kA ˜U − ˜U M kF,

with equality if and only if M = W^HBW. Therefore, if the residual norm kA ˜U − ˜U M k_F is sufficiently small, we may set B˜21and ˜b1to zero to get the approximate decomposition

A U˜ U˜⊥ ≈ ˜U U˜⊥

B˜₁₁ B˜₁₂ 0 B˜₂₂

+ u

0 ˜b^H₂ . Rational Krylov transformations

Shift-and-invert transformations in Arnoldi’s method is to focus the algorithm on the eigenvalues near the shift κ.

How to do when it needs to use more than one shift?

Restart with a new shift and a new vector

Change a Krylov decomposition from one in (A − κ1I)⁻¹ to one in (A − κ2I)⁻¹.

(63)

Suppose we have a Krylov sequence

u, (A − κ₁I)⁻¹u, (A − κ₁I)⁻²u, · · · , (A − κ₁I)^1−ku.

Set v = (A − κ₁I)^1−ku, then the sequence with its terms in reverse order is

v, (A − κ₁I)v, · · · , (A − κ₁I)^k−1v, so that

K_k[(A − κ₁I)⁻¹, u] = K_k[A − κ₁I, v].

By the shift invariance of a Krylov sequence K_k[A − κ1I, v] = Kk[A − κ2I, v].

Set

w = (A − κ₂I)^k−1v, we have

K_k[A − κ₂I, v] = K_k[(A − κ₂I)⁻¹, w].