Download Householder & Givens Rotations: Reflections & Orthogonal Transformations in Linear Algebra and more Study notes Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
Reflections, Rotations and QR Factorization
QR Factorization figures in Least-Squares problems and Singular-Value Decompositions
among other things numerical. These notes explain some reflections and rotations that do it, and
offer MATLAB implementations; in its notation, x ' := (complex conjugate transpose of x).
Householder Reflections
A Householder Reflection is W = I – w·w ' = W ' = W–1^ for any column w satisfying w ' ·w = 2.
If y = W·x then y ' ·y = x ' ·x , and y ' ·x = x ' ·y is real even if x, y and w are complex. W is a
reflection because W·w = –w and W·p = p whenever w ' ·p = 0. Numerical analysts name this
reflection after Alston S. Householder because he introduced it to them in the mid 1950s as part
of an improved way to solve Least-Squares problems.
Given columns x and e := [1, 0, 0, …, 0] ' so that e ' ·x = x 1 , we seek column w so that
w ' ·w = 2 and W := I – w·w ' reflects x to W·x = e·ß for some scalar ß. Its |ß| = ||x|| := √(x ' ·x)
and x ' ·W·x = x ' ·e·ß = x 1 ' ·ß must be real, so ß = ±||x||·x 1 /|x 1 | if x 1 ≠ 0.
Construction: Set Ñ := ||x|| , ß := ±Ñ·x 1 /|x 1 | , d := x – e·ß , and w := d/√(d ' ·d/2) unless
d = o , in which case set w := o. But all bets are off if UNDERFLOW degrades x ' ·x.
Proof: Let p := x + e·ß so that p ' ·d = 0 = p ' ·w. Then W·d = –d and W·p = p , whereupon
2W·x = W·(p+d) = p–d = 2e·ß as desired.
How is the sign ± in ß chosen? The simplest way maximizes Ω^2 := d ' ·d/2 = Ñ·(Ñ – (±)|x 1 |)
by setting ß := –Ñ·x 1 /|x 1 | , as we’ll see. Of course, any ± sign works when x 1 = 0.
Detailed Construction: Let v := x – e·x 1 , so that e ' ·v = v 1 = 0 , and let μ := v ' ·v > 0 , so that
Ñ := √(x ' ·x) = √( μ + |x 1 |^2 ). Next set ç := x 1 /|x 1 | = sign(x 1 ) except that we reset ç := 1 if
x 1 = 0. Next we choose ß := ±ç·Ñ. Numerical stability requires two cases to be distinguished:
If ß = –ç·Ñ set d := x – e·ß = x + e·ç·Ñ by copying x to d and then resetting
d 1 := x 1 – ß = ç·(|x 1 | + Ñ).
If ß = +ç·Ñ set d := x – e·ß = x – e·ç·Ñ by copying x to d and then resetting
d 1 := x 1 + ß = –ç·μ/(|x 1 | + Ñ).
Next Ω := √ ( (|d 1 |^2 + μ)/ 2 ) = √ ( |d 1 |·Ñ ) and w := d/Ω. Return [w, ß] = hshldrw(x).
QR Factorization:
Given an m-by-n matrix F with no fewer rows than columns (so m ≥ n ), we wish to factorize
F = Q·R , with Q ' ·Q = I and R upper-triangular, by using Householder reflections thus:
Wn ·…·W 2 ·W 1 ·F = in which each reflection Wj = Wj ' = Wj–1^ is constructed to annihilate all
subdiagonal elements in column j of Fj–1 := Wj-1 ·…·W 2 ·W 1 ·F. Then Q := W 1 ·W 2 ·…Wn ·.
Each W j = I – wj ·wj ' has w j ' ·wj = 2 (or 0 ) and no nonzero element in w j above row j. Each
Fj = Wj ·Fj–1 has the same first j–1 rows as F j–1 and no nonzero subdiagonal elements in its first
j columns. Each Qj := Wj ·Wj+1 ·…Wn · has ones on the diagonal and zeros elsewhere in its
first j–1 rows and columns, so Qj = Wj ·Qj+1 is obtained by altering only rows and columns of
Qj+1 with indices no less than j.
Detailed Construction in MATLAB:
Start with F 0 := F. For j = 1, 2, …, n in turn get [wj , ßj ] := hshldrw(Fj (j:m, j)) as above, and
store wj in place of Fj (j:m, j) ; then if j < n overwrite Fj (j:m, j+1:n) – wj ·(wj ' ·Fj (j:m, j+1:n))
onto Fj (j:m, j+1:n) to get Fj+1 (j:m, j+1:n).
Next, R := Diag([ß 1 , ß2 , …, ßn ]) + triu(Fn (1:n, 1:n), 1).
Finally, set Gn+1 := Fn and, for j = n, n–1, …, 1 in turn, extract wj from Gj+1 (j:m, j) ,
overwrite column [oj–1 ; 1; om–j ] onto Gj+1 (:, j) , and then onto Gj+1 (j:m, j:n) overwrite
Gj (j:m, j:n) := Gj+1 (j:m, j:n) – wj ·(wj ' ·Gj+1 (j:m, j:n)). Then Q := G 1.
Return [Q, R] = hshldrqr(F).
Numerical experiments indicate that MATLAB uses the same method to get [Q, R] = qr(F, 0).
QR Factorization by Givens Rotations
A Givens Rotation is Q := so chosen that a 2-vector v = is rotated to Q·v =
wherein |r|^2 = v ' ·v , so c^2 + s ' ·s = 1 when (by convention) we choose c ≥ 0. Here v ' is the
complex conjugate transpose of v , and s ' is the complex conjugate of s. The rotation is named
after Wallace Givens who introduced this rotation to numerical analysts in the 1950s while he
was working at Argonne National Labs near Chicago. The rotation is encoded in one complex
number t := (y/x) ' from which are derived c := 1/√(1 + t ' ·t) , s := c·t and r := x/c. In the
special case that t = ∞ (presumably because x = 0 ), we set c := 0 , s := 1 and r := y. In any
event, note that Q–1^ = Q '. Return [c, s, t, r] = givenst(x, y).
R O I O I O
c s
x y
r 0
premultiplications. After the last premultiplication we find R in the product’s first n rows and
columns after ignoring the subdiagonal elements that hold tij s. Then these are used to construct
Q as a product of inverse rotations Qij ' premultiplying in this reverse order (from right to
left):
for j = n down to 1 in turn { for i = m down to j+1 in turn { premultiply by Qij ' }}.
Each premultiplication by Qij ' affects only rows i and j of columns j to n of the product after
tij was extracted from location (i, j) and replaced by 0.
Return [Q, R] = gvnsdnqr(F).
MATLAB appears to use Householder reflections to get its [Q, R] = qr(F, 0). There is no reason
to expect any two of the three different [Q, R] factorizations to agree though they must be related
in the absence of roundoff: R 1 ·R2–1^ = Q 1 ' ·Q 2 must be a diagonal unitary matrix.
MATLAB Programs
function [F, R] = hshldrqr(F) % [Q, R] = hshldrqr(F) uses Householder Reflections to % factorize F = QR so that R is upper-triangular and % Q has orthonormal columns; Q'Q = I. This works only % if F has no more columns than rows, and if underflow % does not degrade F'F. Uses hshldrw.m. [m, n] = size(F) ; if (m < n), error(' F has more columns than rows in hshldrqr(F).'), end z = zeros(1, n) ; w = zeros(m, 1) ; for j = 1:n [w, z(j)] = hshldrw(F(j:m, j)) ; F(j:m, j) = w ; if (j < n), F(j:m, j+1:n) = F(j:m, j+1:n) - w(w'F(j:m, j+1:n)) ; end, end % ... j = 1:n R = diag(z) + triu(F(1:n, 1:n), 1) ; for j = n:-1: w = F(j:m, j) ; F(:, j) = zeros(m,1) ; F(j, j) = 1 ; F(j:m, j:n) = F(j:m, j:n) - w(w'*F(j:m, j:n)) ; end % ... j = n:-1:
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
I O
function [w, z] = hshldrw(x) % [w, z] = hshldrw(x) yields w with w'·w = 2 or 0 , % so W = I - ww' = W' = W^-1 reflects the given column % x to Wx = [z; 0; 0; ...; 0] with |z| = norm(x). % But all bets are off if UNDERFLOW degrades x'x. w = x(:) ; m = length(w) ; x1 = w(1) ; a1 = abs(x1) ; if (m < 2), w = 0 ; z = x1 ; return, end if (a1), s = x1/a1 ; else s = 1 ; end vv = w(2:m)'w(2:m) ; ax = sqrt(a1a1 + vv) ; z = -sax ; a1 = a1 + ax ; w(1) = sa1 ; dd2 = a1ax ; if (dd2), w = w/sqrt(dd2) ; end %... so w'*w = 2 unless w = o.
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
function [c, s, t, r] = givenst(x, y) % [c, s, t, r] = givenst(x, y) satisfies c >= 0 , % c^2 + |s|^2 = 1 , r = cx + sy , t' = x/y = s'/c. % So [c s]* [x] = [r] and c = 1/sqrt(1 + t't) % [-s' c] [y] [0] s = ct , r = x/c. if (x ~= 0) t = conj(y/x) ; u = sqrt(1 + t't) ; r = ux ; c = 1/u ; s = c*t ; else t = inf ; r = y ; c = 0 ; s = 1 ; end
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
function [F, R] = gvnsupqr(F) % [Q, R] = gvnsupqr(F) uses Givens Rotations to % factorize F = QR so that R is upper-triangular and % Q has orthonormal columns; Q'Q = I. This works only % if F has no more columns than rows, and if underflow % does not degrade F'F. Uses givenst.m bottom-up. [m, n] = size(F) ; if (m < n), error(' F has more columns than rows in gvnsupqr(F).'), end for j = 1:n , for i = m-1:-1:j [c, s, F(i+1, j), F(i, j)] = givenst(F(i, j), F(i+1, j)) ; if (j < n), F(i:i+1, j+1:n) = [c, s; -s', c]F(i:i+1, j+1:n) ; end end, end % ... i = m-1:-1:j , j = 1:n R = triu(F(1:n, 1:n)) ; for j = n:-1:1 , F(1:j, j) = zeros(j,1) ; F(j,j) = 1 ; for i = j:m- t = F(i+1, j) ; F(i+1, j) = 0 ; c = 1/sqrt(1 + t't) ; if (c~=0), s = ct ; else s = 1 ; end F(i:i+1, j:n) = [c, -s; s', c]*F(i:i+1, j:n) ; end, end % ... i = j:m-1, j = n:-1:
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
