最小二乘：从伪逆到正规方程

介绍最小二乘问题

伪逆与普通逆矩阵的关系

若有矩阵

$$\mathbf{A} \in \mathbb{R}^{m \times n}$$

我们可以把 $\mathbf{A}$ 看作是 $\mathbb{R}^n \to \mathbb{R}^m$ 的线性变换，它接受 $\mathbb{R}^n$ 的向量输出 $\mathbb{R}^m$ 的向量，我们知道它是 不可逆 的，即 $\mathbf{A}^{-1}$ 不存在，这是因为

• 当 $m < n$

  • 此时 $\mathbf{A}$ 是一个 列数比行数多 的 矮胖 矩阵
  • 此时可以认为 未知数比方程多 (欠定)

  在这种情况下，$\mathbf{A}$ 的 行空间 不足以撑满整个 $\mathbb{R}^n$，这是因为

  $$r = \text{rank}(\mathbf{A}) \leq m < n$$
  • 矩阵的秩不可能超过 $\min{(m,n)}$

  所以 $\mathbf{A}$ 的 零空间 维度是大于零的，因为

  $$\text{rank}(\ker \mathbf{A}) = n - r$$

  由于

  $$r < n$$

  所以

  $$\text{rank}(\ker \mathbf{A}) > 0$$

  即其中一定不止含有 零向量，于是 $\mathbf{A}$ 将 $\mathbb{R}^n$ 中的 多个向量 映射到 $\mathbb{R}^m$ 的 零向量，而我们不可能从 零向量 反向解压出映射前是哪些向量，所以我们无法得到 逆变换

  可以认为 非行满秩 的 方阵 也是这种情况

• 当 $n < m$

  • 此时 $\mathbf{A}$ 是一个 行数比列数多 的 高瘦 矩阵
  • 此时可以认为 未知数比方程少 (过定)

  由于此时

  $$\text{rank}(\mathbf{A}) \leq n < m$$

  因为 列空间 只有 $n$ 维，且小于目标空间 $m$ 维，所以它无法铺满整个目标空间，也就是说 目标空间 中一定存在一些向量是在 源空间 中找不到对应的，自然也就不可能建立一个 逆变换

也可以看看线性代数笔记中关于 伪逆 的详细讨论

但是如果我们无论如何还是想构造 $\mathbf{A}$ 的某种 “逆变换” 该怎么办？显然我们需要做出一些妥协，具体来说，对于 正向变换，也就是 $\mathbb{R}^n \to \mathbb{R}^m$

• 对于 $\mathbb{R}^n$ 内位于 行空间 中的向量，我们将它映射到 $\mathbb{R}^m$ 中的 列空间 内
• 对于 $\mathbb{R}^n$ 内位于 零空间 中的向量，我们将它映射到 零向量

对于 逆向变换，也就是 $\mathbb{R}^m \to \mathbb{R}^n$

• 对于 $\mathbb{R}^m$ 内位于 列空间 中的向量，我们将它映射回 $\mathbb{R}^n$ 中的 行空间 内
• 对于 $\mathbb{R}^m$ 内位于 左零空间 中的向量，我们将它映射回 零向量

我们称完成这种 逆变换 的矩阵为 Moore-Penrose 伪逆 $\mathbf{A}^\dagger$，可以看作 $\mathbf{A}$ 在 可逆部分 上取逆，在其余部分上取零

如果矩阵 $\mathbf{A}$ 本身是可逆的，那么它的 伪逆 就会退化为普通逆矩阵，也就是

$$\mathbf{A}^\dagger = \mathbf{A}^{-1}.$$

原因是当 $\mathbf{A}^{-1}$ 存在时，矩阵满足

$$\mathbf{A}^{-1}\mathbf{A} = \mathbf{A}\mathbf{A}^{-1} = \mathbf{I}.$$

也就是说，不存在 列空间 之外的向量，因此在这种特殊情形下，伪逆 与普通逆是一致的

找到伪逆的方法

我们已经有了 伪逆 应该做什么的指导思想，接下来就是找到一种具体的执行方案来找到它

现代最为流行且数值稳定的方法是通过 奇异值分解 SVD，它对于任意矩阵都是存在的

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

可以得到 伪逆 为

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

这里的 $\mathbf{\Sigma}^\dagger$ 是

