1 引言
2 初等变换矩阵
3 QR 分解
4 奇异值分解
5 线性最小二乘问题的求解方法
6 广义逆与最小二乘
7 最小二乘扰动分析
8 最小二乘问题的推广及其

注脚

展开查看详情

1.第三讲 线性最小二乘问题 1 引言 2 初等变换矩阵 3 QR 分解 4 奇异值分解 5 线性最小二乘问题的求解方法 6 广义逆与最小二乘 7 最小二乘扰动分析 8 最小二乘问题的推广及其应用

2. 最小二乘问题 • 线性最小二乘问题 • 总体最小二乘问题 • 约束最小二乘问题 • ··· ··· 最小二乘问题在统计学, 最优化问题, 材料与结构力学, 信号与图像处理 等方面都有着广泛的应用, 是计算数学的一个重要研究分支, 也是一个活 跃的研究领域. 本讲主要介绍求解线性最小二乘问题的三种直接法. 2/90

3.1 引言 考虑线性最小二乘问题 min ∥Ax − b∥22 (3.1) x∈Rn 其中 A ∈ Rm×n , b ∈ Rm . 问题 (3.1) 的解称为最小二乘解. • 当 m = n 且 A 非奇异时, 这就是一个线性方程组, 解为 x = A−1 b; • 当 m < n 时, 未知量个数大于约束个数, 此时我们称问题 (3.1) 为欠 定 (或亚定) 的 (underdetermined). • 当 m > n 时, 约束个数大于未知量个数, 此时我们称问题 (3.1) 为超 定的 (overdetermined); † 为了讨论方便, 本讲总是假定 A 是满秩的. 3/90

4.1.1 欠定方程组 若 m < n, 则 Ax = b 存在无穷多个解. 这时我们通常寻求最小范数解, 于 是问题就转化为下面的约束优化问题 1 min ∥x∥22 (3.2) Ax=b 2 对应的 Lagrange 函数为 1 L(x, λ) = ∥x∥22 + λ (Ax − b), 2 其中 λ = [λ1 , λ2 , . . . , λm ] 是 Lagrange 乘子. 此时优化问题 (3.2) 的解就是 L(x, λ) 的鞍点, 即下面方程组的解: ∂L ∂L = x + A λ = 0, = Ax − b = 0. ∂x ∂λ 4/90

5.写成矩阵形式为 [ ][ ] [ ] I A x 0 = . A 0 λ b 如果 A (行) 满秩, 即 rank(A) = m, 则系数矩阵非奇异, 上述方程组存在唯 一解. 5/90

6.1.2 超定方程组 当 m > n 时, 线性方程组 Ax = b 的解可能不存在. 记 J(x) = ∥Ax − b∥22 . 易知 J(x) 是关于 x 的二次凸函数. 因此, x∗ 是问题 (3.1) 的解当且仅当 x∗ 是 J(x) 的稳定点. 令其一阶导数为零, 可得 A Ax − A b = 0. 于是 (3.1) 就转化为一个 (半正定) 线性方程组. 本讲我们主要讨论超定线性最小二乘问题的求解. 6/90

7.2 初等变换矩阵 2.1 初等矩阵 2.2 Gauss 变换 2.3 Householder 变换 2.4 Givens 变换 2.5 正交变换的舍入误差分析 矩阵计算的一个基本思想就是把较复杂的问题转化为等价的较简单的, 易于求解问题. 而完成这个转化的基本工具就是初等变换矩阵, 其中常用 的有三个: Gauss 变换, Householder 变换和 Givens 变换. 7/90

8.2.1 初等矩阵 我们考虑初等矩阵 E(u, v, τ ) = I − τ uv ∗ , 其中 u, v ∈ Cn 是非零向量, τ 是一个非零复数. 事实上, E(u, v, τ ) 是单位 矩阵的一个秩 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 (板书) 8/90

9.2.2 Gauss 变换 设 lj = [0, . . . , 0, lj+1,j , . . . , ln,j ] , j = 1, 2, . . . , n, 则 Gauss 变换 定义为   1  .   ..       1  △ L(lj ) = E(lj , ej , −1) = I + lj ej =   ,   lj+1,j 1   ..   .. .   .  ln,j 1 向量 lj 称为 Gauss 向量 [Golub-Van Loan 2013]. 由定理 2.1 可知 det(L(lj )) = 1, (L(lj ))−1 = E(lj , ej , 1) = E(−lj , ej , −1) = L(−lj ). Gauss 变换主要用于矩阵的 LU 分解. 9/90

10.2.3 Householder 变换 定义 我们称矩阵 2 2 H=I− ∗ vv ∗ = I − vv ∗ , 0 ̸= v ∈ Cn , (3.3) v v ∥v∥22 为 Householder 矩阵 (或 Householder 变换, 或 Householder 反射), 向量 v 称为 Householder 向量. 我们通常将矩阵 (3.3) 记为 H(v). 从几何上看, Householder 变换就是一个关于超平面 span{v}⊥ 的反射. 对任意一个向量 x ∈ Cn , 可将其写为 v∗ x △ x= v + y = αv + y, v∗ v 其中 αv ∈ span{v}, y ∈ span{v}⊥ . 则 2 Hx = x − vv ∗ x = x − 2αv = −αv + y, v∗ v 10/90

11.即 Hx 与 x 在 span{v}⊥ 方向有着相同的分量, 而在 v 方向的分量正好相 差一个符号. 也就是说, Hx 是 x 关于超平面 span{v}⊥ 的镜面反射, 见图 2.1. 因此, Householder 矩阵也称为反射矩阵. 图 2.1 Householder 变换的几何意义 11/90

12.Householder 矩阵的几个基本性质 定理 设 H ∈ Cn×n 是一个 Householder 矩阵, 则 (1) H ∗ = H, 即 H Hermite 的; (2) H ∗ H = I, 即 H 是酉矩阵; (3) H 2 = I, 所以 H −1 = H ; (4) det(H) = −1; (5) H 有两个互异的特征值: λ = 1 和 λ = −1, 其中 λ = 1 的代数重 数为 n − 1. 12/90

13.Householder 矩阵的重要应用: 将一个向量除第一个元素以外的所有元素都化为零 引理 设 x, y ∈ Cn 为任意两个互异的向量, 则存在一个 Householder 矩 阵 H(v) 使得 y = H(v)x 的充要条件是 ∥x∥2 = ∥y∥2 且 x∗ y ∈ R. (取 v = x − y 即可) (板书) 定理 设 x = [x1 , x2 , . . . , xn ] ∈ Rn 是一个非零向量, 则存在 House- holder 矩阵 H 使得 Hx = αe1 , 其中 α = ∥x∥2 (或 α = −∥x∥2 ), e1 = [1, 0, . . . , 0] ∈ Rn . 13/90

