最小二乘问题详解2：线性最小二乘求解1. 引言复习上一篇文章《最小二乘问题详解1：线性最小二乘》中的知识，对于一个线性

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$ 的过程还是有点繁琐的，这里就不进一步推导了，免得离题太远。