$$\mathbf{\Sigma}^\dagger = \begin{bmatrix} \frac{1}{\sigma_1} & 0 & \cdots & 0 & \cdots & 0 \\ 0 & \frac{1}{\sigma_2} & \cdots & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & \frac{1}{\sigma_r} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 0 \end{bmatrix}_{n \times m}$$

它乘以 $\mathbf{\Sigma}$ 会得到前 $r$ 个对角线元素为 $1$，其余对角线元素全为 $0$ 的矩阵，这是最接近 单位矩阵 的情况了

可以发现之前讨论的妥协版 正向/逆向 变换，实际上就是向矩阵的 行空间 或 列空间 进行投影

$$\begin{align*} \mathbf{A}\mathbf{A}^\dagger &= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top\mathbf{V}\mathbf{\Sigma}^\dagger\mathbf{U}^\top = \mathbf{\Sigma}\mathbf{\Sigma}^\dagger = [ \ \dagger \ ]_{m\times m} \implies \text{向列空间投影} \\ \mathbf{A}^\dagger\mathbf{A} &= \mathbf{V}\mathbf{\Sigma}^\dagger\mathbf{U}^\top \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top = \mathbf{\Sigma}^\dagger\mathbf{\Sigma} = [ \ \dagger \ ]_{n\times n} \implies \text{向行空间投影} \end{align*}$$

最小二乘解

矩阵 $\mathbf{A}$ 的逆 $\mathbf{A}^{-1}$ 不存在可以理解为矩阵方程

$$\mathbf{A}\mathbf{x} = \mathbf{b}$$

无解，因为

$$\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$$

而 $\mathbf{A}^{-1}$ 不存在，这种问题很常见，但是我们可能仍然想找一个 $\mathbf{x}$ 让 $\mathbf{A}\mathbf{x}$ 尽量靠近 $\mathbf{b}$，用数学语言可以表述成，使得

$$\| \mathbf{A}\mathbf{x} - \mathbf{b} \|$$

尽可能的小，这里 $\|\cdot\|$ 表示的是 “距离” 或者 “模长” 也可以说是 范数，也就是说我们必须找到一种方式来衡量 $\mathbf{A}\mathbf{x}$ 到 $\mathbf{b}$ 的距离，最为常见的方式是 2-范数，对于向量来说是平方和的平方根

$$\|\mathbf{x}\|_2 = \sqrt{x_1^2 + x_2^2 + \cdots + x_n^2} = \sqrt{\mathbf{x}^\top\mathbf{x}}$$

所以我们要使

$$\| \mathbf{A}\mathbf{x} - \mathbf{b} \|_2$$

尽量最小，不过通常为了方便，也可以等价的去最小化它的平方

$$\|\mathbf{A}\mathbf{x} - \mathbf{b} \|_2^2$$

因为平方不会改变最小值点，我们可以继续把它展开成一个关于 $\mathbf{x}$ 的 二次型

$$\begin{align*} \|\mathbf{A}\mathbf{x} - \mathbf{b} \|_2^2 &= (\mathbf{A}\mathbf{x}-\mathbf{b})^\top(\mathbf{A}\mathbf{x}-\mathbf{b}) \\ &= (\mathbf{x}^\top\mathbf{A}^\top - \mathbf{b}^\top)(\mathbf{A}\mathbf{x}-\mathbf{b}) \quad \because (\mathbf{A}\mathbf{x}-\mathbf{b})^\top = \mathbf{x}^\top\mathbf{A}^\top - \mathbf{b}^\top \\ &= \mathbf{x}^\top\mathbf{A}^\top\mathbf{A}\mathbf{x} - \mathbf{x}^\top\mathbf{A}^\top\mathbf{b} - \mathbf{b}^\top\mathbf{A}\mathbf{x} + \mathbf{b}^\top\mathbf{b} \quad \text{展开} \\ &= \mathbf{x}^\top\mathbf{A}^\top\mathbf{A}\mathbf{x} - 2\mathbf{b}^\top\mathbf{A}\mathbf{x} + \mathbf{b}^\top\mathbf{b} \quad \because \mathbf{x}^\top\mathbf{A}^\top\mathbf{b} = \mathbf{b}^\top\mathbf{A}\mathbf{x} \ \ \text{都是标量} \\ \end{align*}$$

所以最终合并展开的式子就可以认为是我们需要进行优化的 目标函数

$$f(\mathbf{x}) = \mathbf{x}^\top\mathbf{A}^\top\mathbf{A}\mathbf{x} - 2\mathbf{b}^\top\mathbf{A}\mathbf{x} + \mathbf{b}^\top\mathbf{b}$$