14.设 x = [x1 , x2 , . . . , xn ] ∈ Rn 是一个实的非零向量, 下面讨论如何计算定 理中的 Householder 矩阵 H(v). 由引理的证明过程可知 v = x − αe1 = [x1 − α, x2 , . . . , xn ] . 在实际计算中, 为了尽可能地减少舍入误差, 我们通常避免两个相近的数 做减法运算, 否则就会损失有效数字. 因此, 我通常取 α = − sign(x1 ) · ∥x∥2 . † 事实上, 我们也可以取 α = sign(x1 )∥x∥2 , 但此时为了减少舍入误差, 我们需要通过下面的公式来计算 v 的第一个分量 v1 α = sign(x1 )∥x∥2 , x21 − ∥x∥22 −(x22 + x23 + · · · + x2n ) v1 = x 1 − α = = . x1 + α x1 + α 14/90

15.   x1 − α, if sign(x1 ) < 0 v1 = −(x2 + x2 + · · · + x2 )   2 3 n , otherwise x1 + α 无论怎样选取 α, 我们都有 H = I − βvv ∗ 其中 2 2 2 1 β= ∗ = = =− . v v (x1 − α) + x2 + · · · + xn 2 2 2 2α − 2αx1 2 αv1 15/90

16.算法 2.1 计算 Householder 向量 % Given x ∈ Rn , compute β, v such that Hx = ∥x∥2 e1 with H = I − βvv ∗ 1: function [β, v] = house(x) 2: n = length(x) (here length(x) denotes the dimension of x) 3: σ = x22 + x23 + · · · + x2n 4: v = x 5: if σ = 0 then 6: if x1 < 0 then 7: v1 = 2x1 , β = 2/v12 8: else 9: v1 = 0, β = 0 10: end if 11: else √ 12: α = x21 + σ % α = ∥x∥2 13: if x1 < 0 then 16/90

17.14: v1 = x1 − α 15: else 16: v1 = −σ/(x1 + α) 17: end if 18: β = 2/(v12 + σ) 19: end if 上 述 算 法 总 运 算 量 大 约 为 3n, 且 具 有 很 好 的 数 值 稳 定 性 ([Wilkinson 1965]). † 在实际计算时, 我们可以将向量 v 单位化, 使得 v1 = 1. 这样, 我们就 无需为 v 另外分配空间, 而是将 v(2 : n) 存放在 x(2 : n) 中, 因为变换 后的向量 x 除第一个分量外, 其它都为零. 17/90

18.† 为了避免可能产生的溢出, 可事先将 x 单位化, 即令 x = x/∥x∥2 . 18/90

19. Householder 变换的运算量 设 A ∈ Rm×n , H = I − βvv ∗ ∈ Rm , 则 HA = (I − βvv ∗ )A = A − βvv ∗ A = A − βv(A∗ v)∗ . 因此, 在做 Householder 变换时, 并不需要生成 Householder 矩阵, 只需要 Householder 向量即可. 上面矩阵相乘的总运算量大约为 4mn. 19/90

20.2.4 Givens 变换   1  ..   .     c s    G(i, j, θ) =   .. .  ∈ Rn×n  (θ ∈ [0, 2π], i ≤ j)    −s c   ..   .  1 Givens 变换 (或 Givens 旋转, 或 Givens 矩阵), 其中 c = cos(θ), s = sin(θ). 定理 G(i, j, θ) 是正交矩阵, 且 det(G(i, j, θ)) = 1. 左乘 Givens 矩阵: 只会影响第 i 行和第 j 行的元素. 右乘 Givens 矩阵: 只会影响第 i 和第 j 列的元素. 20/90

