3.3QR分解.97.则Q也是正交矩阵，且A=(Hn-1...H2Hi)-1R=QR.以上就是基于Householder变换的OR分解的具体实现过程.最后所得到的上三角矩阵R就存放在A的上三角部分如果不需要生成Q，则基于Householder变换的QR分解的运算量大约为2mn?-2n3/3.矩阵Q可通过下面的累积方法来计算：Q= Im,Q=QHk,k=1,2,..,n如果保留了每一步的Householder向量，则Q也可以通过下面的向后累积方法来计算Q= Im,Q=HQ,k=n,n-1,...,1.这样做的好处是一开始Q会非常稀疏（左上角为单位矩阵），随着迭代的进行，Q才会慢慢变满.而前面的计算方法，第一步就将Q变成了一个满矩阵采用向后累积方法计算Q的运算量大约为4mn一4mn2+43/3.如果只需要计算Q的前n列，则运算量大约为2mn2-2n3/3,此时QR分解的总运算量为4mn2-4n3/3.若m=n，则为8n3/3.算法3.5.基于Householder变换的OR分解% Given AE Rmxn, compute Q and R such that A= QR where Q E Rmxm and RE Rmxn% The upper triangular part of R is stored in the upper triangular part of A1: SetQ= Imxm2: for k =1 to n doT= A(k :m,k)3:4:[βk, vs] = House(r)5:A(k:m,k:n)=(Im-k+1-Biu)A(k:m,k:n)6= A(k:m,k :n)- βiuk(TA(k : m,k : n))7:Q(:,k :m) = Q(:,k :m)(Im-k+1 -βkvsuT)8:= Q(:,k: m) - βk(Q(s, k : m)uk)uT9: end for上面的算法只是关于HouseholderQR分解的一个简单描述，并没有考虑运算量问题.在实际计算时，我们通常会保留所有的Householder向量（用于向后累积法计算Q)由于第k步中H，所对应的Householder向量的长度为m－k+1,因此我们先把u单位化，使得uk的第一元素为1，这样就只要存储s(2：end)，共m一k个元素于是我们就可以把所有的
3.3 QR 分解 · 97 · 则 Q 也是正交矩阵, 且 A = (Hn−1 · · · H2H1) −1R = QR. 以上就是基于 Householder 变换的 QR 分解的具体实现过程. 最后所得到的上三角矩阵 R 就存放 在 A 的上三角部分. b 如果不需要生成 Q, 则基于 Householder 变换的 QR 分解的运算量大约为 2mn2 − 2n 3/3. 矩阵 Q 可通过下面的累积方法来计算: Q = Im, Q = QHk, k = 1, 2, . . . , n. 如果保留了每一步的 Householder 向量, 则 Q 也可以通过下面的向后累积方法来计算: Q = Im, Q = HkQ, k = n, n − 1, . . . , 1. 这样做的好处是一开始 Q 会非常稀疏 (左上角为单位矩阵), 随着迭代的进行, Q 才会慢慢变满. 而 前面的计算方法, 第一步就将 Q 变成了一个满矩阵. 采用向后累积方法计算 Q 的运算量大约为 4m2n − 4mn2 + 4n 3/3. 如果只需要计算 Q 的前 n 列, 则运算量大约为 2mn2 − 2n 3/3, 此时 QR 分解的总运算量为 4mn2 − 4n 3/3. 若 m = n, 则为 8n 3/3. 算法 3.5. 基于 Householder 变换的 QR 分解 % Given A ∈ R m×n , compute Q and R such that A = QR where Q ∈ R m×m and R ∈ R m×n % The upper triangular part of R is stored in the upper triangular part of A 1: Set Q = Im×m 2: for k = 1 to n do 3: x = A(k : m, k) 4: [βk, vk] = House(x) 5: A(k : m, k : n) = (Im−k+1 − βkvkv T k )A(k : m, k : n) 6: = A(k : m, k : n) − βkvk v T k A(k : m, k : n) 7: Q(:, k : m) = Q(:, k : m)(Im−k+1 − βkvkv T k ) 8: = Q(:, k : m) − βk Q(:, k : m)vk v T k 9: end for b 上面的算法只是关于 Householder QR 分解的一个简单描述, 并没有考虑运算量问题. 在实 际计算时, 我们通常会保留所有的 Householder 向量 (用于向后累积法计算 Q). 由于第 k 步 中 H˜ k 所对应的 Householder 向量 vk 的长度为 m − k + 1, 因此我们先把 vk 单位化, 使得 vk 的第一元素为 1, 这样就只要存储 vk(2 : end), 共 m − k 个元素. 于是我们就可以把所有的
-98第三讲线性最小二乘问题Householder向量存放在A的严格下三角部分.而A的上三角部分仍然存放R事实上，如果没有明确要求生成Q，通常只需存储U和Bk，即以因式分解（Factored-Form）形式存储 Q [57].我们也可以考虑块HouseholderOR分解，以便充分利用3级BLAS运算，提高计算效率算法3.5针对的是实矩阵，如果是复矩阵则要适当做些修改3.3.4列主元QR分解如果A不是满秩，记1会rank(A)<n，则存在一个置换矩阵P，使得AP的前1列是线性无关的.因此我们可以对AP进行QR分解，于是我们可以得到下面的结论定理3.10 (列主元 QR 分解) 设 A E Cmxn (m ≥ n),且 rank(A)=l < n. 则存在置换矩阵 P 和正交矩阵QECmxm,使得R11R12AP:其中R11EClxl是非奇异上三角矩阵,且对角线元素满足r11≥r22≥≥ru>0.上述结论也可简化为AP=Q1R1R1其中Q1 ECmxl单位列正交(上述结论中 Q的前1列)列主元OR分解的实现过程与OR分解基本类似，只是在第k步时，需要选列主元，同时可能需要做一个列交换假设经过k-1步后，我们得到下面的分解AP(-1) = Q(k-1) [Rt-1) R(-1)](Q(k-1)AP(k-1) = R(k-1)Q(k-1)R(k-1),即R(-1)0其中P(k-1)是置换矩阵,Q(k-1)是正交矩阵，R(h-1)eR(#-1)x(k-1)是非奇异上三角矩阵.下面考虑第k步：(1)首先计算R-1)的所有列的范数,如果范数都为0,则R(-1)=0,此时必有k-1=1,算法结束.(2)当k≤1时,R(-1)≠0,记其范数最大的列为第i列(如果有相等的,取其中一个即可).若i≠1,则交换R(k-1)的第k列与第谁+k-1列,并记相应的置换矩阵为Ps.(3）由于列交换不会影响到R(k-1)的前k－1列,因此列交换后的矩阵可记为[R(-1 (-1)]R(k-1)Pk≤R6-1]0
· 98 · 第三讲 线性最小二乘问题 Householder 向量存放在 A 的严格下三角部分. 而 A 的上三角部分仍然存放 R. b 事实上, 如果没有明确要求生成 Q, 通常只需存储 vk 和 βk, 即以因式分解 (Factored-Form) 形式存储 Q [57]. b 我们也可以考虑块 Householder QR 分解, 以便充分利用 3 级 BLAS 运算, 提高计算效率. b 算法 3.5 针对的是实矩阵, 如果是复矩阵则要适当做些修改. 3.3.4 列主元 QR 分解 如果 A 不是满秩, 记 l ≜ rank(A) < n, 则存在一个置换矩阵 P, 使得 AP 的前 l 列是线性无关 的. 因此我们可以对 AP 进行 QR 分解, 于是我们可以得到下面的结论. 定理 3.10 (列主元 QR 分解) 设 A ∈ C m×n (m ≥ n), 且 rank(A) = l < n. 则存在置换矩阵 P 和 正交矩阵 Q ∈ C m×m, 使得 AP = Q " R11 R12 0 0 # m×n , 其中 R11 ∈ C l×l 是非奇异上三角矩阵, 且对角线元素满足 r11 ≥ r22 ≥ · · · ≥ rll > 0. b 上述结论也可简化为 AP = Q1 h R11 R12i , 其中 Q1 ∈ C m×l 单位列正交 (上述结论中 Q 的前 l 列). 列主元 QR 分解的实现过程与 QR 分解基本类似, 只是在第 k 步时, 需要选列主元, 同时可能 需要做一个列交换. 假设经过 k − 1 步后, 我们得到下面的分解 AP(k−1) = Q (k−1) " R (k−1) 11 R (k−1) 12 0 R (k−1) 22 # ≜ Q (k−1)R (k−1) , 即 Q (k−1)T AP(k−1) = R (k−1) , 其中 P (k−1) 是置换矩阵, Q(k−1) 是正交矩阵, R (k−1) 11 ∈ R (k−1)×(k−1) 是非奇异上三角矩阵. 下面考虑第 k 步: (1) 首先计算 R (k−1) 22 的所有列的范数, 如果范数都为 0, 则 R (k−1) 22 = 0, 此时必有 k − 1 = l, 算 法结束. (2) 当 k ≤ l 时, R (k−1) 22 ̸= 0, 记其范数最大的列为第 ik 列 (如果有相等的, 取其中一个即可). 若 ik ̸= 1, 则交换 R(k−1) 的第 k 列与第 ik + k − 1 列, 并记相应的置换矩阵为 Pk. (3) 由于列交换不会影响到 R(k−1) 的前 k − 1 列, 因此列交换后的矩阵可记为 R (k−1)Pk ≜ " R (k−1) 11 R˜ (k−1) 12 0 R˜ (k−1) 22 #