按照常规方法，我们需要对该函数进行求导，并令导数为 $0$ 来寻找可能的 极值

需要注意的是，这里 $f(\mathbf{x})$ 是一个 标量值 函数，因为它输入一个向量 $\mathbf{x}$，输出一个 标量，可以写成

$$f: \mathbb{R}^n \to \mathbb{R}, \qquad \mathbf{x} \mapsto f(\mathbf{x})$$

或者，我们可以把它直接理解为 多元函数

$$f(x_1, x_2, \cdots, x_n): \mathbb{R}^n \to \mathbb{R}$$

所以对它求导就是在求 $f$ 的 梯度，遵循多元函数求导法则

第一项

对于第一项

$$\mathbf{x}^\top\mathbf{A}^\top\mathbf{A}\mathbf{x}$$

如果记 $\mathbf{M} = \mathbf{A}^\top\mathbf{A}$，由于 $\mathbf{A}^\top\mathbf{A}$ 是 对称矩阵

$$(\mathbf{A}^\top\mathbf{A})^\top = \mathbf{A}^\top\mathbf{A}$$

所以有

$$\begin{align*} \nabla_{\mathbf{x}}(\mathbf{x}^\top\mathbf{M}\mathbf{x}) &= (\mathbf{M} + \mathbf{M}^\top)\mathbf{x} \\ &= (\mathbf{M} + \mathbf{M})\mathbf{x} \\ &= 2\mathbf{M}\mathbf{x} \end{align*}$$

这只有在 $\mathbf{M}$ 是 对称矩阵 时才成立，显然

$$\mathbf{M} = \mathbf{M}^\top$$

而对于

$$\nabla_{\mathbf{x}}(\mathbf{x}^\top\mathbf{M}\mathbf{x}) = (\mathbf{M} + \mathbf{M}^\top)\mathbf{x}$$

我们可以设

$$f(\mathbf{x}) = \mathbf{x}^\top \mathbf{M} \mathbf{x}$$

若 $\mathbf{x} \in \mathbb{R}^n$，则可以写出它的分量形式

$$\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}$$

那么我们可以写出 $f(x)$ 的逐分量计算形式

$$f(x) = \sum_{i=1}^{n}\sum_{j=1}^{n} x_i\mathbf{M}_{ij}\mathbf{x}_j$$

• 注意在进行这个推理的时候，我们只要求 $\mathbf{M}$ 是 $\mathbb{R}^{n\times n}$ 的 一般方阵 并不要求它是 对称 的，所以这里的 求和上限 都是 $n$

• 若 $\mathbf{M}$ 是 非方阵，那

  $$\mathbf{x}^\top \mathbf{M} \mathbf{x}$$

  不一定是合法二次型

现在我们可以对某个分量 $x_k$ 求 偏导

$$\frac{\partial f}{\partial x_k} = \sum_{i,j}^n \frac{\partial}{\partial x_k} (x_i\mathbf{M}_{ij}x_j)$$

因为 $\mathbf{M}_{ij}$ 是常数，所以只有当 $i = k$ 或 $j = k$ 时才会有贡献，于是可以分为两部分讨论

$$\frac{\partial f}{\partial x_k} = \sum_j^n {\mathbf{M}_{kj}x_j} + \sum_i^n {x_i\mathbf{M}_{ik}}$$

• 注意对于 乘积 $\mathbf{x}^\top \mathbf{M} \mathbf{x}$ 求导，需要对左右的 $\mathbf{x}$ 各求一次，然后将它们相加，因为它们都有可能依赖 $x_k$

第一项就是 $\mathbf{M}$ 的第 $k$ 行乘以 $\mathbf{x}$，也就是

$$(\mathbf{M}\mathbf{x})_k$$

而第二项则是 $\mathbf{M}$ 的第 $k$ 列乘以 $\mathbf{x}$，即

$$(\mathbf{M}^\top\mathbf{x})_k$$

所以

$$\frac{\partial f}{\partial x_k} = (\mathbf{M}\mathbf{x})_k + (\mathbf{M}^\top\mathbf{x})_k$$

如果我们对 $\mathbf{x}$ 的所有分量都执行这个过程，而不仅仅是 $x_k$ 那么我们就可以得到

$$\nabla_{\mathbf{x}}(\mathbf{x}^\top\mathbf{M}\mathbf{x}) = (\mathbf{M} + \mathbf{M}^\top)\mathbf{x}$$

