最小二乘问题详解2：线性最小二乘求解
1. 引言
复习上一篇文章《最小二乘问题详解1：线性最小二乘》中的知识，对于一个线性问题模型：
[f(x; theta) = Atheta
]
那么线性最小二乘问题可以表达为求一组待定值(theta)，使得残差的平方和最小：
[min_{theta} |Atheta - b|^2
]
本质上是求解超定线性方程组：
[Atheta = b
]
具体的线性最小二乘解是：
[theta^* = (A^T A)^{-1} A^T b
tag{1}
]
2. 求解
2.1 问题
虽然线性最小二乘解已经给出，但是并不意味着在实际的数值计算中就能按照式(1)来进行求解。一个典型的问题就是求逆矩阵：在工程实践和数值计算中，直接求解逆矩阵通常是一个性能消耗大且可能不精确的操作，应该尽量避免。举例来说，我们按照大学本科《线性代数》课程中的方法写程序来求解一个逆矩阵，假设使用伴随矩阵法：
[A^{-1} = frac{1}{det(A)} cdot text{adj}(A)
]
其中：
- (det(A)) 是矩阵 (A) 的行列式。
- (text{adj}(A)) 是 (A) 的伴随矩阵。
为了求解伴随矩阵 (text{adj}(A))：
-
求代数余子式 (Cofactor)：对于矩阵 (A) 中的每一个元素 (a_{ij})，计算其代数余子式 (C_{ij})。
- 代数余子式 (C_{ij} = (-1)^{i+j} cdot M_{ij})
- (M_{ij}) 是删去 (A) 的第 (i) 行和第 (j) 列后得到的子矩阵的行列式（称为余子式）。
-
构造余子式矩阵：将所有代数余子式 (C_{ij}) 按照原来的位置排列，形成一个新矩阵 (C)（称为余子式矩阵）。
-
转置：将余子式矩阵 (C) 进行转置，得到的矩阵就是伴随矩阵 (text{adj}(A))。
[text{adj}(A) = C^T
] -
代入公式：将 (det(A)) 和 (text{adj}(A)) 代入公式 (A^{-1} = frac{1}{det(A)} cdot text{adj}(A)) 即可。
这里我们大概能估算，使用伴随矩阵法求逆矩阵的理论复杂度是(O(n!))，这是一个阶乘级的增长，算法效率非常低。《线性代数》中介绍的另外一种算法高斯消元法也只能达到(O(n^3))，呈指数级增加。其实效率只是一方面的问题，使用计算机求解的另外一个问题是舍入误差累积：在计算机中，浮点数运算存在固有的舍入误差；求逆过程涉及大量的除法和减法运算，这些误差会在计算过程中不断累积和传播。总而言之，使用通解求解逆矩阵，可能存在不精确且性能消耗大的问题。
2.2 QR分解
那么不使用逆矩阵怎么办呢？我们需要注意的是，最小二乘问题的本质是求解，而不是求逆矩阵，因此关键是要求解正规方程：
[A^T A theta = A^T b
]
对矩阵(A)作QR分解：
[A = Q_1 R
]
其中:
- (Q_1inmathbb R^{mtimes n}) 列正交，满足(Q_1^T Q_1 = I_n)；
- (Rinmathbb R^{ntimes n})是上三角矩阵，如果(A)列满秩，则(R)的对角元均非零，可逆。
那么把(A=Q_1R)代入正规方程，得到：
[(Q_1 R)^T (Q_1 R) x = (Q_1 R)^T b
]
左边整理，因为(Q_1^T Q_1 = I_n)：
[R^T Q_1^T Q_1 R x = R^T R x
]
右边为
[R^T Q_1^T b
]
因此正规方程等价于
[boxed{R^T R x = R^T (Q_1^T b)}
]
若(R)可逆（即(A)满秩，(operatorname{rank}(A)=n)），则(R^T)也可逆。左右两边左乘((R^T)^{-1})，得到：
[R x = Q_1^T b.
]
令(c = Q_1^top b)（这是一个长度为(n)的向量），于是我们得到一个简单的上三角线性系统：
[boxed{R x = c,qquad c = Q_1^T b}
]
这就是QR方法把正规方程化简得到的核心结果：只需解上三角方程(R x = Q_1^T b)。
以上只是对(A)列满秩的情况做了推导，如果(A)列满秩，那么QR分解可以表示为(x = R^{-1}Q_1^top b)；如果(A)列不满秩（(R)奇异），需要使用列主元QR方法对(R^T R x = R^T (Q_1^T b))进行求解，或者干脆使用下面要介绍的SVD分解（奇异值分解）法。
2.3 SVD分解
另外一种求解的方法是SVD分解。对任意矩阵(A)，存在奇异值分解：
[boxed{A = USigma V^T}
]
其中:
-
(Uinmathbb R^{mtimes m})为正交（列为左奇异向量），
-
(Vinmathbb R^{ntimes n})为正交（列为右奇异向量），
-
(Sigmainmathbb R^{mtimes n})为“对角块”矩阵，通常写成
[Sigma=begin{bmatrix}Sigma_r & 0 \ 0 & 0end{bmatrix}
]其中(Sigma_r=operatorname{diag}(sigma_1,dots,sigma_r))，((sigma_1gesigma_2gecdotsgesigma_r>0))，(r=operatorname{rank}(A))。
将SVD代入正规方程，先计算(A^top A)：
[A^T A = (USigma V^T)^T(USigma V^T) = V Sigma^T U^T U Sigma V^T = V (Sigma^T Sigma) V^T.
]
注意(U^T U=I)。而(Sigma^TSigma)是(ntimes n)的对角块矩阵，其非零对角元就是(sigma_i^2(i=1..r))，其余为零。
同样的，计算(A^T b)：
[A^T b = V Sigma^T U^T b.
]
于是正规方程变为：
[V (Sigma^T Sigma) V^T x = V Sigma^T U^T b.
]
两边左乘(V^T)，因为(V)正交，(V^TV=I)，得到：
[(Sigma^T Sigma)(V^T x) = Sigma^T (U^T b)
]
把(y=V^T x)与(c=U^T b)代入，得到更简单的对角方程：
[boxed{(Sigma^TSigma) y = Sigma^T c}
]
接下来按奇异值分块展开对角方程，先写出(Sigma)相关的形状：
[Sigma=begin{bmatrix}Sigma_r & 0 \ 0 & 0end{bmatrix},qquad
Sigma^topSigma = begin{bmatrix}Sigma_r^2 & 0 \ 0 & 0end{bmatrix}
]
对(y)和(c)也做相应分块：
[y=begin{bmatrix}y_1 y_2end{bmatrix},qquad c=begin{bmatrix}c_1 c_2end{bmatrix}
]
其中(y_1,c_1inmathbb R^r)对应非零奇异值，(y_2,c_2)对应奇异值为0的部分（维度 (n-r)）,代入得到分块方程:
[begin{bmatrix}Sigma_r^2 & 0 \ 0 & 0end{bmatrix}
begin{bmatrix}y_1 \ y_2end{bmatrix}=
begin{bmatrix}Sigma_r & 0 \ 0 & 0end{bmatrix}
begin{bmatrix}c_1 \ c_2end{bmatrix}
]
即等价于两组方程：
[Sigma_r^2 y_1 = Sigma_r c_1,qquad 0 = 0cdot c_2 (text{无约束／自由分量})
]
由于(Sigma_r)为对角且可逆，第一式等价于
[Sigma_r y_1 = c_1 quadLongrightarrowquad y_1 = Sigma_r^{-1} c_1.
]
而(y_2)（对应零奇异值的分量）在正规方程中不受约束——这反映了在列秩不足时普通最小二乘解不是唯一的（可以在零空间方向任意加解）。为得到最小范数解（惯常的选择），取 (y_2=0)。
最后回到(x)的求解，对于(y)有：
[y = begin{bmatrix} Sigma_r^{-1} c_1 \ 0 end{bmatrix}
]
将(c_1)与(c=U^top b)关系代回：
[y = begin{bmatrix} Sigma_r^{-1} & 0 \ 0 & 0end{bmatrix} U^T b
]
由于(y=V^T x)，于是：
[x = V y = V begin{bmatrix} Sigma_r^{-1} & 0 \ 0 & 0end{bmatrix} U^T b
]
定义(Sigma^+)为将非零奇异值取倒数后转置得到的伪逆矩阵（对角块为(Sigma_r^{-1})，其余为0)，则
[boxed{x^+ = V Sigma^+ U^T b}
]
这就是 最小二乘的 Moore–Penrose 伪逆解：
- 若(A)列满秩，则为唯一最小二乘解，由于那么(Sigma^+=Sigma^{-1})，SVD求解公式退化为常见的(x = VSigma^{-1}U^T b)
- 若秩亏，它给出 在所有最小二乘解中范数最小的那个（minimum-norm solution）。
2.4 比较
从以上论述可以看到，SVD分解稳定且能处理秩亏的情况，但比QR分解慢，复杂度高，通常(O(mn^2))；QR分解在列满秩、条件数不是太差时更快；若需要判定秩或求最小范数解，SVD是首选。
3. 补充
在最后补充一些基础知识，也是笔者很感兴趣的一点，那就是为什么一个矩阵A可以进行QR分解或者SVD分解呢？
3.1 QR分解
QR分解其实非常好理解，它的本质其实就是大学本科《线性代数》课程中提到的施密特（Gram–Schmidt）正交化。我们先复习一下施密特正交化相关的知识。
设有一组线性无关向量
[a_1, a_2, dots, a_n in mathbb{R}^m
]
我们想把它们变成一组正交（再归一化后变成标准正交）的向量(q_1, q_2, dots, q_n)。具体步骤如下：
-
取第一个向量，归一化：
[q_1 = frac{a_1}{|a_1|}
] -
对第 2 个向量，先减去在(q_1)上的投影：
[tilde{q}_2 = a_2 - (q_1^T a_2) q_1
]然后归一化：
[q_2 = frac{tilde{q}_2}{|tilde{q}_2|}
] -
对第 3 个向量，减去它在前两个正交向量上的投影：
[tilde{q}_3 = a_3 - (q_1^T a_3) q_1 - (q_2^T a_3) q_2
]然后归一化：
[q_3 = frac{tilde{q}_3}{|tilde{q}_3|}
] -
一般地，对第(j)个向量：
[tilde{q}_j = a_j - {sum_{i=1}^{j-1}} (q_i^T a_j) q_i,
quad q_j = frac{tilde{q}_j}{|tilde{q}_j|}
]
这样得到的({q_i})就是标准正交基，且每个(q_j)只用到了前(j-1)个。
现在把矩阵(A)看成由列向量组成：
[A = [a_1, a_2, dots, a_n] in mathbb{R}^{mtimes n}.
]
把施密特正交化写成矩阵形式，我们得到一组正交向量：
[Q_1 = [q_1, q_2, dots, q_n] in mathbb{R}^{mtimes n},
quad Q_1^T Q_1 = I_n.
]
同时，原向量可以写成：
[a_j = sum_{i=1}^j r_{ij} q_i
]
其中：
[r_{ij} = q_i^T a_j
]
把这些关系拼成矩阵形式：
[A = Q_1 R
]
其中：
- (R = (r_{ij}))是(n times n)上三角矩阵，因为第(j)列只用到前(j)个(q_i)。
- (Q_1)的列正交，所以(Q_1^T Q_1 = I)。
3.2 SVD分解
SVD分解其实也非常有意思，同样也可以顺着《线性代数》中基础知识来进行推导。首先复习一下特征值和特征向量。对于一个方阵 $ A in mathbb{R}^{n times n} $，如果存在一个非零向量 $ mathbf{v} in mathbb{R}^n $ 和一个实数 $ lambda $，使得：
[A mathbf{v} = lambda mathbf{v}
]
那么：
- $ lambda $ 称为 特征值（eigenvalue）
- $ mathbf{v} $ 称为对应于 $ lambda $ 的 特征向量（eigenvector）
接下来复习一下什么叫做对角化。如果一个(n times n)矩阵(A)可以写成：
[A = P D P^{-1}
]
其中：
- $ D $ 是一个对角矩阵（只有对角线上有元素）
- $ P $ 是一个可逆矩阵
我们就说 (A) 是 可对角化的。
而且通常：
- $ D $ 的对角元素是 $ A $ 的特征值：$ D = operatorname{diag}(lambda_1, dots, lambda_n) $
- $ P $ 的列是对应的特征向量
即：
[P = [mathbf{v}_1 mathbf{v}_2 cdots mathbf{v}_n],quad
D = begin{bmatrix}
lambda_1 & & \
& ddots & \
& & lambda_n
end{bmatrix}
]
对角化非常重要，因为对角矩阵计算非常简单，比如计算(D^k)只需把对角元各自取(k)次方即可。对角化的本质就是把复杂的线性变换，变成旋转 → 拉伸 → 逆旋转的过程。注意，不是所有矩阵都能对角化，只有当矩阵有(n)个线性无关的特征向量时，才能对角化。但是，所有对称矩阵（如 $ A^T A $）都可以对角化，而且可以使用正交矩阵对角化。
也就是说，存在正交矩阵(V)，使得：
[A^T A = V Lambda V^T,quad Lambda = operatorname{diag}(lambda_1, dots, lambda_n)
]
然后根据这个对角化公式，构造(U)和(Sigma)，最终得到SVD：
[A = U Sigma V^top
]
这里具体构造(U)和(Sigma)的过程还是有点繁琐的，这里就不进一步推导了，免得离题太远。