21. ] [ c s 例 设 x = [x1 , x2 ] ∈ R , 则存在一个 Givens 变换 G = 2 ∈ R2×2 −s c 使得 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 = x21 + x22 . 通过 Givens 变换, 我们可以将向量 x ∈ R2 的第二个分量化为 0. 事实上, 对于任意一个向量 x ∈ Rn , 我们都可以通过 Givens 变换将其任 意一个位置上的分量化为 0. 更进一步, 我们也可以通过若干个 Givens 变 换, 将 x 中除第一个分量外的所有元素都化为 0. 21/90

22.算法 2.2 Givens 变换 % Given x = [a, b] , compute c, 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 22/90

23.2.5 正交变换的舍入误差分析 引理 设 P ∈ Rn×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 变换是向后稳定的. 考虑对矩阵 A 做一系列的正交变换, 则有 fl(P˜k · · · P˜1 AQ ˜1 · · · Q ˜ k ) = Pk · · · P1 (A + E)Q1 · · · Qk , 其中 ∥E∥2 = O(εu ) · (k∥A∥2 ). 这说明整个计算过程是向后稳定的. 23/90

24. ˜ 是其浮点运算近似. 当 X 作 一般地, 假设 X 是一个非奇异的线性变换, X 用到 A 上时, 我们有 = XA + E = X(A + X −1 E) = X(A + F ), △ ˜ fl(XA) 其中 ∥E∥2 = O(εu ) · ∥XA∥2 ≤ O(εu ) · ∥X∥2 · ∥A∥2 , 故 ∥F ∥2 = ∥X −1 E∥2 ≤ O(εu ) · ∥X −1 ∥2 · ∥X∥2 · ∥A∥2 = O(εu ) · κ2 (X) · ∥A∥2 , 因此, 舍入误差将被放大 κ2 (X) 倍. 当 X 是正交变换时, κ2 (X) 达到最小 值 1, 这就是为什么在浮点运算中尽量使用正交变换的原因. 24/90

25.3 QR 分解 3.1 QR 分解的存在唯一性 3.2 基于 MGS 的 QR 分解 3.3 基于 Householder 变换的 QR 分解 3.4 列主元 QR 分解 3.5 基于 Givens 变换的 QR 分解 3.6 QR 分解的稳定性 25/90

26.3.1 QR 分解的存在唯一性 定理 (QR 分解) 设 A ∈ Cm×n (m ≥ n). 则存在一个单位列正交矩阵 Q ∈ Cm×n 和一个上三角矩阵 R ∈ Cn×n , 使得 A = QR. (3.4) 证明. 设 A = [a1 , a2 , . . . , an ] ∈ Cm×n . 若 A 列满秩, 即 rank(A) = n. 则 QR 分解 (3.4) 就是对 A 的列向量组进行 Gram-Schmidt 正交化过程的矩阵描述. 具体过程见下面的算法. 26/90

27.算法 3.1 Gram-Schmidt Process 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 = qi∗ aj % qi∗ 表示共轭转置 7: qj = qj − rij qi 8: end for 9: rjj = ∥qj ∥2 10: qj = qj /rjj 11: end for 27/90

28.如果 A 不是列满秩, 我们可以做类似的正交化过程: • 如果 a1 = 0, 则令 q1 = 0; 否则令 q1 = a1 /∥a1 ∥2 ; ∑ • 对于 j = 2, 3, . . . , n, 计算 q˜j = aj − j−1 ∗ i=1 (qi aj )qi . 如果 q˜j = 0, 则令 qj = 0; 否则令 qj = q˜j /∥˜ qj ∥2 . 于是我们有 A = QR, 其中 Q = [q1 , q2 , . . . , qn ] 列正交 (但不是单位列正交, 列向量中可能有零 向量). 这里的 R = [rij ]n×n 是上三角矩阵, 定义如下   q ∗ a , for i ≤ j i j rij =  0, for i > j 易知, 如果 Q 的某一列 qk = 0, 则 R 中对应的第 k 行就全为 0. 28/90

29.设 rank(A) = l < n, 则 rank(Q) = l, 即 Q 只有 l 个非零列, 不妨设为 qi1 , qi2 , . . . , qil . 它们构成 Cm 中的一个单位正交向量组. 将其扩展成 Cm 中的一组标准正交基, 即 qi1 , qi2 , . . . , qil , q˜1 , . . . , q˜m−l . 然后我们用 q˜1 替换 Q 中的第一个零列, 用 q˜2 替换 Q 中的第二个零列, 依此类推, 将 Q 中的所有零列都替换掉. 将最后得到的矩阵记为 Q, ˜ 于是 Q∈C ˜ m×n ˜ 单位列正交. 由于 Q 中的新添加的列向量正好与 R 中的零行 相对应, 所以 ˜ = QR = A. QR 这就是 A 的 QR 分解. 29/90

30.满秩矩阵 QR 分解的存在唯一性 定理 若 A 列满秩, 并要求 R 的对角线元素都为正, 则 A 的 QR 分解存 在且唯一. (板书) † 若 A 是实矩阵, 则所有运算都是实运算, 因此 Q 和 R 都是实矩阵. † 有时也将 QR 分解定义为: 存在酉矩阵 Q ∈ Cm×m 使得 [ ] R11 A=Q , 0 其中 R11 ∈ Rn×n 为上三角矩阵. 30/90

31.若 A 不满秩, 存在置换矩阵 P , 使得 AP 的前 l 列线性无关, 其中 l = rank(A). 对 AP 进行 QR 分解, 可得: 推论 设 A ∈ Cm×n (m ≥ n). 则存在一个置换矩阵 P , 使得 [ ] R11 R12 AP = Q , 0 0 n×n 其中 Q ∈ Cm×n 单位列正交, R11 ∈ Cl×l 是非奇异上三角矩阵. † 上述结论可简化为 [ ] AP = Q R11 R12 , 其中 Q ∈ Cm×l 单位列正交, R11 ∈ Cl×l 非奇异上三角. 31/90

32.推论 (满秩分解) 设 A ∈ Cm×n 且 rank(A) = l ≤ min{m, n}, 则存在满 秩矩阵 F ∈ Cm×l 和 G ∈ Cl×n , 使得 A = F G. † 如果 A 是非奇异方阵, 则 QR 分解可用来求解线性方程组 Ax = b. † 基于 G-S 正交化的 QR 分解的运算量大约为 2mn2 . 在后面, 我们会介绍基于 Household 变换的 QR 分解, 在不需要计算 Q 的情况下, 运算量大约为 2mn2 − 2n3 /3; 如果需要计算 Q, 则需另 外大约 2mn2 − 2n3 /3 运算量. 32/90

33.3.2 基于 MGS 的 QR 分解 由于数值稳定性方面的原因, 在实际计算中, 我们一般不采用 Gram- Schmidt 过程, 取而代之的是 修正的 Gram-Schmidt 过程 (modified Gram- Schmidt process), 即 MGS . 本算法的运算量大约为 2mn2 . 算法 3.2 基于 MGS 的 QR 分解 % Given A ∈ Rm×n , compute Q = [q1 , . . . , qn ] and R such that A = QR 1: Set R = [rik ] = 0n×n (the n × n zero matrix) 2: if a1 = 0 then 3: q1 = 0 4: else 5: r11 = ∥a1 ∥2 6: q1 = a1 /∥a1 ∥2 7: end if 8: for k = 2 to n do

34. 9: qk = ak 10: for i = 1 to k − 1 do 11: rik = qi qk 12: qk = qk − rik qi 13: end for 14: if qk ̸= 0 then 15: rkk = ∥qk ∥2 16: qk = qk /rkk 17: end if 18: end for † MGS 得到的 QR 分解中, Q ∈ Rm×n , R ∈ Rn×n . 34/90

35.3.3 基于 Householder 变换的 QR 分解 由 Householder 变换的性质可知, 我们可以将任何一个非零变量 x ∈ Rn 转化成 ∥x∥2 e1 , 即除第一个元素外, 其它都为零. 假定 m = n, 即 A ∈ Rn×n , 令 H1 ∈ Rn×n 为 Householder 变换, 满足     a11 r1      21   0  a H1      ..  =  ..  .  .  . an1 0 于是   ˜12 · · · a r1 a ˜1n   0   H1 A =  . , 其中 A˜2 ∈ R(n−1)×(n−1) .   . A˜2  0 35/90

36.同样, 构造 Householder 变换 H ˜ 2 ∈ R(n−1)×(n−1) , 使得   ˜23 · · · a r2 a ˜2n   0  H˜ 2 A˜2 =  .  , 其中 A˜2 ∈ R(n−2)×(n−2)  .   . ˜ A3  0 [ ] 1 0 令 H2 = ∈ Rn×n , 则 ˜2 0 H   r1 ˜12 a a ˜13 ··· a ˜1n   0 r2 a ˜23 ˜2n  ··· a     H2 H1 A =  0 0 .  . ..   .. ˜   . A3  0 0 36/90

37.不断重复上述过程. 这样, 我们就得到一系列的矩阵 [ ] Ik−1 0 Hk = , k = 2, 3, . . . , n − 1 0 H˜k 使得   r1 a˜12 ··· a ˜1n    0 r2 ··· a ˜2n  △ Hn−1 · · · H2 H1 A =   ..  . . ..  = R. . . .  0 0 · · · rn 令 Q = (Hn−1 · · · H2 H1 )−1 = H1 H2 · · · Hn−1 , 则 A = (Hn−1 · · · H2 H1 )−1 R = QR 如果不需要生成 Q, 则运算量大约为 2mn2 − 2/3n3 . 37/90

38. 矩阵 Q 的计算 可通过下面的算法实现 Q = In , Q = QHk , k = 1, 2, . . . , n − 1. 若保留了所有 Householder 向量, 则 Q 可以通过下面的向后累积法实现: Q = In , Q = Hk Q, k = n − 1, n − 2, . . . , 1. 优点:一开始 Q 会比较稀疏, 随着迭代的进行, Q 才会慢慢变满. 运算量大约为 4(m2 n − mn2 + 13 n3 ). † 若 m > n, 则由 Householder 变换得到的 QR 分解中, Q ∈ Rm×m , R ∈ Rm×n . 38/90

39.算法 3.3 基于 Householder 变换的 QR 分解 % Given A ∈ Rm×n , compute Q ∈ Rm×m , R ∈ Rm×n such that A = QR % The upper triangular part of R is stored in the upper triangular part of A 1: Set Q = Im×m 2: for k = 1 to n do 3: x = A(k : m, k) 4: [β, vk ] = house(x) 5: vk = vk /∥vk ∥2 6: A(k : m, k : n) = (Im−k+1 − 2vk vk )A(k : m, k : n) 7: Q(1 : k − 1, k : m) = Q(1 : k − 1, k : m)(Im−k+1 − 2vk vk ) 8: Q(k : m, k : m) = Q(k : m, k : m)(Im−k+1 − 2vk vk ) 9: end for 上面的算法只是一个简单描述, 并没有考虑运算量问题. 39/90

40.† 在实际计算时, 我们通常会保留所有的 Householder 向量. 由于第 k 步中 H ˜ k 所对应的 Householder 向量 vk 的长度为 m − k + 1, 因此我们先把 vk 单位化, 使得 vk 的第一元素为 1, 这样就只要存储 vk (2 : end), 共 m − k 个元素. 这样, 我们就可以把所有的 Householder 向量存放在 A 的严格下三 角部分, 而 A 的上三角部分仍然存放 R. 在计算 Q 时采用向后累积的方法, 所以总的运算量大约为 4m2 n − 2 8 2mn2 + n3 . 若 m = n, 则运算量大约为 n3 . 3 3 † 我们也可以考虑分块 Householder QR 分解, 以便充分利用 3 级 BLAS 运算, 提高计算效率. 40/90

41.3.4 列主元 QR 分解 当 A 不是满秩时, 我们可以进行列主元 QR 分解. 定理 (列主元 QR 分解) 设 A ∈ Cm×n (m ≥ n), 且 rank(A) = k < n. 则 存在置换矩阵 P , 正交矩阵 Q ∈ Cm×m , 使得 [ ] R11 R12 AP = Q , 0 0 m×n 其中 R11 ∈ Ck×k 是非奇异上三角矩阵, 且对角线元素满足 r11 ≥ r22 ≥ · · · ≥ rkk > 0. 列主元 QR 分解的实现过程与 QR 分解基本类似, 只是每一步需要选个列 主元, 同时做一个列交换. 41/90

42.假设经过 l 步后, 我们得到下面的分解 [ ] (l) (l) ( ) (l) (l) R11 R12 △ (l) (l) AP = Q (l) = Q R 或 Q(l) AP (l) = R(l) , 0 R22 (l) 其中 P (l) 置换矩阵, Q(l) 正交矩阵, R12 ∈ R(l)×(l) 非奇异上三角. 下面进行第 l + 1 步: (l) (1) 计算 R22 所有列的范数, 若都为 0, 则算法结束, 此时必有 l = k. (l) (2) 若 l < k, 则 R22 ̸= 0. 设范数最大的列为第 il+1 列 (若有相等的, 取 其中一列即可), 范数为 rl+1,l+1 . 若 il+1 ̸= 1, 则交换 R(l) 的第 l + 1 列与第 il+1 + l 列, 并记相应的置换矩阵为 Pl+1 . (l) ˜ l+1 . (3) 以 R22 的第 1 列构造 Householder 变换 H ˜ 令 Hl+1 = Blkdiag(Il , Hl+1 ), P (l+1) (l) = P Pl+1 , 则 [ ] ( ) (l) ˜ (l) R11 R (l) (l+1) (l) 12 Hl+1 Q AP = Hl+1 R Pl+1 = R(l+1) , 0 R ˜ (l) 22 42/90

43. 其中 R ˜ (l) 的第一列除第一个元素外, 其余都是零, 且该元素等于 22 rl+1,l+1 , 即   (l) R11 ∗ ∗ [ ] (l+1) (l+1)   R11 R12 R (l+1) =  0 rl+1,l+1 ∗  (l+1) , 0 R22 0 0 ∗ (l+1) 其中 R11 ∈ R(l+1)×(l+1) 为非奇异上三角矩阵. 记 Q(l+1) Q(l) Hl+1 , 则 AP (l+1) = Q(l+1) R(l+1) . 依此类推, 直到第 k 步, 即可得 A 的列主元 QR 分解. 矩阵 R11 的对角线元素的递减关系可由列主元的选取方法推出. 43/90

44.3.5 基于 Givens 变换的 QR 分解 我们同样可以利用 Givens 变换来做 QR 分解. 设 A ∈ Rn×n , 构造 Givens 变换 G21 , 作用在 A 的最前面的两行上, 使得     a11 a ˜11 a   0   21        G21  a  a   31  =  31  .  ..   ..   .   .  an1 an1 构造 Givens 变换 G31 , 作用在 G21 A 的第 1 行和第 3 行上, 将 a31 化为零. 由于 G31 只改变第 1 行和第 3 行的值, 所以第二行的零元素维持不变. 44/90

45.以此类推, 我们可以构造一系列的 Givens 变换 G41 , G51 , . . . , Gn1 , 使得 Gn1 · · · G21 A 的第一列中除第一个元素外, 其它元素都化为零, 即   ∗ ∗ ··· ∗   0 ∗ · · · ∗  Gn1 · · · G21 A =  . . ..  .  .. .. . 0 ∗ ··· ∗ 下面我们可以对第二列进行类似的处理. 构造 Givens 变换 G32 , G42 , . . . , Gn2 , 将第二列的第 3 至第 n 个元素全化为零, 同时保持第一列不变. 其他列也做类似的处理. 最后, 通过构造 12 n(n − 1) 个 Givens 变换, 将 A 转化成一个上三角矩阵 R, 即 R = Gn,n−1 · · · G21 A 或 A = (Gn,n−1 · · · G21 ) R QR . 45/90

46.† 与 Householder 变换一样, 在进行 Givens 变换时, 我们不需要显式地 写出 Givens 矩阵. † 对于稠密矩阵而言, 基于 Givens 变换的 QR 分解的运算量比 House- holder 变换要多很多. 基于 Givens 变换的 QR 分解主要用于当矩阵的非零下三角元素相对 较少时的情形, 比如对上 Hessenberg 矩阵进行 QR 分解. † 如果 A ∈ Rm×n , m > n, 仍然可以通过 Givens 变换进行 QR 分解. 46/90

47.算法 3.4 基于 Givens 变换的 QR 分解 % Given A ∈ Rm×n , compute Q ∈ Rm×m and R ∈ Rm×n such that A = QR % The upper triangular part of R is stored in the upper triangular part of A 1: Set Q = Im×m 2: for k = 1 to n do 3: for i = k + 1 to m do 4: [c, s]=givens(akk , aik ) [ ] [ ] [ ] A(k, k : n) A(k, k : n) c s 5: =G where G = A(i, k : n) A(i, k : n) −s c 6: [Q(1 : m, k), Q(1 : m, i)] = [Q(1 : m, k), Q(1 : m, i)]G 7: end for 8: end for 47/90

48.3.6 QR 分解的稳定性 基于 Householder 变换和 Givens 变换的 QR 分解都具有很好的数值稳定 性, 参见 [Wilkinson 1965, Higham 2002]. 基于 MGS 的 QR 分解也是向后 稳定的, 参见 [Paige-Rozloˇ z ník-Strakoˇs 2006]. † 如果需要计算 Q, 则 MGS 的运算量相对较少, 因此当 A 的列向量具 有很好的线性无关性时, 我们可以使用 MGS 来计算 QR 分解. Björck [Björck 1967] 证明了, 通过 MGS 计算的矩阵 Q 满足 Q Q = I + EM GS 其中 ∥EM GS ∥2 ≈ εu κ2 (A). 而 Householder 变换计算的矩阵 Q 满足 Q Q = I + EH 其中 ∥EH ∥2 ≈ εu . 48/90

49.因此, 如果正交性至关重要, 建议使用 Householder 变换, 特别是当 A 的列 向量接近线性相关时. 49/90

50.4 奇异值分解 定理 (SVD) 设 A ∈ Cm×n (m ≥ n), 则 存 在 酉 矩 阵 U ∈ Cm×m 和 V ∈ Cn×n 使得 [ ] [ ] ∗ Σ Σ U AV = 或 A=U V ∗, (3.6) 0 0 其中 Σ = diag(σ1 , σ2 , . . . , σn ) ∈ Rn×n , 且 σ1 ≥ . . . ≥ σn ≥ 0. (板书) † 该定理也可以通过 Hermite 半正定矩阵的特征值分解来证明. 定义 分解 (3.6) 称为 A 的奇异值分解 (SVD ), σ1 , σ2 , . . . , σn 称为 A 的奇 异值, 它们是矩阵 A∗ A 的特征值的平方根. 50/90

51.† 在不做特别说明的情况下, 我们总是假定 σ1 ≥ σ2 ≥ · · · ≥ σn ≥ 0. † 如果 A ∈ Rm×n 是实矩阵, 则 U , V 也都可以是实矩阵 [Horn-Johnson 1985]. 奇异向量 矩阵 U = [u1 , u2 , . . . , um ] 和 V = [v1 , v2 , . . . , vn ] 的列向量分别称为 A 的左奇异向量和右奇异向量, 即存在关系式 Avi = σi ui , i = 1, 2, . . . , n, A∗ ui = σi vi , i = 1, 2, . . . , n, ∗ A ui = 0, i = n + 1, n + 2, . . . , m. 51/90

52. SVD 的其他形式 [ ] Σ A=U V ∗ = σ1 u1 v1∗ + σ2 u2 v2∗ + · · · + σn un vn∗ 0 记 Un = [u1 , u2 , . . . , un ] ∈ Cm×n , 则 Un 是单位列正交矩阵, 且 A = Un ΣV ∗ . (3.8) 这就是所谓的细 SVD (thin SVD ) 或 降阶 SVD (reduced SVD ), 有的文献 将 (3.8) 称为奇异值分解. 当 k < n 时, 我们称 Ak = σ1 u1 v1∗ + σ2 u2 v2∗ + · · · + σk uk vk∗ 为 A 的截断 SVD (truncated SVD ). 52/90

53.奇异值基本性质 [ ] Σ 定理 设 A = U V ∗ 是 A ∈ Cm×n (m ≥ n) 的奇异值分解, 则 0 (1) A∗ A 的特征值是 σi2 , 对应的特征向量是 vi (2) AA∗ 的特征值是 σi2 和 m − n 个零, 对应的特征向量是 ui √ (3) ∥A∥2 = σ1 , ∥A∥F = σ12 + σ22 + · · · + σn2 ; (4) 若 rank(A) = r ≤ n, 则 Ran(A) = span{u1 , u2 , . . . , ur }, Ker(A) = span{vr+1 , vr+2 , . . . , vn } (5) 设 x ∈ Cn 且 ∥x∥2 = 1, 则 σn ≤ ∥Ax∥2 ≤ σ1 ; (6) (酉不变性) 设 X ∈ Cm×m 和 Y ∈ Cn×n 是酉矩阵, 则 σi (X ∗ AY ) = σi (A). 53/90

54.定理 设 A = U ΣV ∗ 是 A ∈ Cn×n 的奇异值分解, 则: (1) | det(A)| = σ1 σ2 · · · σn ; (2) 若 A 非奇异, 则 ∥A−1 ∥2 = σn−1 , κ2 (A) = σ1 /σn ; (3) 若 A 是 Hermite 的, 且 A = U ΛU ∗ 是 A 的特征值分解, 即 U ∗ U = I, Λ = diag(λ1 , λ2 , . . . , λn ), 其中 |λ1 | ≥ |λ2 | ≥ · · · ≥ |λn |, 则 A = U ΣV 是 A 的奇异值分解, 其中 σi = |λi |, vi = sign(λi )ui , 若 λi = 0, 则取 sign(λi ) = 1 ; [ ] 0 A∗ (4) 矩阵 H = 的特征值是 ±σi , 对应的单位特征向量为 A 0 [ ] 1 vi √ . 2 ±ui 54/90

55.低秩逼近 定理 设 A = Un ΣV ∗ 是 A ∈ Cm×n 的细奇异值分解. 令 ∑ k Ak = σi ui vi∗ , i=1 则 Ak 是 min ∥A − B∥2 (3.9) B∈Cm×n , rank(B)=k 的一个解, 且 ∥A − Ak ∥2 = σk+1 . 此时, 我们称 Ak 是 A 的一个秩 k 逼近. (板书) † 对于 Frobenius 范数, 我们有类似的结论. (证明留作练习) 55/90

56.定理 (Weyl) 设 A, B ∈ Cm×n (m ≥ n), 且 rank(B) = k. 则有 max ∥Ax∥2 ≥ σk+1 (A), (3.12) x∈Ker(B), ∥x∥2 =1 和 min ∥Ax∥2 ≤ σn−k (A). (3.13) x∈Ker(B), ∥x∥2 =1 因此, σ1 (A − B) ≥ σk+1 (A), σn (A + B) ≤ σn−k (A) (3.14) 且 σi+j−1 (A) ≤ σi (B)+σj (A−B), i = 1, 2, . . . , n, j = 1, 2, . . . , n−i+1. (3.15) (板书) 56/90

57.根据 Weyl 定理, 我们可以得到矩阵奇异值的最小最大定理. 定理 设 σ1 ≥ σ2 ≥ · · · ≥ σn ≥ 0 是矩阵 A ∈ Cm×n (m ≥ n) 的奇异值, 则 σk (A) = min max ∥Ax∥2 , k = 1, 2, . . . , n dim(L)=n−k+1 x∈L, ∥x∥2 =1 和 σk (A) = max min ∥Ax∥2 , k = 1, 2, . . . , n, dim(L)=k x∈L, ∥x∥2 =1 其中 L 表示 Cn 的一个子空间. (证明留作练习) 57/90

58.引理 (交错不等式) 设 A ∈ Cm×n , Ar 是由 A 去除 r 列 (或 r 行) 后得到 的子矩阵, 则 σi (A) ≥ σi (Ar ) ≥ σi+r (A), i = 1, 2, . . . , min{m, n}. 这里, 当下标 k 大于矩阵的维数时, 我们令 σk = 0. 更进一步, 如果 B ∈ C(m−r)×(n−s) 是 A 的子矩阵, 则 σi (A) ≥ σi (B) ≥ σi+r+s (A). (板书) 58/90

59.引理 设 A ∈ Cn×n , 1 ≤ k ≤ n, 则对任意的单位列正交矩阵 Uk ∈ Cn×k 和 Vk ∈ Cn×k , 有 σi (Uk∗ AVk ) ≤ σi (A), i = 1, 2, . . . , k. (3.16) 因此, | det(Uk∗ AVk )| ≤ σ1 (A)σ2 (A) · · · σk (A). (3.17) (板书) 定理 (Weyl 不等式, 1949) 设 A ∈ Cn×n , 其奇异值和特征值分别为 σ1 ≥ σ2 ≥ · · · ≥ σn ≥ 0 和 |λ1 | ≥ |λ2 | ≥ · · · ≥ |λn |. 则 |λ1 λ2 · · · λk | ≤ σ1 σ2 · · · σk , k = 1, 2, . . . , n 且当 k = n 时, 等号成立. (板书) 59/90

60.奇异值扰动分析 定理 设 A, E ∈ Cm×n (m ≥ n). 则 |σi (A + E) − σi (A)| ≤ ∥E∥2 , i = 1, 2, . . . , n. (证明留作练习) 下面是上述定理的一个推论, 也是 SVD 的一个重要应用. 推论 设 A ∈ Cn×n , ∥ · ∥ 是任一相容矩阵范数, 则对任意的 ε > 0, 总存 在矩阵 Aε 使得 ∥A − Aε ∥ ≤ ε, 其中 Aε 具有互不相同的特征值. 因此, 可对角化矩阵在所有矩阵组成的集合中是稠密的. 60/90

61.5 线性最小二乘问题的求解方法 5.1 正规方程法 5.2 QR 分解法 5.3 奇异值分解法 61/90

62.5.1 正规方程法 定理 设 A ∈ Rm×n (m ≥ n). 则 x∗ ∈ Rn 是线性最小二乘问题 (3.1) 的 解当且仅当残量 r = b − Ax∗ 与 Ran(A) 正交, 即 x∗ 是下面的正规方 程的解 A (b − Ax) = 0 或 A Ax = A b. (3.18) (板书) 由于 A b ∈ Ran(A ) = Ran(A A) 因此正规方程 A Ax = A b 总是相容的, 即最小二乘解总是存在的,且当 A 满秩时, 这个解也是唯一的. 62/90

63.定理 设 A ∈ Rm×n (m > n). 则 A A 对称正定当且仅当 A 是列满秩 的, 即 rank(A) = n. 此时, 线性最小二乘问题 (3.1) 的解是唯一的, 其表 达式为 x = (A A)−1 A b. 当 A 列满秩时, 我们就可以使用 Cholesky 分解来求解正规方程, 总的运算 1 量大约为 mn2 + n3 + O(n2 ), 其中大部分的运算量 (mn2 ) 是用来计算 3 A A (由于 A A 对称, 因此只需计算其下三角部分). 63/90

64. 最小二乘解的几何含义 我们可以把 b 写成 b = Ax∗ + r, 其中 r⊥Ran(A). (3.19) 所以 Ax∗ 就是 b 在 Ran(A) 上的正交投影: † 最小二乘解可能并不唯一, 但分解 (3.19) 总是唯一的. 64/90

65. 最小二乘与鞍点问题 线性最小二乘问题 (3.1) 等价于 A r = 0, r = b − Ax, 即 [ ][ ] [ ] I A r b = . (3.20) A 0 x 0 这就是线性最小二乘问题 (3.1) 的增广方程 (augmented system). 事实上, 方程组 (3.20) 是下面方程组的一种特殊情形 [ ][ ] [ ] B A r f = , A 0 x g 其中 B ∈ Rm×m 对称半正定. 这就是通常所说的鞍点问题. 这个方程组 存在唯一解当且仅当 A 列满秩且矩阵 [B, A] 行满秩. 65/90

66.如果 B 对称正定, 则 r = B −1 (f − Ax), 代入第二个方程可得 A B −1 Ax = A B −1 f − g. 这就是广义的正规方程. 其所对应的广义线性最小二乘问题是 1 2 minn ∥Ax − f ∥B −1 + g x, x∈R 2 其中范数 ∥ · ∥B −1 的定义是 ∥x∥2B −1 = x B −1 x. 66/90

67.5.2 QR 分解法 假定 A ∈ Rm×n (m ≥ n) 满秩, 用三种不同方法来推导. ˆ ∈ Rm×m . 于是有 (1) 将 Q 的扩充成一个正交矩阵, 记为 [Q, Q] 2 ∥Ax − b∥22 = [Q, Q] ˆ (Ax − b) 2 2 ˆ (QRx − b) = [Q, Q] 2 [ ] 2 Rx − Q b = −Qˆ b 2 2 2 2 = ∥Rx − Q b∥2 ˆ b + Q ≥ Q ˆ b , 2 2 等号成立当且仅当 Rx = Q b. 所以最小二乘解为 x∗ = R−1 Q b. 67/90

68.(2) 将 b 写成 b = (QQ + I − QQ )b, 则 Ax − b = Ax − (QQ + I − QQ )b = (Ax − QQ b) − (I − QQ )b. 由于 QQ 是 Ran(A) 上的正交投影, 因此 Ax − QQ b 与 (I − QQ )b 正交. 所以 ∥Ax − b∥22 = ∥Ax − QQ b∥22 + ∥(I − QQ )b∥22 = ∥Rx − Q b∥22 + ∥(I − QQ )b∥22 ≥ ∥(I − QQ )b∥22 , 等号成立当且仅当 Rx = Q b. 所以最小二乘解为 x∗ = R−1 Q b. 68/90

69.(3) 解正规方程. 由定理 5.1 可知, 最小二乘解为 x∗ = (A A)−1 A b = (R Q QR)−1 R Q b = (R R)−1 R Q b = R−1 Q b. † 用 QR 分解来求最小二乘解的运算量大约为 2mn2 (如果采用 House- holder 变换的话, 运算量大约为 2mn2 − 32 n3 ). 当 m ≫ n 时, 大约为 正规方程的两倍. 当 m = n 时, 几乎相同. † 通常 QR 算法比较稳定, 是求解最小二乘问题的首选方法, 特别是当 A 条件数较大 (病态) 时. 69/90

70.5.3 奇异值分解法 [ ] Σ 设A∈R m×n 列满秩, A = U ˜ ], 则 V 是 A 的奇异值分解, U = [Un , U 0 [ ] 2 [ ] 2 Σ Σ ∥Ax − b∥22 = U V x−b = V x − [Un , U ˜] b 0 0 2 2 [ ] 2 ΣV x − Un b = −U˜ b 2 = ∥ΣV x − Un b∥22 + ∥U ˜ b∥22 ≥ ∥U ˜ b∥22 , 等号当且仅当 ΣV x − Un b = 0 时成立, 即 x = (ΣV )−1 Un b = V Σ−1 Un b. 这就是线性最小二乘问题 (3.1) 的解. 70/90

71.6 广义逆与最小二乘 6.1 广义逆的定义 6.2 广义逆的基本性质 6.3 广义逆的计算 6.4 广义逆与线性最小二乘 71/90

72.6.1 广义逆的定义 广义逆的概念最早由 Moore 于 1920 年提出: 设 A ∈ Cm×n , 若 X ∈ Cn×m 满足 AX = PRan(A) , XA = PRan(X) , (3.21) 即 AX 和 XA 分别为 Ran(A) 和 Ran(X) 上的正交投影, 则称 X 是 A 的 广义逆. 72/90

73.1955 年, Penrose 利用四个矩阵方程给出了广义逆的定义. 定义 设 A ∈ Cm×n , 若 X ∈ Cn×m 满足 AXA = A (3.22) XAX = X (3.23) (AX)∗ = AX (3.24) ∗ (XA) = XA. (3.25) 则称 X 为 A 的广义逆 (或 Moore-Penrose 逆, 简称 MP 逆), 记为 A† . † 可以证明, 以上两个定义是等价的. † 若 A ∈ Cn×n 非奇异, 则 A† = A−1 . 73/90

74. 定理 设 A ∈ Cm×n , 则满足矩阵方程 (3.22)-(3.25) 的矩阵 X ∈ Cn×m 存 在且唯一, 即广义逆存在且惟一. 广义逆与 SVD 设 A ∈ Cm×n , rank(A) = r > 0, A 的 SVD 为 [ ] Σ1 0 A=U V ∗ , Σ1 = diag(σ1 , σ2 , . . . , σr ) ∈ Rr×r . 0 0 则容易验证 [ ] −1 Σ 0 A† = V 1 U ∗. 0 0 74/90

75.6.2 广义逆的基本性质 定理 设 A ∈ Cm×n , 则 (1) (A† )† = A; (2) (A )† = (A† ) , (A∗ )† = (A† )∗ ; (3) rank(A) = rank(A† ) = rank(A† A); (4) (AA∗ )† = (A∗ )† A† , (A∗ A)† = A† (A∗ )† ; (5) (AA∗ )† AA∗ = AA† , (A∗ A)† A∗ A = A† A; (6) A† = (A∗ A)† A∗ = A∗ (AA∗ )† , 特别地, 若 A 列满秩, 则 A† = (A∗ A)−1 A∗ , 若 A 行满秩, 则 A† = A∗ (AA∗ )−1 ; (7) 若 U, V 是酉矩阵, 则 (U AV )† = V ∗ A† U ∗ . 75/90

76.† 一般来说, 当 A, B 是方阵时, • (AB)† ̸= B † A† ; • AA† ̸= A† A; • (Ak )† ̸= (A† )k ; • A 和 A† 的非零特征值并不是互为倒数. 设 A ∈ Cm×n , 则 Ran(AA† ) = Ran(AA∗ ) = Ran(A), Ran(A† A) = Ran(A∗ A) = Ran(A∗ ) = Ran(A† ), Ker(AA† ) = Ker(AA∗ ) = Ker(A∗ ) = Ker(A† ), Ker(A† A) = Ker(A∗ A) = Ker(A). 76/90

77.6.3 广义逆的计算 • 利用 SVD, 但运算量较大. • 利用满秩分解 定理 设 A ∈ Rm×n . (1) 若 A 是列满秩矩阵, 则 A† = (A∗ A)−1 A∗ ; (2) 若 A 是行满秩矩阵, 则 A† = A∗ (AA∗ )−1 ; (3) 若 A 的秩是 r ≤ min{m, n}, 且其满秩分解为 A = F G, 其中 F ∈ Rm×r , G ∈ Rr×n , 则 A† = G† F † = G∗ (GG∗ )−1 (F ∗ F )−1 F ∗ . 77/90

78.• 利用 QR 分解 定理 设 A ∈ Rm×n (m ≥ n) 是列满秩矩阵, 其 QR 分解为 A = QR, 其中 Q ∈ Rm×n , R ∈ Rn×n , 则 A† = R−1 Q∗ . • 其他算法 其他比较重要的算法有: Greville 递推算法, Cline 算法等. 78/90

79.6.4 广义逆与线性最小二乘 定理 设 A ∈ Rm×n , 则线性最小二乘问题 (3.1) 的解为 x = A† b + (I − PA )z, ∀ z ∈ Rn . (3.26) † 通常, 线性最小二乘问题的解 (3.26) 不是唯一的. 但当 A 列满秩时, PA = I, 此时解唯一. 定理 设 A ∈ Rm×n 的解集为 S, 则 min ∥x∥2 (3.27) x∈S 存在唯一解, 即满足 (3.27) 的线性最小二乘问题 (3.1) 的解存在且唯一. 79/90

80.7 最小二乘扰动分析 定理 设 A ∈ Rm×n (m ≥ n) 且 rank(A) = n. 设 x 是线性最小二乘问题 ˜ 极小化 ∥(A + δA)˜ (3.1) 的解, x x − (b + δb)∥2 , 则 { } ∥˜x − x∥2 2κ2 (A) ≤ε· + κ22 (A) tan θ + O(ε2 ), ∥x∥2 cos θ 其中 κ2 (A) = σ1 (A)/σn (A), θ 为 b 与 Ran(A) 的夹角, { } △ ∥δA∥2 ∥δb∥2 ε = max , , ∥A∥2 ∥b∥2 并假定 ε · κ2 (A) < 1 (确保 A + δA 满秩, 从而 x ˜ 唯一确定). 我们记 △ 2κ2 (A) κLS = + κ22 (A) tan θ, cos θ 80/90

81.这就是最小二乘问题的条件数. 当 θ = 0 时, b ∈ Ran(A), 此时 κLS = 2κ2 (A); 当 θ = π/2 时, b⊥Ran(A), 此时最小二乘解为 x = 0, 而 κLS = ∞; 当 0 < θ < π/2 时, κLS = O(κ22 (A)). 定义残量 r = b−Ax, r˜ = (b+δb)−(A+δA)˜ x, 我们有下面的性质 [Higham 2002] ∥˜ r − r∥2 ≤ 1 + 2ε · κ2 (A). ∥r∥2 当我们使用 QR 分解或 SVD 分解求解最小二乘问题时, 由于采用的是正 交变换, 它们都是数值稳定的. 而正规方程涉及求解方程组 A Ax = A b, 其精度依赖于条件数 κ2 (A A) = κ22 (A), 因为其误差是以 κ22 (A) 倍数增 长. 因此当 A 的条件数较大时, 正规方程法的精度会大大降低. 81/90

82.8 最小二乘问题的推广及其应用 8.1 正则化 8.2 加权正则化 8.3 约束最小二乘问题 8.4 应用: 多项式数据拟合 8.5 应用: 线性预测 8.6 其它应用 82/90

83.8.1 正则化 在求解超定线性方程组时, 我们极小化 ∥Ax − b∥22 , 而对于欠定线性方程 组, 我们极小化的是 ∥x∥22 . 两者合起来就是 1 α min ∥Ax − b∥22 + ∥x∥22 , (3.28) x∈Rn×n 2 2 其中 α > 0 是权系数. 对应的目标函数记为 J(x) = ∥Ax − b∥22 + α∥x∥22 . 当 α > 0 时, J(x) 是一个严格凸的二次函数, 因此存在唯一的最小值点. 关于 x 求导后, 令其等于零, 可得 (A A + αI)x = A b. 由于 A A + αI 对称正定, 故非奇异. 所以问题 (3.28) 的唯一解为 x = (A A + αI)−1 AT b. 83/90

84.8.2 加权正则化 一类应用更广泛的问题是下面的加权正则化问题 1 α min ∥Ax − b∥22 + ∥W x∥22 , (3.29) x∈Rn×n 2 2 其中 W ∈ Rp×n 是广义加权矩阵, 可以是非负对角矩阵, 对称正定矩阵, 也可以是一般矩阵 (如一阶差分算子, 二阶差分算子等). 注意, W 不一定 要求是方阵. 经过类似的推导, 问题 (3.29) 的解满足 (A A + αW W )x = A b. 如果 A A + αW W 非奇异, 则存在唯一解 x = (A A + αW W )−1 AT b. 84/90

85.8.3 约束最小二乘问题 考虑带有约束的最小二乘问题 1 min ∥Ax − b∥22 (3.30) x∈Rn×n 2 s.t. Bx = f 其中 Bx = f 是约束条件. 对应的 Lagrange 函数为 1 J(x) = ∥Ax − b∥22 + λ (Bx − f ). 2 分别对 x 和 λ 求一阶导数, 并令其等于零, 可得 [ ][ ] [ ] A AB x A b = . B 0 λ f 85/90

86.如果 A A 非奇异, 且 B 行满秩, 则存在唯一解 [ ]−1 [ ] λ = B(A A)−1 B B(A A)−1 A b − f x = (A A)−1 (A b − B λ). 86/90

87.8.4 应用: 多项式数据拟合 已知 n 个点 {(ti , fi )}ni=1 , 寻找一个低次多项式来拟合这些数据. 设拟合多项式为 p(x) = a0 + a1 t + a2 t2 + · · · + am tm . 通常 m ≪ n. 将上述 n 个点代入可得       1 t1 t21 · · · tm 1 a0 f1       1 t2 t22 · · · tm   a1   f2  . . . ..  2  . = .  或 Ax = f, . . .   .   .  . . . .   .   .  1 tn tn · · · tn n×(m+1) am 2 m fn 其中 A ∈ Rn×(m+1) , x = [a0 , a1 , . . . , am ] . 由于 m ≪ n, 该方程组是超定 的, 解通常是不存在的. 因此, 我们寻找一个近似解, 使得残量 ∥f − Ax∥2 最小, 即求解最小二乘问题 min ∥f − Ax∥22 x∈Rm+1 87/90

88.8.5 应用: 线性预测 预测一个时间序列的未来走向, 一个常用方法就是线性预测. 假定在 tk 时 刻的值 fk 线性依赖于其前 m 个时刻的值 fk−1 , fk−2 , . . . , fk−m , 即 fk = a1 fk−1 + a2 fk−2 + · · · + am fk−m . (3.31) 现 在 已 经 测 得 该 时 间 序 列 的 前 n 个 值 fi , i = 0, 1, 2, . . . , n − 1, 其 中 n ≫ m. 如何预测其未来的取值? 将现有的数据代入关系式 (3.31) 可得       fm−1 fm−2 fm−3 · · · f0 a1 fm        fm fm−1 fm−2 · · · f1   a2  fm+1   .       . .. .. ..   ..   ..  或 Ax = f, =  . . . .   .   .  fn−2 fn−3 fn−4 · · · fn−m−1 (n−m)×m am fn−1 这也是一个超定问题, 因此需要求解最小二乘 minx∈Rm ∥f − Ax∥22 88/90

89.8.6 其它应用 • 信号光滑 : 在获取数字信号时, 由于各种各样的原因, 最后得到的信 号总会带有一定的噪声. 去噪是数字信号和图像处理中的一个基本 问题. 其中一个有效方法就是加权最小二乘法. • 反卷积: 一个典型例子是图像处理去模糊 • 信号补全 89/90

90.课后习题 见讲义