ReflRotn March 13, 2006 7:43 amProf. W. Kahan Math. 128B Page 1/6 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 M ATLAB 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) .ReflRotn March 13, 2006 7:43 amProf. W. Kahan Math. 128B Page 2/6 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·…·W2·W1·F = in which each reflection Wj = Wj' = Wj–1 is constructed to annihilate all subdiagonal elements in column j of Fj–1 := Wj-1·…·W2·W1·F . Then Q := W1·W2·…Wn· . Each Wj = I – wj·wj' has wj'·wj = 2 (or 0 ) and no nonzero element in wj above row j . Each Fj = Wj·Fj–1 has the same first j–1 rows as Fj–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 F0 := 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 := G1 .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 RotationsA Givens Rotation is Q := so chosen that a 2-vector v = is rotated to Q·v = wherein |r|2 = v'·v , so c2 + 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) .ROIOIOcss'–cxyr0ReflRotn March 13, 2006 7:43 amProf. W. Kahan Math. 128B Page 3/6Bottom-Up 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 Givens rotations thus:For 1 ≤ i ≤ m–1 and 1 ≤ j ≤ n let Qij be the Givens rotation that acts upon an m-by-n matrix Z to overwrite Qij· = onto . We shall premultiply F by a sequence of rotations Qij in this order (from right to left):for j = 1 up to n in turn