1. 引言
复习上一篇文章《最小二乘问题详解1：线性最小二乘》中的知识，对于一个线性问题模型：
f(x; \theta) = A\theta
那么线性最小二乘问题可以表达为求一组待定值
\theta，使得残差的平方和最小：
\min_{\theta} \|A\theta - b\|^2
本质上是求解超定线性方程组：
A\theta = 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_1\in\mathbb R^{m\times n} 列正交，满足
Q_1^T Q_1 = I_n；
R\in\mathbb R^{n\times 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 = U\Sigma V^T}
其中:
U\in\mathbb R^{m\times m}为正交（列为左奇异向量），
V\in\mathbb R^{n\times n}为正交（列为右奇异向量），
\Sigma\in\mathbb R^{m\times n}为“对角块”矩阵，通常写成
\Sigma=\begin{bmatrix}\Sigma_r & 0 \\ 0 & 0\end{bmatrix}
其中
\Sigma_r=\operatorname{diag}(\sigma_1,\dots,\sigma_r)，
(\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0)，
r=\operatorname{rank}(A)。
将SVD代入正规方程，先计算
A^\top A：
A^T A = (U\Sigma V^T)^T(U\Sigma V^T) = V \Sigma^T U^T U \Sigma V^T = V (\Sigma^T \Sigma) V^T.
注意
U^T U=I。而
\Sigma^T\Sigma是
n\times 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^T\Sigma) y = \Sigma^T c}
接下来按奇异值分块展开对角方程，先写出
\Sigma相关的形状：
\Sigma=\begin{bmatrix}\Sigma_r & 0 \\ 0 & 0\end{bmatrix},\qquad
\Sigma^\top\Sigma = \begin{bmatrix}\Sigma_r^2 & 0 \\ 0 & 0\end{bmatrix}
对
y和
c也做相应分块：
y=\begin{bmatrix}y_1\ y_2\end{bmatrix},\qquad c=\begin{bmatrix}c_1\ c_2\end{bmatrix}
其中
y_1,c_1\in\mathbb R^r对应非零奇异值，
y_2,c_2对应奇异值为0的部分（维度
n-r）,代入得到分块方程:
\begin{bmatrix}\Sigma_r^2 & 0 \\ 0 & 0\end{bmatrix}
\begin{bmatrix}y_1 \\ y_2\end{bmatrix}=
\begin{bmatrix}\Sigma_r & 0 \\ 0 & 0\end{bmatrix}
\begin{bmatrix}c_1 \\ c_2\end{bmatrix}
即等价于两组方程：
\Sigma_r^2 y_1 = \Sigma_r c_1,\qquad 0 = 0\cdot c_2 \ (\text{无约束／自由分量})
由于
\Sigma_r为对角且可逆，第一式等价于
\Sigma_r y_1 = c_1 \quad\Longrightarrow\quad 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 & 0\end{bmatrix} U^T b
由于
y=V^T x，于是：
x = V y = V \begin{bmatrix} \Sigma_r^{-1} & 0 \\ 0 & 0\end{bmatrix} U^T b
定义
\Sigma^+为将非零奇异值取倒数后转置得到的伪逆矩阵（对角块为
\Sigma_r^{-1}，其余为0)，则
\boxed{x^+ = V \Sigma^+ U^T b}
这就是 最小二乘的 Moore–Penrose 伪逆解：
A列满秩，则为唯一最小二乘解，由于那么
\Sigma^+=\Sigma^{-1}，SVD求解公式退化为常见的
x = V\Sigma^{-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}^{m\times n}.
把施密特正交化写成矩阵形式，我们得到一组正交向量：
Q_1 = [q_1, q_2, \dots, q_n] \in \mathbb{R}^{m\times 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的过程还是有点繁琐的，这里就不进一步推导了，免得离题太远。