更直观一点的理解，也可以把

$$\mathbf{x}^\top \mathbf{M} \mathbf{x}$$

看成“左边一个 $\mathbf{x}$”和“右边一个 $\mathbf{x}$”各贡献一次导数：

• 左边那个 $\mathbf{x}$ 求导，给出 $\mathbf{M}\mathbf{x}$
• 右边那个 $\mathbf{x}$ 求导，给出 $\mathbf{M}^\top\mathbf{x}$

所以总共就是

$$\mathbf{M}\mathbf{x} + \mathbf{M}^\top\mathbf{x}$$

如果 $\mathbf{M}$ 对称，两边贡献一样，于是就是

$$2\mathbf{M}\mathbf{x}$$

第二项

$$-2\mathbf{b}^\top\mathbf{A}\mathbf{x}$$

我们要求它对于 $\mathbf{x}$ 的梯度

可以发现 $\mathbf{b}^\top\mathbf{A}$ 是一个 行向量 也可以写成 列向量 转置的表示方式

$$(\mathbf{A}^\top\mathbf{b})^\top$$

因为这一项对 $\mathbf{x}$ 来说其实就是一个 线性函数

先把它改写一下

$$\mathbf{b}^\top\mathbf{A}\mathbf{x}$$

注意矩阵乘法结合后

$$\mathbf{b}^\top\mathbf{A}$$

是一个行向量，所以整个式子就是

$$(\mathbf{b}^\top\mathbf{A})\mathbf{x}$$

也可以写成

$$(\mathbf{A}^\top\mathbf{b})^\top\mathbf{x}$$

因为

$$(\mathbf{b}^\top\mathbf{A})^\top=\mathbf{A}^\top\mathbf{b}$$

所以它本质上就是标准形式

$$\mathbf{c}^\top\mathbf{x} \quad \text{其中} \quad \mathbf{c}=\mathbf{A}^\top\mathbf{b}$$

而我们知道

$$\nabla_{\mathbf{x}}(\mathbf{c}^\top\mathbf{x})=\mathbf{c}$$

因此

$$\nabla_{\mathbf{x}}(\mathbf{b}^\top\mathbf{A}\mathbf{x})=\mathbf{A}^\top\mathbf{b}$$

再乘上前面的 $-2$，就得到

$$\nabla_{\mathbf{x}}(-2\mathbf{b}^\top\mathbf{A}\mathbf{x})=-2\mathbf{A}^\top\mathbf{b}$$

到这里第二项的梯度就已经出来了，接下来把它和第一项合并起来，就能得到最小二乘问题的一阶最优条件

从目标函数到法方程

我们前面把目标函数写成了

$$f(\mathbf{x})=\mathbf{x}^\top\mathbf{A}^\top\mathbf{A}\mathbf{x}-2\mathbf{b}^\top\mathbf{A}\mathbf{x}+\mathbf{b}^\top\mathbf{b}$$

因此它的梯度就是

$$\nabla_{\mathbf{x}} f(\mathbf{x}) =2\mathbf{A}^\top\mathbf{A}\mathbf{x}-2\mathbf{A}^\top\mathbf{b}$$

这里最后一项 $\mathbf{b}^\top\mathbf{b}$ 与 $\mathbf{x}$ 无关，所以导数为 $0$

令梯度为零，就得到

$$2\mathbf{A}^\top\mathbf{A}\mathbf{x}-2\mathbf{A}^\top\mathbf{b}=0$$

两边同除以 $2$

$$\mathbf{A}^\top\mathbf{A}\mathbf{x}=\mathbf{A}^\top\mathbf{b}$$

把最优解记成 $\hat{\mathbf{x}}$，就得到

$$\mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}}=\mathbf{A}^\top\mathbf{b}$$

这就是 normal equation，也就是 正规方程

直观理解

这件事本质上是在解

$$\mathbf{A}\mathbf{x}\approx \mathbf{b}$$

但很多时候这个方程组并没有精确解，比如方程数太多，或者数据本身带有噪声

于是我们退一步，不再要求 $\mathbf{A}\mathbf{x}=\mathbf{b}$ 完全成立，而是让误差向量

$$\mathbf{A}\mathbf{x}-\mathbf{b}$$

尽量小，这就是 最小二乘 的意思

从法方程的形式也能看出这一点：它并不是直接在解原方程，而是在寻找让残差最小的那个 $\mathbf{x}$