:78.第三讲线性最小二乘问题Ui=2r1,β=2/v?7:else8:9:U1 =0,β= 0endif10:11:else%=l212:Vri+aα=13:if <0 then14:=-15:else16:ui=-/(ri+α)endif17:18:β= 2/(u2 +o)19:endif+在实际计算时，我们可以将向量u单位化，使得=1.这样，我们就无需为另外分配空间，而是将u(2:n)存放在r(2:n)中，因为经过Householder变换后，向量除第个分量外,其它都为零.同时，为了避免可能产生的溢出，我们也可以事先将单位化，即令=/川2MATLAB源代码3.1.House向量function [beta,v] = House(x)2n = length(x);3sigma = dot(x(2:n),x(2:n));4V=x;5if sigma ==06if x(1)<07v(1) - 2*x(1);8beta =2/(v(1)*v(1));9else10v(1) = 1;11beta =e;12end13else14alpha = sqrt(x(1)*x(1)+sigma);15if x(1) < 016v(1) =x(1)alpha;17else18v(1)=-sigma/(x(1)+alpha);19end20beta = 2/(v(1)*v(1)+sigma);21end
· 78 · 第三讲 线性最小二乘问题 7: v1 = 2x1, β = 2/v 2 1 8: else 9: v1 = 0, β = 0 10: end if 11: else 12: α = √ x 2 1 + σ % α = ∥x∥2 13: if x1 < 0 then 14: v1 = x1 − α 15: else 16: v1 = −σ/(x1 + α) 17: end if 18: β = 2/(v 2 1 + σ) 19: end if † 在实际计算时, 我们可以将向量 v 单位化, 使得 v1 = 1. 这样, 我们就无需为 v 另外分配空间, 而 是将 v(2 :n) 存放在 x(2 :n) 中, 因为经过 Householder 变换后, 向量 x 除第一个分量外, 其它都为 零. 同时, 为了避免可能产生的溢出, 我们也可以事先将 x 单位化, 即令 x = x/∥x∥2 ✞ MATLAB 源代码 3.1. House 向量 ☎ 1 function [beta,v] = House(x) 2 n = length(x); 3 sigma = dot(x(2:n),x(2:n)); 4 v = x; 5 if sigma == 0 6 if x(1) < 0 7 v(1) = 2*x(1); 8 beta = 2/(v(1)*v(1)); 9 else 10 v(1) = 1; 11 beta = 0; 12 end 13 else 14 alpha = sqrt(x(1)*x(1)+sigma); 15 if x(1) < 0 16 v(1) = x(1) - alpha; 17 else 18 v(1) = -sigma/(x(1)+alpha); 19 end 20 beta = 2/(v(1)*v(1)+sigma); 21 end ✝ ✆
3.2初等变换矩阵. 79 .可以证明，上述算法具有很好的数值稳定性[75]，即H-H|2=O(eu),其中H是由上述算法计算得到的近似矩阵，eu是机器精度Houscholder变换运算量设AERmxn,H=I-Buu*ERm,则HA= (I - βuu*)A= A- Bu(u*A)因此，在做Householder变换时，并不需要生成Householder矩阵，只需要Householder向量即可.上面矩阵相乘的总运算量大约为4mn.3.2.4Givens变换为简单起见，我们这里这讨论实数域中的Givens变换.设e[0,2元]，我们称矩阵5ERnxn,G(i, j,0) =(i<j)为Givens变换(或Givens旋转，或Givens矩阵),其中c=cos(0)，s=sin(0).即将单位矩阵的(i,i）和(j.i）位置上的元素用c代替，而（i,）和（G,i）位置上的元素分别用s和一s代替，所得到的矩阵就是G(i,j,0).定理3.5G(i,j,0)是正交矩阵，且det(G(i,j,0))=1我们需要注意的是，当一个矩阵左乘一个Givens矩阵时，只会影响其第i行和第i行的元素.而当个矩阵右乘一个Givens矩阵时，只会影响其第i和第i列的元素ER2x2使得GaER2,则存在一个Givens变换G其中c,s和例3.1设0T2r的值如下：若1=2=0,则c=1,S=0,r=0;若1=0但20,则c=08=2//2，=2；若1±0但2=0,则c=sign(r1),s=0,r=ril;.若≠0且2≠0,则c=/r，8=2/r,=V+也就是说，通过Givens变换，我们可以将向量工ER2的第二个分量化为0事实上，对于任意一个向量rERn，我们都可以通过Givens变换将其任意一个位置上的分量化为
3.2 初等变换矩阵 · 79 · 可以证明, 上述算法具有很好的数值稳定性 [75], 即 ∥H˜ − H∥2 = O(εu), 其中 H˜ 是由上述算法计算得到的近似矩阵, εu 是机器精度. Householder 变换运算量 设 A ∈ R m×n, H = I − βvv∗ ∈ R m, 则 HA = (I − βvv∗ )A = A − βv(v ∗A). 因此, 在做 Householder 变换时, 并不需要生成 Householder 矩阵, 只需要 Householder 向量即可. 上面矩 阵相乘的总运算量大约为 4mn. 3.2.4 Givens 变换 为简单起见, 我们这里这讨论实数域中的 Givens 变换. 设 θ ∈ [0, 2π], 我们称矩阵 G(i, j, θ) = 1 . . . c s . . . −s c . . . 1 ∈ R n×n , (i ≤ j) 为 Givens 变换 (或 Givens 旋转, 或 Givens 矩阵), 其中 c = cos(θ), s = sin(θ). 即将单位矩阵的 (i, i) 和 (j, j) 位置上的元素用 c 代替, 而 (i, j) 和 (j, i) 位置上的元素分别用 s 和 −s 代替, 所得到的矩阵就是 G(i, j, θ). 定理 3.5 G(i, j, θ) 是正交矩阵, 且 det(G(i, j, θ)) = 1. 我们需要注意的是, 当一个矩阵左乘一个 Givens 矩阵时, 只会影响其第 i 行和第 j 行的元素. 而当 一个矩阵右乘一个 Givens 矩阵时, 只会影响其第 i 和第 j 列的元素. 例 3.1 设 x = [ x1 x2 ] ∈ R 2 , 则存在一个 Givens 变换 G = [ c s −s c] ∈ R 2×2 使得 Gx = [ r 0 ] , 其中 c, s 和 r 的值如下: • 若 x1 = x2 = 0, 则 c = 1, s = 0, r = 0 ; • 若 x1 = 0 但 x2 ̸= 0, 则 c = 0, s = x2/|x2|, r = |x2| ; • 若 x1 ̸= 0 但 x2 = 0, 则 c = sign(x1), s = 0, r = |x1| ; • 若 x1 ̸= 0 且 x2 ̸= 0, 则 c = x1/r, s = x2/r, r = √ x 2 1 + x 2 2 . 也就是说, 通过 Givens 变换, 我们可以将向量 x ∈ R 2 的第二个分量化为 0. 事实上, 对于任意一个向量 x ∈ R n, 我们都可以通过 Givens 变换将其任意一个位置上的分量化为
第三讲线性最小二乘问题.80.0.更进一步，我们也可以通过若干个Givens变换，将中除第一个分量外的所有元素都化为0.算法3.2.Givens变换% Given =[a,b]T R?, compute c and s such that Ga =[r,O]T wherer =r21: function [c, s] -givens(a,b)2: if b =0 then3:if a≥ 0 then4:c=l,s=05:else6:S=0-C=endif7:8:else9:if [6] > [a| thensign(b)a10:?ST6'VI+T2else11:6sign(a)12:T=CTVI+T2aend if13:14:endif3.2.5正交变换的舍人误差分析引理3.2设PERnxn是一个精确的Householder或Givens变换，P是其浮点运算近似，则f(PA)=P(A+E),f(AP) =(A+F)P其中E|2=O(u)·A2,F2=O(u)-A|2.这说明对一个矩阵做Householder变换或Givens变换是向后稳定的定理3.6考虑对矩阵A做一系列的正交变换，则有f(P...PAQ....Q)=P....P(A+E)Q...Qk,其中E2=O(eu)·(A|2).这说明整个计算过程是向后稳定的一般地，假设X是一个非奇异的线性变换，X是其浮点运算近似.当X作用到A上时，我们有f(XA)=XA+E=X(A+X-1E)≤ X(A+F),其中Ell2=O(u)-XAll2≤O(eu)·X|l2·/All2,故F2=X-1E2≤O(eu) -1X-12-IXI2-IAll2=O(u) -K2(X)-IIA2因此，舍入误差将被放大K2(X）倍.当X是正交变换时，k2(X）达到最小值1，这就是为什么在浮点运算中尽量使用正交变换的原因
· 80 · 第三讲 线性最小二乘问题 0. 更进一步, 我们也可以通过若干个 Givens 变换, 将 x 中除第一个分量外的所有元素都化为 0. 算法 3.2. Givens 变换 % Given x = [a, b] ⊺ ∈ R 2 , compute c and s such that Gx = [r, 0]⊺ where r = ∥x∥2 1: function [c, s] = givens(a, b) 2: if b = 0 then 3: if a ≥ 0 then 4: c = 1, s = 0 5: else 6: c = −1, s = 0 7: end if 8: else 9: if |b| > |a| then 10: τ = a b , s = sign(b) √ 1 + τ 2 , c = sτ 11: else 12: τ = b a , c = sign(a) √ 1 + τ 2 , s = cτ 13: end if 14: end if 3.2.5 正交变换的舍入误差分析 引理 3.2 设 P ∈ R n×n 是一个精确的 Householder 或 Givens 变换, P˜ 是其浮点运算近似, 则 fl(P A˜ ) = P(A + E), fl(AP˜) = (A + F)P, 其中 ∥E∥2 = O(εu) · ∥A∥2, ∥F∥2 = O(εu) · ∥A∥2. 这说明对一个矩阵做 Householder 变换或 Givens 变换是向后稳定的. 定理 3.6 考虑对矩阵 A 做一系列的正交变换, 则有 fl(P˜ k · · · P˜ 1AQ˜ 1 · · · Q˜ k) = Pk · · · P1(A + E)Q1 · · · Qk, 其中 ∥E∥2 = O(εu) · (k∥A∥2). 这说明整个计算过程是向后稳定的. 一般地, 假设 X 是一个非奇异的线性变换, X˜ 是其浮点运算近似. 当 X 作用到 A 上时, 我们有 fl(XA˜ ) = XA + E = X(A + X−1E) ≜ X(A + F), 其中 ∥E∥2 = O(εu) · ∥XA∥2 ≤ O(εu) · ∥X∥2 · ∥A∥2, 故 ∥F∥2 = ∥X−1E∥2 ≤ O(εu) · ∥X−1 ∥2 · ∥X∥2 · ∥A∥2 = O(εu) · κ2(X) · ∥A∥2, 因此, 舍入误差将被放大 κ2(X) 倍. 当 X 是正交变换时, κ2(X) 达到最小值 1, 这就是为什么在浮点运 算中尽量使用正交变换的原因
3.3QR分解·81.3.3QR分解QR分解是将一个矩阵分解一个正交矩阵（酉矩阵）和一个三角矩阵的乘积.QR分解被广泛应用于线性最小二乘问题的求解和矩阵特征值的计算，3.3.1QR分解的存在性与唯一性定理3.7(QR分解）[43]设A ECmxn(m≥n).则存在一个单位列正交矩阵QECmxn(即Q*QInxn）和一个上三角矩阵ReCnxn，使得(3.4)A=QR.若A列满秩，则存在一个具有正对角线元素的上三角矩阵R使得（3.4)成立，且此时QR分解唯一，即Q和R都唯一证明.设A=[a1,a2，...an)ECmxn.若A列满秩，即rank(A)=n.则QR分解(3.4)就是对A的列向量组进行Gram-Schmidt正交化过程的矩阵描述（见算法3.3)算法3.3.Gram-Schmidt过程1: r11 = a122:qi=a1/r113: for j=2 to n do4:qj = aj5:fori=lto j-ldo6:%9表示共轭转置rij=q,aj7:qi=q-rqiend for8:9:Tjj = Ilq;ll210:qj = q;/rijIl:end for如果A不是列满秩，我们可以做类似的正交化过程：·如果a1=0,则令q1=0;否则令q1=a1/la1l2;对于j=23,..,n计算=aj-=(qa)gi如果=0,则表明a可以由a1,2,j-1线性表出，令9=0.否则令9;=;/l;ll2于是我们有A=QR,其中Q=[91,92….，9n]列正交(但不是单位列正交),其列向量要么是单位向量，要么就是零向量.这里的R=[ri]nxn是上三角矩阵,定义如下fori<jqiaj,Ti0.fori>j
3.3 QR 分解 · 81 · 3.3 QR 分解 QR 分解是将一个矩阵分解一个正交矩阵 (酉矩阵) 和一个三角矩阵的乘积. QR 分解被广泛应用于 线性最小二乘问题的求解和矩阵特征值的计算. 3.3.1 QR 分解的存在性与唯一性 定理 3.7 (QR 分解) [43] 设 A ∈ C m×n (m ≥ n). 则存在一个单位列正交矩阵 Q ∈ C m×n (即 Q∗Q = In×n) 和一个上三角矩阵 R ∈ C n×n, 使得 A = QR. (3.4) 若 A 列满秩, 则存在一个具有正对角线元素的上三角矩阵 R 使得 (3.4) 成立, 且此时 QR 分解唯一, 即 Q 和 R 都唯一. 证明. 设 A = [a1, a2, . . . , an] ∈ C m×n. 若 A 列满秩, 即 rank(A) = n. 则 QR 分解 (3.4) 就是对 A 的列向 量组进行 Gram-Schmidt 正交化过程的矩阵描述 (见算法 3.3). 算法 3.3. Gram-Schmidt 过程 1: r11 = ∥a1∥2 2: q1 = a1/r11 3: for j = 2 to n do 4: qj = aj 5: for i = 1 to j − 1 do 6: rij = q ∗ i aj % q ∗ i 表示共轭转置 7: qj = qj − rij qi 8: end for 9: rjj = ∥qj∥2 10: qj = qj/rjj 11: end for 如果 A 不是列满秩, 我们可以做类似的正交化过程: • 如果 a1 = 0, 则令 q1 = 0; 否则令 q1 = a1/∥a1∥2 ; • 对于 j = 2, 3, . . . , n, 计算 q˜j = aj − ∑j−1 i=1 (q ∗ i aj )qi . 如果 q˜j = 0, 则表明 aj 可以由 a1, a2, . . . , aj−1 线性表出, 令 qj = 0. 否则令 qj = ˜qj/∥q˜j∥2 . 于是我们有 A = QR, 其中 Q = [q1, q2, . . . , qn] 列正交 (但不是单位列正交), 其列向量要么是单位向量, 要么就是零向量. 这里 的 R = [rij ]n×n 是上三角矩阵, 定义如下 rij = q ∗ i aj , for i ≤ j 0, for i > j
第三讲线性最小二乘问题.82.如果Q的某一列9k=0,则R中对应的第k行就全为0设rank(A)=1<n,则Q有1个非零列，设为qi,qi2...,qi它们形成Cm中的一个单位正交向量组，所以我们可以将其扩展成Cm中的一组标准正交基，即qii,qia....,qiu,q1,...,qml然后我们用91替换Q中的第一个零列，用92替换Q中的第二个零列，依此类推，将Q中的所有零列都替换掉.将最后得到的矩阵记为Q,则QeCmxn单位列正交，且QR=QR这是由于Q中的新添加的列向量正好与R中的零行相对应.所以我们有OR分解A=QR.下面证明满秩矩阵QR分解的存在唯一性存在性：由于A列满秩，由Gram-Schmidt正交化过程(算法3.3)可知,存在上三角矩阵R=[ri]nxn满足rji>0,使得A=QR,其中Q单位列正交唯一性：假设A存在QR分解A=Q,Ri=Q2R2,其中Q1,Q2ECmxn单位列正交，R1,R2ECnxn为具有正对角元素的上三角矩阵.则有Q1= Q2R2R-1.(3.5)由于R1,R2均为上三角矩阵，所以R2R-1也是上三角矩阵，且其对角线元素为R2(i,)/R1(i,i)，=1,2....n.由(3.5)可得1=Q1ll2=IQ2R2R-1l2=R2R-112所以R2(i,i)≤1,i=1,2,...,n.Ri(i,i)同理可证Ri(i,)/R2(i,）≤1.所以Ri(i,i) = R2(i,i), i=l,2,...,n.又IlQi=tr(QiQ1）=n，所以由(3.5)可知R2R-:=Q2R2R-=Q1/=n由于R2R-1的对角线元素都是1,所以R2R-1只能是单位矩阵，即R2=R1.因此Q2=AR-1=口AR-1=Q1,即A的QR分解是唯一的+有时也将QR分解定义为：存在西矩阵QECmxm使得Ri1A=Q
· 82 · 第三讲 线性最小二乘问题 如果 Q 的某一列 qk = 0, 则 R 中对应的第 k 行就全为 0. 设 rank(A) = l < n, 则 Q 有 l 个非零列, 设为 qi1 , qi2 , . . . , qil . 它们形成 C m 中的一个单位正交向量 组, 所以我们可以将其扩展成 C m 中的一组标准正交基, 即 qi1 , qi2 , . . . , qil , q˜1, . . . , q˜m−l . 然后我们用 q˜1 替换 Q 中的第一个零列, 用 q˜2 替换 Q 中的第二个零列, 依此类推, 将 Q 中的所有零列都 替换掉. 将最后得到的矩阵记为 Q˜, 则 Q˜ ∈ C m×n 单位列正交, 且 QR˜ = QR. 这是由于 Q˜ 中的新添加的列向量正好与 R 中的零行相对应. 所以我们有 QR 分解 A = QR. ˜ 下面证明满秩矩阵 QR 分解的存在唯一性. 存在性: 由于 A 列满秩, 由 Gram-Schmidt 正交化过程 (算法 3.3) 可知, 存在上三角矩阵 R = [rij ]n×n 满足 rjj > 0, 使得 A = QR, 其中 Q 单位列正交. 唯一性: 假设 A 存在 QR 分解 A = Q1R1 = Q2R2, 其中 Q1, Q2 ∈ C m×n 单位列正交, R1, R2 ∈ C n×n 为具有正对角元素的上三角矩阵. 则有 Q1 = Q2R2R −1 1 . (3.5) 由于 R1, R2 均为上三角矩阵, 所以 R2R −1 1 也是上三角矩阵, 且其对角线元素为 R2(i, i)/R1(i, i), i = 1, 2, . . . , n. 由 (3.5) 可得 1 = ∥Q1∥2 = ∥Q2R2R −1 1 ∥2 = ∥R2R −1 1 ∥2. 所以 R2(i, i) R1(i, i) ≤ 1, i = 1, 2, . . . , n. 同理可证 R1(i, i)/R2(i, i) ≤ 1. 所以 R1(i, i) = R2(i, i), i = 1, 2, . . . , n. 又 ∥Q1∥ 2 F = tr(Q∗ 1Q1) = n, 所以由 (3.5) 可知 ∥R2R −1 1 ∥ 2 F = ∥Q2R2R −1 1 ∥ 2 F = ∥Q1∥ 2 F = n. 由于 R2R −1 1 的对角线元素都是 1, 所以 R2R −1 1 只能是单位矩阵, 即 R2 = R1. 因此 Q2 = AR−1 2 = AR−1 1 = Q1, 即 A 的 QR 分解是唯一的. □ † 有时也将 QR 分解定义为: 存在酉矩阵 Q ∈ C m×m 使得 A = Q [ R11 0 ]