最小二乘问题详解2：线性最小二乘求解

1. 引言

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

$$f(x;\theta )=A\theta$$

那么线性最小二乘问题可以表达为求一组待定值$\theta$，使得残差的平方和最小：

$$\underset{\theta }{\min }\Vert A\theta -b{\Vert }^2$$

本质上是求解超定线性方程组：

$$A\theta =b$$

具体的线性最小二乘解是：

$${\theta }^{\ast }=(A^{T}A{)}^{-1}{A}^{T}b\text{\hspace{1em}(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)$：

1. 求代数余子式 (Cofactor)：对于矩阵 $A$ 中的每一个元素 $a_{ij}$，计算其代数余子式 $C_{ij}$。

   • 代数余子式 $C_{ij}=(-1{)}^{i+j}\cdot M_{ij}$
   • $M_{ij}$ 是删去 $A$ 的第 $i$ 行和第 $j$ 列后得到的子矩阵的行列式（称为余子式）。

2. 构造余子式矩阵：将所有代数余子式 $C_{ij}$ 按照原来的位置排列，形成一个新矩阵 $C$（称为余子式矩阵）。

3. 转置：将余子式矩阵 $C$ 进行转置，得到的矩阵就是伴随矩阵 $\text{adj}(A)$。

   $$\text{adj}(A)=C^T$$

4. 代入公式：将 $\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_{1}R$代入正规方程，得到：

$$(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}Rx=R^{T}Rx$$

右边为

$$R^T{Q}_1^{T}b$$

因此正规方程等价于

$$\boxed{R^{T}Rx=R^T(Q_1^{T}b)}$$

若$R$可逆（即$A$满秩，$\mathrm{rank}(A)=n$），则$R^T$也可逆。左右两边左乘$(R^T{)}^{-1}$，得到：

$$Rx=Q_1^{T}b.$$

令$c=Q_1^{\top }b$（这是一个长度为$n$的向量），于是我们得到一个简单的上三角线性系统：

$$\boxed{Rx=c,c=Q_1^{T}b}$$

这就是QR方法把正规方程化简得到的核心结果：只需解上三角方程$Rx=Q_1^{T}b$。

以上只是对$A$列满秩的情况做了推导，如果$A$列满秩，那么QR分解可以表示为$x=R^{-1}{Q}_1^{\top }b$；如果$A$列不满秩（$R$奇异），需要使用列主元QR方法对$R^{T}Rx=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 =\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right]$$

  其中${\Sigma }_r=\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_r)$，$({\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0)$，$r=\mathrm{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=\mathrm{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^{T}V=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 =\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right],{\Sigma }^{\top }\Sigma =\left[\begin{array}{cc}{\Sigma }_r^2& 0\\ 0& 0\end{array}\right]$$

对$y$和$c$也做相应分块：

$$y=\left[\begin{array}{c}{y}_1{y}_2\end{array}\right],c=\left[\begin{array}{c}{c}_1{c}_2\end{array}\right]$$

其中$y_1,c_1\in {\mathbb{R}}^r$对应非零奇异值，$y_2,c_2$对应奇异值为0的部分（维度 $n-r$）,代入得到分块方程:

$$\left[\begin{array}{cc}{\Sigma }_r^2& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]=\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]$$

即等价于两组方程：

$${\Sigma }_r^2{y}_1={\Sigma }_r{c}_1,0=0\cdot c_2(\text{无约束／自由分量})$$

由于${\Sigma }_r$为对角且可逆，第一式等价于

$${\Sigma }_r{y}_1=c_1⟹y_1={\Sigma }_r^{-1}{c}_1.$$

而$y_2$（对应零奇异值的分量）在正规方程中不受约束——这反映了在列秩不足时普通最小二乘解不是唯一的（可以在零空间方向任意加解）。为得到最小范数解（惯常的选择），取 $y_2=0$。

最后回到$x$的求解，对于$y$有：

$$y=\left[\begin{array}{c}{\Sigma }_r^{-1}{c}_1\\ 0\end{array}\right]$$

将$c_1$与$c=U^{\top }b$关系代回：

$$y=\left[\begin{array}{cc}{\Sigma }_r^{-1}& 0\\ 0& 0\end{array}\right]U^{T}b$$

由于$y=V^{T}x$，于是：

$$x=Vy=V\left[\begin{array}{cc}{\Sigma }_r^{-1}& 0\\ 0& 0\end{array}\right]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(m{n}^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$。具体步骤如下：

1. 取第一个向量，归一化：

   $$q_1=\frac{a_1}{|a_1|}$$

2. 对第 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. 对第 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|}$$

4. 一般地，对第$j$个向量：

   $${\tilde{q}}_j=a_j-\sum _{i=1}^{j-1}(q_i^T{a}_j)q_i,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},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=PD{P}^{-1}$$

其中：

• $D$ 是一个对角矩阵（只有对角线上有元素）
• $P$ 是一个可逆矩阵

我们就说 $A$ 是 可对角化的。

而且通常：

• $D$ 的对角元素是 $A$ 的特征值：$D=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_n)$
• $P$ 的列是对应的特征向量

即：

$$P=[{\mathbf{v}}_1{\mathbf{v}}_2\cdots {\mathbf{v}}_n],D=\left[\begin{array}{ccc}{\lambda }_1& & \\ & \ddots & \\ & & {\lambda }_n\end{array}\right]$$

对角化非常重要，因为对角矩阵计算非常简单，比如计算$D^k$只需把对角元各自取$k$次方即可。对角化的本质就是把复杂的线性变换，变成旋转 → 拉伸 → 逆旋转的过程。注意，不是所有矩阵都能对角化，只有当矩阵有$n$个线性无关的特征向量时，才能对角化。但是，所有对称矩阵（如 $A^{T}A$）都可以对角化，而且可以使用正交矩阵对角化。

也就是说，存在正交矩阵$V$，使得：

$$A^{T}A=V\Lambda V^T,\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_n)$$

然后根据这个对角化公式，构造$U$和$\Sigma$，最终得到SVD：

$$A=U\Sigma V^{\top }$$

这里具体构造$U$和$\Sigma$的过程还是有点繁琐的，这里就不进一步推导了，免得离题太远。