第三讲线性最小二乘问题广义的最小二乘问题（LeastSquares）包括线性最小二乘问题，总体最小二乘问题，等式约束最小二乘问题，刚性加权最小二乘问题等等，它在统计学，最优化问题，材料与结构力学，信号与图像处理等方面都有着广泛的应用，是计算数学的一个重要研究分支，也是一个活跃的研究领域本讲主要介绍求解线性最小二乘问题的三种常用方法：正规方程法（也称法方程法），QR分解法和SVD分解法.一般来说，正规方程法是最快的，特别是当A的条件数较小时，正规方程法几乎与其他方法一样精确.而SVD分解法是最慢的，但结果最可靠为了方便起见，我们有时将线性最小二乘问题（3.1）简称为最小二乘问题关于最小二乘问题的更多求解方法，请参考[12,88]3.1引信考虑线性最小二乘问题minl/Ar-bll2(3.1)其中AERmxn,bERm.问题(3.1)的解称为最小二乘解·当m=n且A非奇异时，这就是一个线性方程组，解为a=A-1b；·当m>n时，约束个数大于未知量个数，此时我们称问题（(3.1）为超定的（overdetermined);.当mn时，线性方程组Ar=b的解可能不存在.此时一般考虑求解最小二乘问题（3.1).记J()=Ar-b2易知J(）是关于的二次函数，而且是凸函数(Hessen阵半正定).因此，由凸函数的性质可知，工是问题（3.1）的解当且仅当工，是J（）的稳定点（但解可能不唯一）.令其一阶导数为零，可得ATAr-ATb=O.于是将最小二乘问题转化为一个（半正定）线性方程组，这就是后面的正规方程73
第三讲 线性最小二乘问题 广义的最小二乘问题 (Least Squares ) 包括线性最小二乘问题, 总体最小二乘问题, 等式约束最小二 乘问题, 刚性加权最小二乘问题等等. 它在统计学, 最优化问题, 材料与结构力学, 信号与图像处理等方 面都有着广泛的应用, 是计算数学的一个重要研究分支, 也是一个活跃的研究领域. 本讲主要介绍求解线性最小二乘问题的三种常用方法: 正规方程法 (也称法方程法), QR 分解法和 SVD 分解法. 一般来说, 正规方程法是最快的, 特别是当 A 的条件数较小时, 正规方程法几乎与其他方 法一样精确. 而 SVD 分解法是最慢的, 但结果最可靠. 为了方便起见, 我们有时将线性最小二乘问题 (3.1) 简称为最小二乘问题. 关于最小二乘问题的更多求解方法, 请参考 [12, 88]. 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) 为超定的 (overdetermined); • 当 m n 时, 线性方程组 Ax = b 的解可能不存在. 此时一般考虑求解最小二乘问题 (3.1). 记 J(x) = ∥Ax − b∥ 2 2 . 易知 J(x) 是关于 x 的二次函数, 而且是凸函数 (Hessen 阵半正定). 因此, 由凸函数的性质可知, x∗ 是问 题 (3.1) 的解当且仅当 x∗ 是 J(x) 的稳定点 (但解可能不唯一). 令其一阶导数为零, 可得 A ⊺Ax − A ⊺ b = 0. 于是将最小二乘问题转化为一个 (半正定) 线性方程组, 这就是后面的正规方程. 73
: 74.第三讲线性最小二乘问题3.1.2欠定方程组若m<n，则线性方程组Ar=b存在无穷多个解.这时我们通常寻求最小范数解、即所有解中范数最小的解，于是原问题就转化为下面的约束优化问题minlll2(3.2)Ar=h'对应的Lagrange函数为C(,) = l + (A -b),其中入=[1，入2...*,入m]T是Lagrange乘子此时优化问题（3.2)的解就是C(r，入)的鞍点，即下面方程组的解：aoc=Ar-6=0=T+AT=0Orn写成矩阵形式为IA如果A（行）满秩，即rank(A）=m，则系数矩阵非奇异（见练习3.2），上述方程组存在唯一解本讲主要讨论超定线性最小二乘问题的求解。3.2初等变换矩阵矩阵计算的一个基本思想就是把较复杂的问题转化为等价的较简单的，易于求解问题.而完成这个转化的基本工具就是初等变换矩阵，其中常用的有三个：Gauss变换,Householder变换和Givens变换3.2.1初等矩阵我们考虑初等矩阵E(u, U, T) = I - Tuv*,其中u,ECn是非零向量，T是一个非零复数.事实上，E(u,U,T)是单位矩阵的一个秩1扰动）定理3.1设E(u,U,T）是一个初等矩阵，我们有(I) det(E(u,,T)) =1-Tu*u;(2)若1-T*u≠0,则E(u,,T)非奇异,且(E(u,U,T))-l =E(u,U,),其中=TU*u-1证明.(1)易知Tuu*0TT701-Ty'u0
· 74 · 第三讲 线性最小二乘问题 3.1.2 欠定方程组 若 m < n, 则线性方程组 Ax = b 存在无穷多个解. 这时我们通常寻求最小范数解, 即所有解中范数 最小的解. 于是原问题就转化为下面的约束优化问题 min Ax=b 1 2 ∥x∥ 2 2 (3.2) 对应的 Lagrange 函数为 L(x, λ) = 1 2 ∥x∥ 2 2 + λ ⊺ (Ax − b), 其中 λ = [λ1, λ2, . . . , λm] ⊺ 是 Lagrange 乘子. 此时优化问题 (3.2) 的解就是 L(x, λ) 的鞍点, 即下面方程 组的解: ∂L ∂x = x + A ⊺λ = 0, ∂L ∂λ = Ax − b = 0. 写成矩阵形式为 [ I A⊺ A 0 ] [x λ ] = [ 0 b ] . 如果 A (行) 满秩, 即 rank(A) = m, 则系数矩阵非奇异 (见练习 3.2), 上述方程组存在唯一解. 本讲主要讨论超定线性最小二乘问题的求解. 3.2 初等变换矩阵 矩阵计算的一个基本思想就是把较复杂的问题转化为等价的较简单的, 易于求解问题. 而完成这个 转化的基本工具就是初等变换矩阵, 其中常用的有三个: Gauss 变换, Householder 变换和 Givens 变换. 3.2.1 初等矩阵 我们考虑初等矩阵 E(u, v, τ ) = I − τuv∗ , 其中 u, v ∈ C n 是非零向量, τ 是一个非零复数. 事实上, E(u, v, τ ) 是单位矩阵的一个秩 1 扰动. 定理 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 . 证明. (1) 易知 [ I 0 v ∗ 1 ] [I − τuv∗ −τu 0 1 ] [ I 0 −v ∗ 1 ] = [ I −τu 0 1 − τv∗u ]
3.2初等变换矩阵· 75 .由行列式的乘法可知d () ([-- ]) ([) -(b )所以det (E(u, , T)) = det(I - Tuu*) = 1 - Tu*u.(2)若1-Tu*u0,则det(E(u,u,T))+0,所以E(u,u,T)非奇异.通过直接计算可知E(u,,T)E(u,v,)=I -T*-(1-Tu*u)uw*=I -TuU*(1 - T*u)uv*Tu*u-= I.口定理3.2设AECnXn，则A非奇异当且仅当A可以分解成若千个初等矩阵的乘积3.2.2Gauss变换设lj=[0,..,0,li+1,..,n.],j=1,2...,n,则Gauss变换定义为1.L(t,) ≤ E(lj,ej, -1) = I + lje] =lj+1.3:In.j1向量l,称为Gauss向量[33].由定理3.1可知det(L(,) = 1, (L(,))-1 = E(j,ej,1) = E(-lj,ej, -1) = L(-l)Gauss变换主要用于矩阵的LU分解，3.2.3Householder变换定义3.1我们称矩阵22(3.3)H=IVU*，0+Ecn110112y*v为Householder矩阵（或Householder变换，或Householder反射)，向量u称为Householder向量我们通常将矩阵(3.3)记为H(u)
3.2 初等变换矩阵 · 75 · 由行列式的乘法可知 det ([ I 0 v ∗ 1 ]) · det ([I − τuv∗ −τu 0 1 ]) · det ([ I 0 −v ∗ 1 ]) = det ([I −τu 0 1 − τv∗u ]) . 所以 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.2 设 A ∈ C n×n, 则 A 非奇异当且仅当 A 可以分解成若干个初等矩阵的乘积. 3.2.2 Gauss 变换 设 lj = [0, . . . , 0, lj+1,j , . . . , ln,j ] ⊺ , j = 1, 2, . . . , n, 则 Gauss 变换 定义为 L(lj ) ≜ E(lj , ej , −1) = I + lj e ⊺ j = 1 . . . 1 lj+1,j 1 . . . . . . ln,j 1 , 向量 lj 称为 Gauss 向量 [33]. 由定理 3.1 可知 det(L(lj )) = 1, (L(lj ))−1 = E(lj , ej , 1) = E(−lj , ej , −1) = L(−lj ). Gauss 变换主要用于矩阵的 LU 分解. 3.2.3 Householder 变换 定义 3.1 我们称矩阵 H = I − 2 v ∗v vv∗ = I − 2 ∥v∥ 2 2 vv∗ , 0 ̸= v ∈ C n , (3.3) 为 Householder 矩阵 (或 Householder 变换, 或 Householder 反射), 向量 v 称为 Householder 向量. 我们通 常将矩阵 (3.3) 记为 H(v)
.76.第三讲线性最小一乘问题从几何上看，一个Householder变换就是一个关于超平面spanfu+的反射.对任意一个向量ECn可将其写为u*ru+y≤qu+y,2*7其中 au Espan[u],y Espanu].则2Hr=r.=-20=-av+yU*U即H与在span[u]方向有着相同的分量，而在u方向的分量正好相差一个符号.也就是说,Ha是r关于超平面span[u]+的镜面反射，见图3.1.因此，Householder矩阵也称为反射矩阵.span(u)xHx图3.1.Householder变换的几何意义下面是关于Householder矩阵的几个基本性质定理3.3设HECnxn是一个Householder矩阵，则(1)H*=H,即HHermite的;(2)H*H=I,即H是酉矩阵;(3) H? = I,所以H-1 =H;(4) det(H) = -1;(5)H有两个互异的特征值：入=1和入=-1,其中入=1的代数重数为n-1.Householder矩阵的一个非常重要的应用就是可以将一个向量除第一个元素以外的所有元素都化为零我们首先给出一个引理引理3.1设工,y ECn为任意两个互异的向量,则存在一个Householder矩阵H(u)使得y=H(u)r的充要条件是2=yl2且yER.证明若llrll2=llyll2且a*yER,则y*y=a*且a*y=yr.于是[r -yl =(r -y)*(r -y) = r*r-y*r -r*y +y*y = 2(r*r -y*r).令u=T-则有H(o)=-2(-)(-)=--2(±-)(-m)yx-yl22(r*r-y*r)
· 76 · 第三讲 线性最小二乘问题 从几何上看, 一个 Householder 变换就是一个关于超平面 span{v} ⊥ 的反射. 对任意一个向量 x ∈ C n, 可将其写为 x = v ∗x v ∗v v + y ≜ αv + y, 其中 αv ∈ span{v}, y ∈ span{v} ⊥. 则 Hx = x − 2 v ∗v vv∗x = x − 2αv = −αv + y, 即 Hx 与 x 在 span{v} ⊥ 方向有着相同的分量, 而在 v 方向的分量正好相差一个符号. 也就是说, Hx 是 x 关于超平面 span{v} ⊥ 的镜面反射, 见图 3.1. 因此, 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. Householder 矩阵的一个非常重要的应用就是可以将一个向量除第一个元素以外的所有元素都化 为零. 我们首先给出一个引理. 引理 3.1 设 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
3.2初等变换矩阵· 77.即存在Householder矩阵H(u）使得y=H(）r反之，如果存在Householder矩阵H使得y=Ha，由于H是Hermite的，所以r*y=r*HrER.又口因为H是酉矩阵，所以yll2=H2=2由引理3.1，我们可以立即得到下面的定理定理3.4设=[1,2,.,n]Te R"是一个非零向量,则存在Householder矩阵H(u)使得H(u)r=ae1,其中α=2(或α=-l2),e1=[1,0,...,0]TERn设=[r1,2，...，n]TER”是一个实的非零向量，下面讨论如何计算定理3.4中Householder矩阵H(u）所对应的Householder向量u.由引理3.1的证明过程可知u=r-aei =[r1-a,a2,...,n]在实际计算中，为了尽可能地减少舍人误差，我们通常避免两个相近的数做减运算，否则就会损失有效数字，因此、我通常取=-sign(1)·l2事实上，我们也可以取α=sign(1)lrll2，但此时为了减少舍人误差，我们需要通过下面的公式来计算的第一个分量1--+路++)α=sign(1)2i=11 +α + αif sign(r1)<0(++...+)otherwiseTi+a无论怎样选取α,我们都有H=I-Buu*其中2223202-20r1*y(1-α)2+r +...+r2aU1在实数域中计算Householder向量v的算法如下，总运算量大约为3n如果是复向量，则有没有相应的结论？如果有的话，其对应的Householder向量是什么算法3.1.计算Householder向量% Given E R", compute u E R" such that Hr = ll2ei, where H = I - Buy*I:function [β,u] -House()2:n=length()(herelength()denotes thedimension of)30=3++.+4:U=T5: if g=0 then6:if <0 then
3.2 初等变换矩阵 · 77 · 即存在 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. □ 由引理 3.1, 我们可以立即得到下面的定理. 定理 3.4 设 x = [x1, x2, . . . , xn] ⊺ ∈ R n 是一个非零向量, 则存在 Householder 矩阵 H(v) 使得 H(v)x = αe1, 其中 α = ∥x∥2 (或 α = −∥x∥2), e1 = [1, 0, . . . , 0]⊺ ∈ R n. 设 x = [x1, x2, . . . , xn] ⊺ ∈ R n 是一个实的非零向量, 下面讨论如何计算定理 3.4 中 Householder 矩 阵 H(v) 所对应的 Householder 向量 v. 由引理 3.1 的证明过程可知 v = x − αe1 = [x1 − α, x2, . . . , xn] ⊺ . 在实际计算中, 为了尽可能地减少舍入误差, 我们通常避免两个相近的数做减运算, 否则就会损失有效 数字. 因此, 我通常取 α = − sign(x1) · ∥x∥2. 事实上, 我们也可以取 α = 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 + α . 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 . 在实数域中计算 Householder 向量 v 的算法如下, 总运算量大约为 3n. ♣ 如果 x 是复向量, 则有没有相应的结论? 如果有的话, 其对应的 Householder 向量是什么? 算法 3.1. 计算 Householder 向量 % Given x ∈ R n, compute v ∈ R n such that Hx = ∥x∥2e1, where H = I − βvv∗ 1: function [β, v] = House(x) 2: n = length(x) (here length(x) denotes the dimension of x) 3: σ = x 2 2 + x 2 3 + · · · + x 2 n 4: v = x 5: if σ = 0 then 6: if x1 < 0 then
: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是上三角矩阵,定义如下forij
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)=10,使得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 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 ]