SLIDE 1
Finding Eigenvalues: Arnoldi Iteration and the QR Algorithm Will - - PowerPoint PPT Presentation
Finding Eigenvalues: Arnoldi Iteration and the QR Algorithm Will - - PowerPoint PPT Presentation
Finding Eigenvalues: Arnoldi Iteration and the QR Algorithm Will Wheeler Algorithms Interest Group Jul 13, 2016 Outline Finding the largest eigenvalue Largest eigenvalue power method So much work for one eigenvector. What about
SLIDE 2
SLIDE 3
Power Method
How do we find the largest eigenvalue of an m × m matrix A?
◮ Start with a vector b and make a power sequence:
b, Ab, A2b, . . .
◮ Higher eigenvalues dominate:
An(v1 + v2) = λn
1v1 + λn 2v2 ≈ λn 1v1
Vector sequence (normalized) converged?
◮ Yes: Eigenvector ◮ No: Iterate some more
SLIDE 4
Krylov Matrix
Power method throws away information along the way.
◮ Put the sequence into a matrix
K = [b, Ab, A2b, . . . , An−1b] (1) = [xn, xn−1, xn−2, . . . , x1] (2)
◮ Gram-Schmidt approximates first n eigenvectors
v1 = x1 (3) v2 = x2 − x2 · x1 x1 · x1 x1 (4) v3 = x3 − x3 · x1 x1 · x1 x1 − x3 · x2 x2 · x2 x2 (5) . . . (6)
◮ Why not orthogonalize as we build the list?
SLIDE 5
Arnoldi Iteration
◮ Start with a normalized vector q1
xk = Aqk−1 (next power) (7) yk = xk −
k−1
- j=1
(qj · xk) qj (Gram-Schmidt) (8) qk = yk/ |yk| (normalize) (9)
◮ Orthonormal vectors {q1, . . . , qn} span Krylov subspace
(Kn = span{q1, Aq1, . . . , An−1q1})
◮ These are not the eigenvectors of A, but make a similarity
(unitary) transformation Hn = Q†
nAQn
SLIDE 6
Arnoldi Iteration
◮ Start with a normalized vector q1
xk = Aqk−1 (next power) (10) yk = xk −
k−1
- j=1
(qj · xk) qj (Gram-Schmidt) (11) qk = yk/|yk| (normalize) (12)
◮ Orthonormal vectors {q1, . . . , qn} span Krylov subspace
(Kn = span{q1, Aq1, . . . , An−1q1})
◮ These are not the eigenvectors of A, but make a similarity
(unitary) transformation Hn = Q†
nAQn
SLIDE 7
Arnoldi Iteration - Construct Hn
We can construct the matrix Hn along the way.
◮ For k ≥ j: hj,k−1 = qj · xk = q† j Aqk−1 ◮ For k = j: hj,k−1 = |yk|
(equivalent)
◮ For k < j: hj,k−1 = 0 ◮ ⇒ Hn is upper Hessenberg
Hn = h1,1 h1,2 h1,3 h1,4 h1,5 h1,6 h2,1 h2,2 h2,3 h2,4 h2,5 h2,6 h3,2 h3,3 h3,4 h3,5 h3,6 h4,3 h4,4 h4,5 h4,6 h5,4 h5,5 h5,6 h6,5 h6,6 (13)
◮ Why is this useful? Hn = Q† nAQn → same eigenvalues as A
SLIDE 8
Finding Eigenvalues
How do we find eigenvalues of a matrix?
◮ Start with a triangular matrix A ◮ The diagonal elements are the eigenvalues
A = a11 a12 a13 a14 a15 a16 a22 a23 a24 a25 a26 a33 a34 a35 a36 a44 a45 a46 a55 a56 a66 (14) What if A isn’t triangular?
SLIDE 9
QR algorithm
Factorize A = QR
◮ Q is orthonormal (sorry this is a different Q) ◮ R is upper triangular (aka right triangular)
Algorithm:
- 1. Start A1 = A
- 2. Factorize Ak = QkRk
- 3. Construct Ak+1 = RkQk
= Q−1
k QkRkQk = Q−1 k AkQk
← similarity transform
- 4. Ak converges to upper triangular
SLIDE 10
QR Algorithm Example (numpy.linalg.qr)
A1 = 0.313 0.106 0.899 0.381 0.979 0.375 0.399 0.488 0.876 A2 = 1.649 0.05 −0.256 0.035 0.502 0.534 −0.014 0.004 0.017 A3 = 1.653e + 00 2.290e − 02 2.305e − 01 5.966e − 03 5.063e − 01 −5.342e − 01 8.325e − 05 −9.429e − 05 9.666e − 03 A4 = 1.653e + 00 1.872e − 02 −2.285e − 01 1.800e − 03 5.063e − 01 5.349e − 01 −4.812e − 07 1.803e − 06 9.554e − 03
SLIDE 11
QR Decomposition - Gram-Schmidt
How do we get the matrices Q and R? http://www.seas.ucla.edu/ vandenbe/103/lectures/qr.pdf Recursively solve for the first column of Q and the first row of R:
◮ A = QR ◮
a1 A2
- =
- q1
Q2 r11 R12 R22
- ◮
a1 A2
- =
- r11q1
q1R12 + Q2R22
- ◮ r11 = ||a1||, q1 = a1/r11
◮ Note q⊺ 1q1 = 1 and q⊺ 1Q2 = 0 ◮ R12 = q⊺ 1A2 ◮ Solve A2 − q1R12 = Q2R22
SLIDE 12
QR Decomposition of Hessenberg Matrix
What if we have a matrix in upper Hessenberg form? Hn = h1,1 h1,2 h1,3 h1,4 h1,5 h1,6 h2,1 h2,2 h2,3 h2,4 h2,5 h2,6 h3,2 h3,3 h3,4 h3,5 h3,6 h4,3 h4,4 h4,5 h4,6 h5,4 h5,5 h5,6 h6,5 h6,6 (15) We just need to remove the numbers below the diagonal by combining rows → Givens rotation.
SLIDE 13
Givens Rotation
A Givens rotation acts only on two rows and leaves the others unchanged. G1 = γ −σ σ γ 1 1 1 (16) Want G a b
- =
r
- ,
r = √ a2 + b2
◮ γ = a/
√ a2 + b2
◮ σ = −b/
√ a2 + b2
SLIDE 14
QR Decomposition of Hessenberg Matrix
G2G1Hn =
- h11
- h12
- h13
- h14
- h15
- h22
- h23
- h24
- h25
- h33
- h34
- h35
h43 h44 h45 h54 h55 (17) So Hn = G⊺
1G⊺ 2 . . . G⊺ n−1R = QR
SLIDE 15
QR Decomposition of Hessenberg Matrix
QR algorithm is iterative; is RQ also upper Hessenberg? Yes: Acting on the right, the Givens rotations mix two columns instead of rows, but change the same zeros. G2G1Hn =
- r11
- r12
- r13
r14 r15
- r21
- r22
- r23
r24 r25
- r32
- r33
r34 r35 r44 r45 r55 (18)
SLIDE 16
Householder Reflections
◮ Arnoldi finds a few eigenvalues of a large matrix ◮ In practice, QR uses Householder reflections ◮ www.math.usm.edu/lambers/mat610/sum10/lecture9.pdf
Reflection across the plane perpendicular to a unit vector v: P = I − 2vv⊺ (19) Px = x − 2vv⊺x (20) Px − x = −v(2v⊺x) (21) Want to arrange first column into zeros: find v so that Px = αe1.
SLIDE 17
Householder Reflections
Want to arrange first column into zeros: find v so that Px = αe1. ||x|| = α (22) x − 2vv⊺x = αe1 (23) 1 2(x − αe1) = v(v⊺x) (24) So we know
◮ v ∝ x − αe1 ◮ v is a unit vector
v = x − αe1 ||x − αe1|| (25)
SLIDE 18
Lanczos
What if we apply the Arnoldi algorithm to a Hermitian matrix?
◮ Start with a normalized vector q1
xk = Aqk−1 (next power) (26) yk = xk −
k−1
- j=k−2
(qj · xk) qj (Gram-Schmidt) (27) qk = yk/ |yk| (normalize) (28) Since Hn = Q†
nAQn, ◮ Hn is also Hermitian ◮ Entries of Hn above the first diagonal are 0 (tridiagonal) ◮ We don’t need to calculate them!
SLIDE 19