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

第三讲线性最小二乘问题最小二乘问题（LeastSquares）包括线性最小二乘问题，总体最小二乘问题，等式约束最小二乘问题，刚性加权最小二乘问题等等.它在统计学，最优化问题，材料与结构力学，信号与图像处理等方面都有着广泛的应用，是计算数学的一个重要研究分支，也是一个活跃的研究领域本讲主要介绍求解线性最小二乘问题的三种常用方法：正规方程法（也称法方程法），QR分解法和SVD分解法，一般来说，正规方程法是最快的，特别是当A的条件数较小时，正规方程法几乎与其他方法一样精确.而SVD分解法是最慢的，但结果最可靠为了方便起见，我们有时将线性最小二乘问题（3.1）简称为最小二乘问题关于最小二乘问题的相关参考资料AkeBjorck,Numerical MethodsforLeastSquaresProblems,1996[16]魏木生，李莹，赵建立，广义最小二乘问题的理论与计算（第二版),2020[149]3.1问题介绍考虑线性最小二乘问题min I/Ar - bll2(3.1)其中ARmxn，beRm.问题(3.1)的解称为最小二乘解当m=n且A非奇异时，这就是一个线性方程组，解为r*=A-1b;·当m>n时，约束个数大于未知量个数，此时我们称问题（3.1)为超定的·当mn时，线性方程组Ar=b的解可能不存在.此时一般考虑求解最小二乘问题（3.1).记J(r)≤Ar-b易知J）是关于的二次函数，而且是凸函数（当A满秩时，J(a）的Hessen阵是正定的).因此由凸函数的性质可知，是问题(3.1)的解当且仅当是J(）的稳定点。令其一阶导数为零，可得AT AT - ATb = 0.于是将最小二乘问题转化为一个线性方程组，这就是后面的正规方程82

第三讲 线性最小二乘问题 最小二乘问题 (Least Squares ) 包括线性最小二乘问题, 总体最小二乘问题, 等式约束最小二 乘问题, 刚性加权最小二乘问题等等. 它在统计学, 最优化问题, 材料与结构力学, 信号与图像处理 等方面都有着广泛的应用, 是计算数学的一个重要研究分支, 也是一个活跃的研究领域. 本讲主要介绍求解线性最小二乘问题的三种常用方法: 正规方程法 (也称法方程法), QR 分解 法和 SVD 分解法. 一般来说, 正规方程法是最快的, 特别是当 A 的条件数较小时, 正规方程法几 乎与其他方法一样精确. 而 SVD 分解法是最慢的, 但结果最可靠. 为了方便起见, 我们有时将线性最小二乘问题 (3.1) 简称为最小二乘问题. 关于最小二乘问题的相关参考资料 ▶ Åke Björck, Numerical Methods for Least Squares Problems, 1996 [16] ▶ 魏木生, 李莹, 赵建立, 广义最小二乘问题的理论与计算 (第二版), 2020 [149] 3.1 问题介绍 考虑线性最小二乘问题 min x∈Rn ∥Ax − b∥ 2 2 , (3.1) 其中 A ∈ R m×n , b ∈ R m. 问题 (3.1) 的解称为最小二乘解. • 当 m = n 且 A 非奇异时, 这就是一个线性方程组, 解为 x∗ = A−1 b; • 当 m > n 时, 约束个数大于未知量个数, 此时我们称问题 (3.1) 为超定的; • 当 m n 时, 线性方程组 Ax = b 的解可能不存在. 此时一般考虑求解最小二乘问题 (3.1). 记 J(x) ≜ ∥Ax − b∥ 2 2 . 易知 J(x) 是关于 x 的二次函数, 而且是凸函数 (当 A 满秩时, J(x) 的 Hessen 阵是正定的). 因此, 由凸函数的性质可知, x∗ 是问题 (3.1) 的解当且仅当 x∗ 是 J(x) 的稳定点. 令其一阶导数为零, 可 得 A TAx − A T b = 0. 于是将最小二乘问题转化为一个线性方程组, 这就是后面的正规方程. 82

3.1问题介绍·83.如果 A不是满秩,则 ATA半正定,此时解不唯一3.1.2欠定方程组若m<n，则线性方程组Ar=b存在无穷多个解（假定A满秩）.这时我们通常寻求最小范数解，即所有解中范数最小的解.于是原问题就转化为下面的约束优化问题(3.2)对应的Lagrange函数为C(a, ) =,llal2 + T(Ar - b),其中入=[△1,入2,….,入m]T是Lagrange乘子.此时优化问题(3.2)的解就是C(α,)的鞍点，即下面方程组的解：acac=T+AT=0,=Ar-b=0.drx写成矩阵形式为[ -如果A满秩，即rank(A)=m，则系数矩阵非奇异（见练习3.1)，上述方程组存在唯一解本讲主要讨论超定线性最小二乘问题的求解

3.1 问题介绍 · 83 · b 如果 A 不是满秩, 则 ATA 半正定, 此时解不唯一. 3.1.2 欠定方程组 若 m < n, 则线性方程组 Ax = b 存在无穷多个解 (假定 A 满秩). 这时我们通常寻求最小范 数解, 即所有解中范数最小的解. 于是原问题就转化为下面的约束优化问题 min Ax=b 1 2 ∥x∥ 2 2 (3.2) 对应的 Lagrange 函数为 L(x, λ) = 1 2 ∥x∥ 2 2 + λ T (Ax − b), 其中 λ = [λ1, λ2, . . . , λm] T 是 Lagrange 乘子. 此时优化问题 (3.2) 的解就是 L(x, λ) 的鞍点, 即下 面方程组的解: ∂L ∂x = x + A Tλ = 0, ∂L ∂λ = Ax − b = 0. 写成矩阵形式为 " I AT A 0 # "x λ # = " 0 b # . 如果 A 满秩, 即 rank(A) = m, 则系数矩阵非奇异 (见练习 3.1), 上述方程组存在唯一解. 本讲主要讨论超定线性最小二乘问题的求解

-84第三讲线性最小二乘问题3.2几类重要的矩阵变换矩阵计算的一个基本思想就是把复杂的问题转化为等价的且易于求解的问题.完成这个转化的基本工具就是矩阵变换，除了高等代数（或线性代数）中提到的三类初等变换以外，在矩阵计算中常用的矩阵变换还有：Gauss变换，Householder变换和Givens变换，其中Gauss变换主要用于矩阵的LU分解，Householder变换和Givens变换是正交变换，可用于计算线性最小二乘问题、矩阵特征值和奇异值问题等3.2.1初等矩阵变换本科高等代数教材中介绍了三类初等矩阵变换，这里将其推广到更一般的情形我们称E(u, V, T) = I - Tuu*(3.3)为初等矩阵变换或初等矩阵，其中u.uECn是非零向量，T是一个非零复数.因此.初等矩阵是单位矩阵的一个秩1修正高等代数中介绍的三类初等变换(elementarytransformation）分别是（以行变换为例)：（1）交换两行；(2)某行乘以一个非零常数；（3）某行乘以一个常数后加到另外一行.这三类变换所对应的矩阵为（简单示例）110Ei =, E2 =E3=1可以验证，它们都可以表示为（3.3）形式，留作练习.下面的定理给出了初等矩阵的基本性质定理3.1设E(u,U,T）是一个初等矩阵，我们有(1) det(E(u,v, T)) = 1 - Tu*u;(2)若1-Tu*u≠0,则E(u,U,T)非奇异,且(E(u,U, T)-1=E(u, U,), 其中 =TU*u-1(3)对任意非零向量a,yECn,存在u,ECn和TEC,使得E(u,u, T)r = y(板书)证明.(1)易知[3 [- ][ [ ]由行列式的乘法可知d( ) ([- ]) ([) ([β - )

· 84 · 第三讲 线性最小二乘问题 3.2 几类重要的矩阵变换 矩阵计算的一个基本思想就是把复杂的问题转化为等价的且易于求解的问题. 完成这个转化 的基本工具就是矩阵变换, 除了高等代数 (或线性代数) 中提到的三类初等变换以外, 在矩阵计算 中常用的矩阵变换还有: Gauss 变换, Householder 变换和 Givens 变换, 其中 Gauss 变换主要用于矩 阵的 LU 分解, Householder 变换和 Givens 变换是正交变换, 可用于计算线性最小二乘问题、矩阵 特征值和奇异值问题等. 3.2.1 初等矩阵变换 本科高等代数教材中介绍了三类初等矩阵变换, 这里将其推广到更一般的情形. 我们称 E(u, v, τ ) = I − τuv∗ (3.3) 为初等矩阵变换或初等矩阵, 其中 u, v ∈ C n 是非零向量, τ 是一个非零复数. 因此, 初等矩阵是单 位矩阵的一个秩 1 修正. b 高等代数中介绍的三类初等变换 (elementary transformation) 分别是 (以行变换为例): (1) 交 换两行; (2) 某行乘以一个非零常数; (3) 某行乘以一个常数后加到另外一行. 这三类变换所 对应的矩阵为 (简单示例) E1 =        1 . . . 0 1 . . . 1 0 . . . 1        , E2 =        1 . . . 1 c 1 . . . 1        , E3 =        1 . . . 1 . . . c 1 . . . 1        . 可以验证, 它们都可以表示为 (3.3) 形式, 留作练习. 下面的定理给出了初等矩阵的基本性质. 定理 3.1 设 E(u, v, τ ) 是一个初等矩阵, 我们有 (1) det(E(u, v, τ )) = 1 − τv∗u; (2) 若 1 − τv∗u ̸= 0, 则 E(u, v, τ ) 非奇异, 且 (E(u, v, τ ))−1 = E(u, v, γ), 其中 γ = τ τv∗u − 1 . (3) 对任意非零向量 x, y ∈ C n , 存在 u, v ∈ C n 和 τ ∈ C, 使得 E(u, v, τ )x = y. (板书) 证明. (1) 易知 " I 0 v ∗ 1 # "I − τuv∗ −τu 0 1 # " I 0 −v ∗ 1 # = " I −τu 0 1 − τv∗u # . 由行列式的乘法可知 det " I 0 v ∗ 1 #! · det "I − τuv∗ −τu 0 1 #! · det " I 0 −v ∗ 1 #! = det "I −τu 0 1 − τv∗u #!

3.2几类重要的矩阵变换·85.所以det (E(u, v, T) = det(I - Tuu*) = 1 - Tu*u(2)若1-Tu*u≠0,则 det(E(u,U,T)≠0,所以E(u,U,T)非奇异.通过直接计算可知E(u,V,T)E(u, V,) = I - Tu* -(1 - Tu*u)u*(1 -Tu*u)u*= I - Tuv*-TU*u-1= 1.(3)留作练习.更一般地，我们有下面的结论Sylverster降幂公式设 AE Cmxn,BE Cnxm,其中m≥n,则有det(^I - AB) = ^m-n det(^I - BA)(3.4)因此AB和BA具有相同的非零特征值（留作练习）定理3.2设AECnxn，则A非奇异当且仅当A可以分解成若干个初等矩阵的乘积3.2.2Gauss变换设lj=[0,..,0,lj+1j,.,lnj]j=1,2,n,则Gauss变换定义为1L(t) ≤ E(lj,ej, -1) = I + lje, :li+1j1..In.j1显然,Gauss变换也属于初等矩阵变换.向量l,称为Gauss向量[57].由定理3.1可知det(L(,)) = 1, (L(,)-1 = E(lj,ej,1) = E(-lj,ej -1) = L(-lj)Gauss变换是对单位矩阵的秩1修正，主要功能是消零，即将指定的某些元素变为零Gauss变换有时也称为初等下三角矩阵3.2.3Householder变换

3.2 几类重要的矩阵变换 · 85 · 所以 det ￾ E(u, v, τ )
= det(I − τuv∗ ) = 1 − τv∗u. (2) 若 1 − τv∗u ̸= 0, 则 det ￾ E(u, v, τ )
̸= 0, 所以 E(u, v, τ ) 非奇异. 通过直接计算可知 E(u, v, τ )E(u, v, γ) = I − τuv∗ − γ(1 − τv∗u)uv∗ = I − τuv∗ − τ τv∗u − 1 (1 − τv∗u)uv∗ = I. (3) 留作练习. □ 更一般地, 我们有下面的结论. Sylverster 降幂公式 设 A ∈ C m×n , B ∈ C n×m, 其中 m ≥ n, 则有 det(λI − AB) = λ m−n det(λI − BA). (3.4) 因此 AB 和 BA 具有相同的非零特征值. (留作练习) 定理 3.2 设 A ∈ C n×n , 则 A 非奇异当且仅当 A 可以分解成若干个初等矩阵的乘积. 3.2.2 Gauss 变换 设 lj = [0, . . . , 0, lj+1,j , . . . , ln,j ] T, j = 1, 2, . . . , n, 则 Gauss 变换 定义为 L(lj ) ≜ E(lj , ej , −1) = I + lje T j =              1 . . . 1 lj+1,j 1 . . . . . . ln,j 1              . 显然, Gauss 变换也属于初等矩阵变换. 向量 lj 称为 Gauss 向量 [57]. 由定理 3.1 可知 det(L(lj )) = 1, (L(lj ))−1 = E(lj , ej , 1) = E(−lj , ej , −1) = L(−lj ). b Gauss 变换是对单位矩阵的秩 1 修正, 主要功能是消零, 即将指定的某些元素变为零. b Gauss 变换有时也称为初等下三角矩阵. 3.2.3 Householder 变换

-86第三讲线性最小二乘问题定义3.1我们称矩阵H=1-βot, 0#vec", β=2(3.5)*1为Householder矩阵或Householder变换，向量称为Householder向量.我们通常将矩阵(3.5)记为 H(u).Householder矩阵也是初等矩阵有时也将Householder变换定义为H=I-2v*,wCn且lvl2=1白Householder矩阵有时也称为初等Hermite矩阵[151]从几何上看，一个Householder变换是一个关于超平面span[u]+的反射.由于Cn=span[u]④span[u]+,因此对任意向量rECn,都可写为r=au+au+,其中auEspan[u],u+Espan(u].易知=于是20*rU*rHa=T-Buu=U=U+T-ru+KlU*yU*y即Ha与&在span[u}方向有着相同的分量，而在方向的分量正好相差一个符号.也就是说，Ha是关于超平面span[u]+的镜面反射，见图3.1.因此，Householder变换也称为Householder反射或反射变换spanv)1Hx图3.1.Householder变换的几何意义下面是关于Householder矩阵的几个基本性质定理3.3设HECnxn是一个Householder矩阵，则(1)H*=H,即H是Hermite的;(2)H*H=I,即H是酉矩阵；(3)H2=I,所以H-1=H;(4) det(H) = -1;(5)H有两个互异的特征值：入=1和入=-1,其中入=1的代数重数为n-1

· 86 · 第三讲 线性最小二乘问题 定义 3.1 我们称矩阵 H = I − βvv∗ , 0 ̸= v ∈ C n , β = 2 v ∗v (3.5) 为 Householder 矩阵 或 Householder 变换, 向量 v 称为 Householder 向量. 我们通常将矩阵 (3.5) 记为 H(v). b Householder 矩阵也是初等矩阵. b 有时也将 Householder 变换定义为 H = I − 2vv∗ , v ∈ C n 且 ∥v∥2 = 1. b Householder 矩阵有时也称为初等 Hermite 矩阵 [151]. 从几何上看, 一个 Householder 变换是一个关于超平面 span{v} ⊥ 的反射. 由于 C n = span{v}⊕ span{v} ⊥, 因此对任意向量 x ∈ C n , 都可写为 x = xv + xv⊥ , 其中 xv ∈ span{v}, xv⊥ ∈ span{v} ⊥. 易知 xv = v ∗x v ∗v v. 于是 Hx = x − βvv∗x = x − 2v ∗x v ∗v v = − v ∗x v ∗v v + xv⊥ = −xv + xv⊥ , 即 Hx 与 x 在 span{v} ⊥ 方向有着相同的分量, 而在 v 方向的分量正好相差一个符号. 也就是说, Hx 是 x 关于超平面 span{v} ⊥ 的镜面反射, 见图 3.1. 因此, Householder 变换也称为Householder 反射或反射变换. 图 3.1. Householder 变换的几何意义 下面是关于 Householder 矩阵的几个基本性质. 定理 3.3 设 H ∈ C n×n 是一个 Householder 矩阵, 则 (1) H∗ = H, 即 H 是 Hermite 的; (2) H∗H = I, 即 H 是酉矩阵; (3) H2 = I, 所以 H−1 = H ; (4) det(H) = −1; (5) H 有两个互异的特征值: λ = 1 和 λ = −1, 其中 λ = 1 的代数重数为 n − 1

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) =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 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, 这就是为什么在 浮点运算中尽量使用正交变换的原因