第三讲线性最小二乘问题广义的最小二乘问题（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);.当m<n时，未知量个数大于约束个数，此时我们称问题(3.1)为欠定（或亚定）的（underdetermined)十为了讨论方便，这里假定A是满秩的3.1.1超定方程组当m>n时，线性方程组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 时, 未知量个数大于约束个数, 此时我们称问题 (3.1) 为欠定 (或亚定) 的 (underdetermined). † 为了讨论方便, 这里假定 A 是满秩的. 3.1.1 超定方程组 当 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