Download Gram-Schmidt Algorithm and QR Factorization: Different Views - Prof. P. Blomgren and more Study notes Mathematics in PDF only on Docsity!
RecapGram-Schmidt
Gram-Schmidt and Householder: Different Views of QR^ Numerical Matrix Analysis^ Lecture Notes #07 — The QR-Factorization and Least SquaresProblems: Gram-Schmidt and Householder
Peter Blomgren, 〈[email protected]
Department of Mathematics and Statistics^ Dynamical Systems GroupComputational Sciences Research Center^ San Diego State UniversitySan Diego, CA 92182-7720^ http://terminus.sdsu.edu/^ Spring 2010 Peter Blomgren, 〈[email protected]〉^ QR & LSQ: Gram-Schmidt and Householder
— (1/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR Outline^1 Recap^ Projectors: Orthogonal and Non-Orthogonal^ Classical Gram-Schmidt^2 Gram-Schmidt^ Bad News for the Classical Version^ Improving Gram-Schmidt^ Feel the Need for Speed?^3 Gram-Schmidt and Householder: Different Views of QR^ Gram-Schmidt — Triangular Orthogonalization^ Householder — Orthogonal Triangularization^ Householder vs. Gram-Schmidt^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (2/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Projectors: Orthogonal and Non-Orthogonal Classical Gram-Schmidt Last Time^ Orthogonal and non-orthogonal projectors
[ 2 P = P,^ P^ =^ P
]∗^.
Projection with an orthonormal, and arbitrary, basis^ P^ =
∗∗̂̂ Q Q,^ P^ =^ A(AA
−^1 ∗)A.
Rank-one projections, rank-(
m^ −^ 1) complementary projections∗∗ P = ˜q˜q, P=^ I^ −^ ˜q˜q.⊥^ QR-Factorization, using classical Gram-Schmidt orthogonalization.^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (3/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Projectors: Orthogonal and Non-Orthogonal Classical Gram-Schmidt Algorithm:^ Classical Gram-Schmidt^ Algorithm (Classical Gram-Schmidt)^ for j = 1:n^ ˜v=^ ˜aj^ j^ for i=1:(j-1)∗^ r=^ ˜q˜aij^ j^ i^ ˜v=^ ˜v−^ r˜qj^ j^ ij^ i^ endfor^ r=^ ‖˜v‖jj^ j^2 ˜q=^ ˜v/rj^ j^ jj^ endfor^ Mathematically, we are done. Numerically, however, we can runinto trouble due to roundoff errors.^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (4/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Projectors: Orthogonal and Non-Orthogonal Classical Gram-Schmidt Classical Gram-Schmidt: Revisitedm×n^ Let^ A^ ∈^ C,^ m^
≥^ n, be a full-rank matrix with columns
˜a.j^ With orthogonal projectors
Pwe can express the Gram-Schmidtj^ orthogonalization using the formulas^ ˜qj^
P˜aj^ j= ,^ j^ = 1,... , ‖P˜a‖j^ j^2 n The projector^ Pmust be anj^
m^ ×^ m-matrix of rank
m^ −^ (j^ −^ 1) which projects the space
m^ Corthogonally onto the space orthogonal to^ 〈˜q,... ,^1
˜q〉. (P=^ I^ ).j−^11 Note that^ ˜q∈ 〈˜a,... ,j^1
˜a〉^ and^ ˜q⊥ 〈˜q,... ,j^ j^1
˜q〉; therefore thisj−^1 description is equivalent to the algorithm on slide 4.We can represent the projector
∗̂̂ P= I − QQwherej j− 1 j−^1 ̂ Qisj−^1 the^ m^ ×^ (j^ −^ 1)-matrix [
˜q˜q...^ ˜q].^1 2 j−^1 Peter Blomgren, 〈[email protected]〉^ QR & LSQ: Gram-Schmidt and Householder
— (5/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical Version Improving Gram-SchmidtFeel the Need for Speed? Classical Gram-Schmidt: The Bad News
1 of 2
Unfortunately, classical Gram-Schmidt is not numerically stable —in finite precision, the vectors
˜qmay lose orthogonality...j^ 10 20 30 40 50 60 70 8010 20 30 40 50 60 70 80
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1^10 20 30 40
1 0.9 10 0.8 20 0.7 30 0.6 40 0.5 0.4 50 0.3 60 0.2 70 0.1 80 060 70 80 ∗ Figure: Comparing QQ^ (which should be the identity matrix) for classical (left) andmodified (right) Gram-Schmidt on a particularly hard problem where
−^1 λ= 2and 1 −^80 λ= 2.^ We see that C-GS completely loses orthogonality after 20-some steps; (^80) whereas GSdoes not suffer this catastrophic breakdown.mod^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (6/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical Version Improving Gram-SchmidtFeel the Need for Speed? Classical Gram-Schmidt: The Bad News
2 of 2 0 10 20 30 40
50 60 70 80 010 −5 10 −10 10 −15 10 −20 10 Classical Gram−SchmidtModified Gram−SchmidtCorrect Eigenvalues−25 (^10) Figure: The performance of classical (blue ’o’s
) and modified (red ’x’s) Gram-Schmidt on a particularly hard problem where
(^11) λ= and^ λ=^. C-GS identifies the first 1 80 80 2 2
∼
26 eigenvalues (down to the size
√ ∼ ǫ), whereas M-CG identifiesmach
∼^ 54 eigenvalues (down to the size^ ∼^ ǫ).mach^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (7/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical Version Improving Gram-SchmidtFeel the Need for Speed? Description: The Test Problem
Matlab-centric First we let^ U^ and^ V^ be two randomly selected 80
×^ 80 unitary matrices [U,X] = qr(randn(80));
[V,X] = qr(randn(80)); Then we build a matrix
−^1 A with eigenvalues 2 −^2 −^80 , 2 ,... ,^2 : S = diag(2.̂^ (-1:-1:-N));
A = U * S * V’; Finally we compute the QR-factorization using both classical andmodified Gram-Schmidt^ [QC,RC] = qr
cgs(A);^ [QM,RM] = qr
mgs(A); Now, the diagonals of^ RM
and^ RC^ contain the recovered eigenvalues. — This is what was plotted on the previous slide. Burning Question:^ What is the modified Gram-Schmidt method?!?^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (8/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical VersionImproving Gram-Schmidt Feel the Need for Speed? Counting Work: The Memory Access Latency Model^ If we have three cache-levels (L1, L2, and L3), some averagehit-rate (and hence miss-rate) for each level and the time it takesto access that cache-level (the hit-cycle-time), then we end upwith a measure for the average memory access latency per memoryaccess^ T^ ∼^ +^ (L
hit^ rate^ ∗^ L1^ hit^ cycle
time)
∗^ L2^ hit^ cycle^ time
∗^ L3^ hit^ cycle^ time
- (L3^ miss^ rate^ ∗^ [S]DRAM
latency) If this does not scare you, please get a Ph.D. in algorithm designon the compiler / silicon level !!!Meanwhile, the rest of us will count
“flops”,^ i.e.^ floating-point operations (multiplications and additions)!^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (13/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical VersionImproving Gram-Schmidt Feel the Need for Speed? Counting Work: Gram-Schmidt Orthogonalization^ Theorem (Computational Complexity of Modified Gram-Schmidt)^ The modified Gram-Schmidt orthogonalization algorithm requires
(^2) ∼ 2 mnflops to compute the QR-factorization of an m
×^ n matrix. Here we have assumed that complex arithmetic is just as fast asreal arithmetic. This is not true in general.^ c·^ c=^1
[r·^ r−^ i·^ i] +^ i^ [^1 2 1
r·^ i+^ r·^ i]^1 2 2 c+^ c=^ [r+^ r] +^1 2 1
i^ [i+^ i]^1 Hence, the complex multiplication consists of 4 real multiplicationsand 2 real additions; and the complex addition consists of 2 realadditions. Also, we need
at least^ double the amount of memory accesses.^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (14/31)
Recap Gram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Bad News for the Classical VersionImproving Gram-Schmidt Feel the Need for Speed? Counting Flops^ The Outer Loop:
for i=1:n MMThe Inner Loop:^ for j=(i+1):n MMMMris formed by an^ m-inner product – requiringij
m^ multipli- cations and (m^ −^ 1) additions vrequires^ m^ multiplications andj^
m^ subtractions MMEnd Inner LoopEnd Outer Loop^ Work
n^ n∑∑ ∼^4 m i=1^ j=i+1 n∑ ∼^4 m(n^ −^ i) i=1^2 2 ∼ 4 mn−^4 mn/^22 ∼ 2 mn Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (15/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal TriangularizationHouseholder vs. Gram-Schmidt Gram-Schmidt as Triangular Orthogonalization
1 of 3
Each outer loop in the modified Gram-Schmidt algorithm can beseen as a right-multiplication by a square upper triangular matrix. E.g. Iteration#1^223 |^ |^ ||^ |^ ||^ |^ |^66 |^ |^ |^76676 ˜v˜v...^ ˜v^61 2 n^76454 |^ |^ ||^ |^ ||^ |^ ||^ |^ |
(^31) rr (^1213) − −... r r r 11 11 11 7177 = 175.. .| {z } R 1 23 |^ |^ ||^ |^ ||^ |^ ||^ |^ | 6767 (2)(2) 6 ˜q˜v...^ ˜v^71 n^2 45 |^ |^ ||^ |^ ||^ |^ ||^ |^ | Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (16/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal TriangularizationHouseholder vs. Gram-Schmidt Gram-Schmidt as Triangular Orthogonalization
2 of 3
E.g. Iteration#2^2 |^ |^ ||^ |^ ||^ |^ ||^ |^ |^66 (2)(2)^6 ˜q˜v...^ ˜v^1 n^2 4 |^ |^ ||^ |^ ||^ |^ ||^ |^ |
2331 r^123677 −^... 6 r^ r^7722 22 677167545.. .|^ {z^ }^ R^2
23 |^ |^ ||^ |^ ||^ |^ ||^ |^ | 6767 (3) = 6 ˜q˜q...^ ˜v^71 2 n 45 |^ |^ ||^ |^ ||^ |^ ||^ |^ | When we are done we have^ A R^1
̂R... R = Q^ ⇔^ A 2 n ︸ ︷︷ ︸ b− (^1) R ̂̂ = Q^ R Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (17/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal TriangularizationHouseholder vs. Gram-Schmidt Gram-Schmidt as Triangular Orthogonalization
3 of 3
This formulation of the QR-factorization shows that we can
think of the modified Gram-Schmidt algorithm as a method of triangular orthogonalization
We apply a sequence of triangular operations from the right of thematrix^ A^ in order to reduce it to a matrix
̂ Q^ with orthonormal columns.In practice we^ do not
explicitly form the matrices
Rand multiplyi^ them together.However, this view tells us something about the structure ofmodified Gram-Schmidt. Note:^ From now on when we say “Gram-Schmidt” we mean themodified version.^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (18/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal Triangularization Householder vs. Gram-Schmidt Householder Triangularization^ Householder triangularization is another way of computing theQR-factorization:^ Gram-Schmidt^
Householder Numerically stable^
Even better stability Useful for iterative methods
Not as useful for iterative methods “Triangular Orthogonalization”
“Orthogonal Triangularization” ̂ ARR... R= Q^12 n Q...^ QQA^ =^ Rn^21 Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (19/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal Triangularization Householder vs. Gram-Schmidt Householder Triangularization
By Picture 232 ×^ ×^ ×^ ∗^ ×^ ×^ × 676 Q 67616 ×^ ×^ ×^7 →^6454 ×^ ×^ ××^ ×^ ×|^ {z^ }^ A
(^32) ∗ ∗^ ×^ ×^ × 0 ∗ ∗^ ∗^ ∗ 76 Q (^7620) ∗ ∗ 7 →^60 ∗ (^540) ∗ ∗^0 ∗ 0 ∗ ∗^0 ∗| {z } QA 1 323 ×^ ×^ ××^ × 767 Q 76737 →^6 ∗^754500 | {z }|^ {z^ } QQA^ QQQA 21321 represents a new zero. 0 ∗ represents a modified entry. × represents an unchanged entry. The Big Question:^ How do we find the unitary matrices
Q?!?i^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (20/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal Triangularization Householder vs. Gram-Schmidt Algorithm:^ Householder QR-Factorization^ Algorithm (Householder QR-Factorization)^ for k = 1:n^ ˜x^ =^ A(k:m,k)^ ˜v=^ sign(x)‖˜x‖k^1
˜e+^ ˜x (^21) ˜v= ˜v/‖˜v‖k k k 2 ∗ A(k:m,k:n) = A(k:m,k:n)^ −^2 ˜v(˜vA(k:m,k:nk^ k^
)) endfor A(k:m,k)^ Denotes the
kth thru^ mth rows, in the
kth column of^ A^ — a vector quantity. A(k:m,k:n) Denotes the^ kth thru
mth rows, in the^ kth thru
nth columns of^ A^ — a matrix quantity. Peter Blomgren, 〈[email protected]〉^ QR & LSQ: Gram-Schmidt and Householder
— (25/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal Triangularization Householder vs. Gram-Schmidt Householder-QR: Where is the Q ?!?
1 of 2
At the completion of the Householder QR-factorization, themodified matrix^ A^ contains
R^ (of the full QR-factorization), but
Q
is nowhere to be found.Often, we only need^ Q
implicitly, as in the^ action
of^ Q^ on something.^ I.e.^ if we need
∗˜ Qb, we can add the line ˜˜b(k:m) =^ b(k:m)^ −^2 ˜
∗˜v(˜vb(k:m))k k^ to the loop, or store the generated vectors
˜v, and^ a posteriorik^ compute for k = 1:n^ ˜˜b(k:m) =^ b(k:m)^ −
∗˜ 2 ˜v(˜vb(k:m))k k^ endfor^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (26/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular Orthogonalization Householder — Orthogonal Triangularization Householder vs. Gram-Schmidt Householder-QR: Where is the Q ?!?
2 of 2
If we need^ Q˜x, then we must store the generated vectors
˜v, andk^ compute for k = n:-1:1^ ˜x(k:m) =^ ˜x(k:m)^ −
∗ 2 ˜v(˜v˜x(k:m))k k^ endfor We can also use this algorithm to explicitly generate
Q
Q^ =^ In×n for k = n:-1:1^ Q(k:m,k:n) =^ Q(k:m,k:n
∗) − 2 ˜v(˜vQ(k:m,k:nk k^
endfor^ Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (27/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular OrthogonalizationHouseholder — Orthogonal Triangularization Householder vs. Gram-Schmidt Comparison: Householder vs. Gram-Schmidt (modified)^0
20 30 40 50
60 70 80 010 −5 10 −10 10 −15 10 −20 10 HouseholderModified Gram−SchmidtCorrect Eigenvalues−25 (^10) Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (28/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular OrthogonalizationHouseholder — Orthogonal Triangularization Householder vs. Gram-Schmidt Q-Orthogonality: Householder, Modified-GS, and Classical-GS^10 20
1 0.9 10 0.8 20 0.7 30 0.6 40 0.5 0.4 50 0.3 60 0.2 70 0.1 80 040 50 60 70 80 10 20 30 40 50 60 70 80 10 20 30 40 50 60 70 80
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 10 20 30 40 50 60 70 8010 20 30 40 50 60 70 80
1 0.9 0.8 0.7^ Figure:^ Entries of^ Q 0.6 0.5 0.4 0.3 0.2 0.
∗Q^ for House-holder (top-left), GS(top-right)mod^ and classical (left) Gram-Schmidt. Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (29/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular OrthogonalizationHouseholder — Orthogonal Triangularization Householder vs. Gram-Schmidt Householder-QR: Work^ The dominating work is done in the operation^ A(k:m,k:n) =^
∗A(k:m,k:n) − 2 ˜v(˜vk k A(k:m,k:n)) Each entry in^ A(k:m,k:n
) is “touched” by 4 flops per iteration ( from the inner product, 1 scalar multiplication, and 1 subtraction).The size of the sub-matrix
A(k:m,k:n) is (m^ −^ k^ + 1)^ ×^ (n^ −^ k
- 1), so we get n ∑^ (m−k+1)(n−k+1) k=
n∑ ∼^ (m−k)(n−k)^ ∼ k= n∑^ (^2 mn^ +^ k−^ k(m^ + k=
) n) (^32) nn (^2) ∼ mn+ −^ (m^3
(^23) mnn+ n) ∼ −^2 Hence, the work of Householder-QR is
(^3) 2n (^2) ∼ 2mn− flops. 3 Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (30/31)
RecapGram-Schmidt Gram-Schmidt and Householder: Different Views of QR
Gram-Schmidt — Triangular OrthogonalizationHouseholder — Orthogonal Triangularization Householder vs. Gram-Schmidt The Gram-Schmidt Song
c©2000 Rebecca Hartman-Baker (to the tune of ”When Johnny Comes Marching Home”)The wonderful Gram-Schmidt algorithm, hurrah, hurrah!Is something you compute QR with, hurrah, hurrah!Triangular orthogonalization,This is the way Gram-Schmidt is doneand both Q and R are explicitly formed and stored.The classical Gram-Schmidt incarnate, hurrah, hurrah!It really isn’t all that great, hurrah, hurrah!For its leftward look and dependence on Alose orthogonality in computed q
k And you need another matrix-space to store Q.The modified Gram-Schmidt incarnate, hurrah, hurrah!Is better than its Classical mate, hurrah, hurrah!For advantages of looking rightare numerical accuracy and an overwrite-ing of A with Q, so bully for modified!Now we know Gram-Schmidt algorithms, hurrah, hurrah!We know the good and bad of them, hurrah, hurrah!But unless I need an explicit Q,Gram-Schmidt is not the best to dofor Householder is much superior.From^ http://www.cse.uiuc.edu/
∼rjhartma/gssong.html^ without permission. Peter Blomgren,^ 〈[email protected]
〉^ QR & LSQ: Gram-Schmidt and Householder
— (31/31)
