《矩阵计算》课程教学资源（讲义）第三讲 线性最小二乘问题

3.2几类重要的矩阵变换·87.Householder矩阵的一个非常重要的应用就是可以将一个向量除第一个元素以外的所有元素都化为零.我们首先给出一个引理引理3.4 设 a,y E Cn 为任意两个互异的向量,则存在一个 Householder矩阵 H(u)使得 y =(板书)H(u)的充要条件是lll2=llyll2且a*yER证明.若all2=yll2且ryER，则y*y=r*a且a*y=y*a.于是Ilr -yllz=(r -y)*(α -y) =a*a -y*r -r*y+yy=2(a*a -y*a).令=-y则有2(-y)(r-y)*r2(r-y)(r-yr)H(v)r=rya-yl22(r*r -y*r)即存在Householder矩阵H(u)使得y=H(u)r.反之,如果存在Householder矩阵H使得y=Hr,由于H是Hermite的，所以r*y=r*HaER.又因为H是酉矩阵，所以Ilyll2=Hrll2=cl2.Householder向量也可以取=ei（-y),0≤?<2元.口如果a,y都是实向量，则条件r*yER自然成立，此时充要条件就是rll2=llyll2由引理3.4,我们可以立即得到下面的结论.定理3.5设=[1,2,..,nJTERn是一个非零向量，则存在Householder矩阵H()使得H(u) =Qei,其中 α=ll2 (或 α =-[ll2),e1 =[1,0,...,o]T e Rn在后面的讨论中,我们将定理中的向量称为对应的Householder向量设r=[1,2.….,n]TERn是一个实的非零向量，下面讨论如何计算定理3.5中House-holder矩阵H(u）所对应的Householder向量w.由引理3.4的证明过程可知u= a -Qei = [e1 - a, x2,..,an]T.在实际计算中，为了尽可能地减少舍入误差，我们通常避免两个相近的数做减运算，否则就会损失有效数字.因此，我通常取(3.6)α =- sign(r1) lll2.事实上，我们也可以取α=sign(r1)lll2，但此时为了减少舍人误差，我们需要通过下面的公式来计算的第一个分量u1α=sign()2, = - α=±-_(+晚++2)(3.7)T1+QTi +Q

3.2 几类重要的矩阵变换 · 87 · Householder 矩阵的一个非常重要的应用就是可以将一个向量除第一个元素以外的所有元素 都化为零. 我们首先给出一个引理. 引理 3.4 设 x, y ∈ C n 为任意两个互异的向量, 则存在一个 Householder 矩阵 H(v) 使得 y = H(v)x 的充要条件是 ∥x∥2 = ∥y∥2 且 x ∗y ∈ R. (板书) 证明. 若 ∥x∥2 = ∥y∥2 且 x ∗y ∈ R, 则 y ∗y = x ∗x 且 x ∗y = y ∗x. 于是 ∥x − y∥ 2 2 = (x − y) ∗ (x − y) = x ∗x − y ∗x − x ∗ y + y ∗ y = 2(x ∗x − y ∗x). 令 v = x − y, 则有 H(v)x = x − 2(x − y)(x − y) ∗x ∥x − y∥ 2 2 = x − 2(x − y)(x ∗x − y ∗x) 2(x ∗x − y ∗x) = y, 即存在 Householder 矩阵 H(v) 使得 y = H(v)x. 反之, 如果存在 Householder 矩阵 H 使得 y = Hx, 由于 H 是 Hermite 的, 所以 x ∗y = x ∗Hx ∈ R. 又因为 H 是酉矩阵, 所以 ∥y∥2 = ∥Hx∥2 = ∥x∥2. b Householder 向量也可以取 v = e i θ (x − y), 0 ≤ θ < 2π. □ b 如果 x, y 都是实向量, 则条件 x ∗y ∈ R 自然成立, 此时充要条件就是 ∥x∥2 = ∥y∥2. 由引理 3.4, 我们可以立即得到下面的结论. 定理 3.5 设 x = [x1, x2, . . . , xn] T ∈ R n 是一个非零向量, 则存在 Householder 矩阵 H(v) 使得 H(v)x = αe1, 其中 α = ∥x∥2 (或 α = −∥x∥2), e1 = [1, 0, . . . , 0]T ∈ R n . b 在后面的讨论中, 我们将定理中的向量 v 称为 x 对应的 Householder 向量. 设 x = [x1, x2, . . . , xn] T ∈ R n 是一个实的非零向量, 下面讨论如何计算定理 3.5 中 House￾holder 矩阵 H(v) 所对应的 Householder 向量 v. 由引理 3.4 的证明过程可知 v = x − αe1 = [x1 − α, x2, . . . , xn] T . 在实际计算中, 为了尽可能地减少舍入误差, 我们通常避免两个相近的数做减运算, 否则就会损失 有效数字. 因此, 我通常取 α = − sign(x1) · ∥x∥2. (3.6) 事实上, 我们也可以取 α = sign(x1)∥x∥2, 但此时为了减少舍入误差, 我们需要通过下面的公式来 计算 v 的第一个分量 v1 α = sign(x1)∥x∥2, v1 = x1 − α = x 2 1 − ∥x∥ 2 2 x1 + α = −(x 2 2 + x 2 3 + · · · + x 2 n ) x1 + α . (3.7)

-88.第三讲线性最小二乘问题在U1的两种计算方法(3.6)和（3.7)中，α的取值都与1的符号有关.但在某些应用中，我们需要确保α非负，此时我们可以将这两种方法结合起来使用，即：if sign(r1) < 0T1- α,U1(+++)otherwise1+α无论怎样选取α我们都有H=I-βBv*，其中2221B=2220r1U*aU1(i-α)2 +r +..+r2思考：如果是复向量，有没有相应的结论在实数域中计算Householder向量的算法如下，总运算量大约为2n，如果加上赋值运算的话则为3n.算法3.1.计算Householder向量% Given r e R"n, compute w E Rn such that Hr = lrllze1, where H = I - Buu*(House.m)1: function [β, u] -House(c)2: n =length(α)(here length(r) denotes the dimension of c)3: compute o = r +a + ...+r and let u =4: if α=0 and 21>=0 thenβ=0 %U=05:6:elseif g=0and i<0 then7:U1 = 2r1,β= 2/u28: else9:Q=V+% Q= l210:if 1<= 0 then11: = 1 - α12:else13:Ui=-o/(ri+α)end if14:15:β=2/(u+o)16: end if可以证明，上述算法具有很好的数值稳定性[135]，即H - Hl2 = O(eu),其中H是由上述算法计算得到的近似矩阵，eu是机器精度

· 88 · 第三讲 线性最小二乘问题 在 v1 的两种计算方法 (3.6) 和 (3.7) 中, α 的取值都与 x1 的符号有关. 但在某些应用中, 我们需 要确保 α 非负, 此时我们可以将这两种方法结合起来使用, 即: v1 =    x1 − α, if sign(x1) < 0 −(x 2 2 + x 2 3 + · · · + x 2 n ) x1 + α , otherwise 无论怎样选取 α, 我们都有 H = I − βvv∗ , 其中 β = 2 v ∗v = 2 (x1 − α) 2 + x 2 2 + · · · + x 2 n = 2 2α2 − 2αx1 = − 1 αv1 . 思考：如果 x 是复向量, 有没有相应的结论? 在实数域中计算 Householder 向量 v 的算法如下, 总运算量大约为 2n, 如果加上赋值运算的 话则为 3n. 算法 3.1. 计算 Householder 向量 % Given x ∈ R n , compute v ∈ R n such that Hx = ∥x∥2e1, where H = I − βvv∗ (House.m) 1: function [β, v] = House(x) 2: n = length(x) (here length(x) denotes the dimension of x) 3: compute σ = x 2 2 + x 2 3 + · · · + x 2 n and let v = x 4: if σ = 0 and x1 >= 0 then 5: β = 0 % v1 = 0 6: else if σ = 0 and x1 < 0 then 7: v1 = 2x1, β = 2/v 2 1 8: else 9: α = p x 2 1 + σ % α = ∥x∥2 10: if x1 <= 0 then 11: v1 = x1 − α 12: else 13: v1 = −σ/(x1 + α) 14: end if 15: β = 2/(v 2 1 + σ) 16: end if 可以证明, 上述算法具有很好的数值稳定性 [135], 即 ∥H˜ − H∥2 = O(εu), 其中 H˜ 是由上述算法计算得到的近似矩阵, εu 是机器精度

3.2几类重要的矩阵变换·89.在实际计算时，我们可以将向量单位化，使得1=1.这样，我们就无需为另外分配空间，而是将u(2：n）存放在r(2：n）中，因为经过Householder变换后，向量r除第一个分量外，其它都为零思考：这里要求U1≠0,那么什么情况下U1=0？为了避免可能产生的溢出，我们也可以事先将单位化，即令=/2Householder变换与矩阵的乘积设AeRmxn,H=I-Buu*ERmxm,则HA=(I-Bvu*)A=A-βv(u*A)因此，在做Householder变换时，并不需要生成Householder矩阵，只需要Householder向量即可.上面矩阵相乘的总运算量大约为4mn.注记·Householder矩阵早在1932年就出现了[128],但直到1958年才由Householder用于矩阵的三角化[75].·关于其他类型的Householder变换（包括复数情形)，可以参见[36]3.2.4Givens变换为简单起见，我们这里这讨论实数域中的Givens变换.设e[0,2元],我们称矩阵1SERnxnG(i, j, 0) =(i≤).1为Givens变换(或Givens旋转，或Givens矩阵),其中c=cos(0),8=sin(0).也就是说，将单位矩阵的（i,i）和（i.i）位置上的元素用c代替，而（i,）和（j,i）位置上的元素分别用s和-s代替，所得到的矩阵就是G(i,j,の).定理3.6 G(i,j,0)是正交矩阵,且det(G(i,j,0))=1.Givens变换不是初等矩阵变换,事实上，它是单位矩阵的一个秩2修正，即[ei,e]TG(i, j, 0) = I + [ei,ej如果是定义在复数域上的Givens变换，则c=eiacos0,s=eipsin0,0≤α,β,0<2

3.2 几类重要的矩阵变换 · 89 · b 在实际计算时, 我们可以将向量 v 单位化, 使得 v1 = 1. 这样, 我们就无需为 v 另外分配空 间, 而是将 v(2 : n) 存放在 x(2 : n) 中, 因为经过 Householder 变换后, 向量 x 除第一个分量 外, 其它都为零. 思考：这里要求 v1 ̸= 0, 那么什么情况下 v1 = 0? b 为了避免可能产生的溢出, 我们也可以事先将 x 单位化, 即令 x = x/∥x∥2 Householder 变换与矩阵的乘积 设 A ∈ R m×n , H = I − βvv∗ ∈ R m×m, 则 HA = (I − βvv∗ )A = A − βv(v ∗A). 因此, 在做 Householder 变换时, 并不需要生成 Householder 矩阵, 只需要 Householder 向量即可. 上 面矩阵相乘的总运算量大约为 4mn. 注记 • Householder 矩阵早在 1932 年就出现了 [128], 但直到 1958 年才由 Householder 用于矩阵 的三角化 [75]. • 关于其他类型的 Householder 变换 (包括复数情形), 可以参见 [36]. 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.6 G(i, j, θ) 是正交矩阵, 且 det(G(i, j, θ)) = 1. b Givens 变换不是初等矩阵变换, 事实上, 它是单位矩阵的一个秩 2 修正, 即 G(i, j, θ) = I + [ei , ej ] " c − 1 s −s c − 1 # [ei , ej ] T . b 如果是定义在复数域上的 Givens 变换, 则 c = e i α cos θ, s = e i β sin θ, 0 ≤ α, β, θ < 2π

-90第三讲线性最小二乘问题21E IR2x2 使得 Ga:其中ER2,则存在一个Givens变换G例3.1设c,s和r的值如下：38-22+，即GC=也就是说，通过Givens变换，我们可以将向量ER?的第二个分量化为0.事实上，对任意一个向量ER"，都可以通过Givens变换将其任意一个位置上的分量化为0.更进一步，我们也可以通过若干个Givens变换，将r中除第一个分量外的所有元素都化为0算法3.2.Givens变换% Given r = [ri, r2]T e R?, compute c and s such that Gr = [r, o]T where r = r2 (Givens .m)1: function [c, s] = givens(1, 2)2: if 2=0then3:c=sign(1), 8=04:else5:if |2] >[1] thensign(r2)216:T=ST2V1+ +27:elsesign(r1)328:V1+T2C1end if9:10: end if例3.2(Givens变换的几何意义)T2,则ER2,其中r=Vi+2,α=arctanTir(coscosα-sinsina)rcos(α+0)Gr(sincosa+cossina)r sin(Q +0)sinrsinacos6也就是说，GTa相当于将向量按逆时针旋转Q角度.因此当6=-α时，Ga：

· 90 · 第三讲 线性最小二乘问题 例 3.1 设 x = " x1 x2 # ∈ R 2 , 则存在一个 Givens 变换 G = " c s −s c# ∈ R 2×2 使得 Gx = " r 0 # , 其中 c, s 和 r 的值如下: c = x1 r , s = x2 r , r = q x 2 1 + x 2 2 , 即 G = 1 r " x1 x2 −x2 x1 # . 也就是说, 通过 Givens 变换, 我们可以将向量 x ∈ R 2 的第二个分量化为 0. 事实上, 对任意一个向量 x ∈ R n , 都可以通过 Givens 变换将其任意一个位置上的分量化为 0. 更进一步, 我们也可以通过若干个 Givens 变换, 将 x 中除第一个分量外的所有元素都化为 0. 算法 3.2. Givens 变换 % Given x = [x1, x2] T ∈ R 2 , compute c and s such that Gx = [r, 0]T where r = ∥x∥2 (Givens.m) 1: function [c, s] = givens(x1, x2) 2: if x2 = 0 then 3: c = sign(x1), s = 0 4: else 5: if |x2| > |x1| then 6: τ = x1 x2 , s = sign(x2) √ 1 + τ 2 , c = sτ 7: else 8: τ = x2 x1 , c = sign(x1) √ 1 + τ 2 , s = cτ 9: end if 10: end if 例 3.2 (Givens 变换的几何意义) 设 x = " x1 x2 # = " r cos α r sin α # ∈ R 2 , 其中 r = p x 2 1 + x 2 2 , α = arctan x2 x1 , 则 G Tx = " cos θ sin θ − sin θ cos θ #T " r cos α r sin α # = " r(cos θ cos α − sin θ sin α) r(sin θ cos α + cos θ sin α) # = " r cos(α + θ) r sin(α + θ) # , 也就是说, GTx 相当于将向量 x 按逆时针旋转 θ 角度. 因此当 θ = −α 时, Gx = " r 0 #

3.2几类重要的矩阵变换.91.Givens变换与矩阵的乘积设AeRmxn,G=G(i，j,の)eRm，则GTA只会影响其第i行和第行的元素，也就是说，只需对A的第i行和第行做Givens变换即可，运算量大约为6n.同样地，如果是右乘Givens变换，则只会影响其第i列和第i列的元素，运算量约为6m.任何一个正交矩阵都可以写成若干个Householder矩阵或Givens矩阵的乘积（见习题3.5和3.23)，所以正交矩阵所对应的线性变换可以看作是反射变换和旋转变换的推广，因此它不会改变向量的长度与(不同向量之间的)角度.3.2.5正交变换的舍人误差分析引理3.7设PERnxn是一个精确的Householder或Givens变换，P是其浮点运算近似，则f(PA) =P(A+E),f(AP) =(A+F)P,其中Ell2=O(eu）1A2,Fll2=O(e）·All2这说明对一个矩阵做Householder变换或Givens变换是向后稳定的定理3.8考虑对矩阵A做一系列的正交变换，则有f(P.... PAQ1..Q) = P...P(A+ E)Q1-Qk其中IEll2=O(eu)·(kilAll2).这说明整个计算过程是向后稳定的一般地，假设X是一个非奇异的线性变换，X是其浮点运算近似.当X作用到A上时，有f(XA) =XA+E= X(A+X-1E) X(A+F),其中E|l2=O(eu)XAl2≤O(eu)X|l2/A2,故IFll2=IX-1Ell2≤O(eu)X-1l21X/2llAll2=O(eu)2(X)lAl2因此，舍入误差可能会放大Kk2(X）倍.当X是正交变换时，K2(X）达到最小值1，这就是为什么在浮点运算中尽量使用正交变换的原因

3.2 几类重要的矩阵变换 · 91 · Givens 变换与矩阵的乘积 设 A ∈ R m×n , G = G(i, j, θ) ∈ R m, 则 GTA 只会影响其第 i 行和第 j 行的元素, 也就是说, 只 需对 A 的第 i 行和第 j 行做 Givens 变换即可, 运算量大约为 6n. 同样地, 如果是右乘 Givens 变换, 则只会影响其第 i 列和第 j 列的元素, 运算量约为 6m. b 任何一个正交矩阵都可以写成若干个 Householder 矩阵或 Givens 矩阵的乘积 (见习题 3.5 和 3.23), 所以正交矩阵所对应的线性变换可以看作是反射变换和旋转变换的推广, 因此它不 会改变向量的长度与 (不同向量之间的) 角度. 3.2.5 正交变换的舍入误差分析 引理 3.7 设 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.8 考虑对矩阵 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, 这就是为什么在 浮点运算中尽量使用正交变换的原因