[人工智能数学基础] 最小二乘法

1630 字

8 分钟

[人工智能数学基础] 最小二乘法

2026-04-21

人工智能数学基础

/

数学

/

线性代数

最小二乘法的线代表示#

给定 $m$ 个数据点 $(x_i, y_i)$ ，我们希望找到一个函数 $f(x)$ 来拟合这些数据点

设拟合目标函数为:

\begin{aligned} f(x) &= a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n \\ &= \begin{bmatrix} 1 & x & x^2 & \cdots & x^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} \end{aligned}

所以我们的目标就是求解参数向量 $\mathbf{a} = [a_0, a_1, a_2, \cdots, a_n]^T$

于是可以将该拟合问题转化成一个解线性方程组问题：

将每个数据点 $(x_i, y_i)$ 代入拟合函数 $f(x)$ 中，我们得到以下超定方程组( $m > n$ )：

\begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ 1 & x_3 & x_3^2 & \cdots & x_3^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \cdots & x_m^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_m \end{bmatrix}

简写为：

\mathbf{X} \mathbf{a} = \mathbf{y}

可以看到 $X$ 的维度是 $m \times (n+1)$ ，这个方程肯定是无精确解的，所以我们需要定义残差向量(损失) $\mathbf{r}$ 为真实值与预测值之差：

\mathbf{r} = \mathbf{y} - \mathbf{X} \mathbf{a}

但 $r$ 是一个向量，为了直观且便于计算，通常使用残差平方和(RSS)来衡量拟合的好坏：

RSS(\mathbf{a}) = \Vert \mathbf{r} \Vert^2 = (\mathbf{y} - \mathbf{X} \mathbf{a})^T (\mathbf{y} - \mathbf{X} \mathbf{a})

为了找到 $RSS(\mathbf{a})$ 的最小值，我们对 $\mathbf{a}$ 求导并设置为零：

\begin{aligned} \frac{\partial RSS}{\partial \mathbf{a}} &= -2 \mathbf{X}^T (\mathbf{y} - \mathbf{X} \mathbf{a}) \\ &= -2 \mathbf{X}^T \mathbf{y} + 2 \mathbf{X}^T \mathbf{X} \mathbf{a} \end{aligned}

设置导数为零：

-2 \mathbf{X}^T \mathbf{y} + 2 \mathbf{X}^T \mathbf{X} \mathbf{a} = 0

简化得到了最小二乘法的正规方程：

\mathbf{a} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}

如果此时 $\mathbf{X}^T \mathbf{X}$ 恰好满秩(可逆)且是良态，那就直接解吧

什么是良态&病态&条件数

如果 $\mathbf{X}^T \mathbf{X}$ 不可逆或病态怎么办#

SVD奇异值分解(伪逆)#

其实用这个方法的话不到正规方程那步就可以

\begin{aligned} \mathbf{X} &= \mathbf{U} \Sigma \mathbf{V}^T \\ \mathbf{a} &= X^{-1}\mathbf{y} \\ &= X^\dagger \mathbf{y} \\ &= \mathbf{V} \Sigma^{\dagger} \mathbf{U}^T \mathbf{y} \end{aligned}

$\Sigma^{\dagger}$ 是 $\Sigma$ 的伪逆，计算方法是将 $\Sigma$ 中非零奇异值取倒数，0奇异值保持为0，然后转置

除了又重又慢以外毫无缺点

L2正则化(岭回归, Ridge, Tikhonov正则化)#

正则化就是在残差基础上加上一个惩罚项，来限制模型的复杂度，防止数值不稳定

\begin{aligned} J(\mathbf{a}) &= RSS(\mathbf{a}) + \lambda \Vert \mathbf{a} \Vert^2_2 \\ &= (\mathbf{y} - \mathbf{X} \mathbf{a})^T (\mathbf{y} - \mathbf{X} \mathbf{a}) + \lambda \mathbf{a}^T \mathbf{a} \end{aligned}

直观上理解就像用一个圆形围墙将 $\mathbf{a}$ 限制在一个范围内，当 $\Vert \mathbf{a} \Vert_2$ 变大，损失会迅速增大，所以岭回归这个名字很形象

对 $J(\mathbf{a})$ 求导并设置为零：

\begin{aligned} \frac{\partial J}{\partial \mathbf{a}} &= -2 \mathbf{X}^T (\mathbf{y} - \mathbf{X} \mathbf{a}) + 2 \lambda \mathbf{a} \\ &= -2 \mathbf{X}^T \mathbf{y} + 2 \mathbf{X}^T \mathbf{X} \mathbf{a} + 2 \lambda \mathbf{a} \end{aligned}

设置导数为零得到带L2正则化的正规方程：

\mathbf{a} = (\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X}^T \mathbf{y}

其中 $\mathbf{I}$ 是单位矩阵， $\lambda$ 是正则化强度的超参数，通常需要通过交叉验证来选择合适的 $\lambda$ 值

L1正则化(套索回归, Lasso)#

L1正则化的损失函数定义为：

\begin{aligned} J(\mathbf{a}) &= RSS(\mathbf{a}) + \lambda \Vert \mathbf{a} \Vert_1 \\ &= (\mathbf{y} - \mathbf{X} \mathbf{a})^T (\mathbf{y} - \mathbf{X} \mathbf{a}) + \lambda \sum_{i=1}^{n} |a_i| \end{aligned}

这东西有绝对值，不可导，也就没有对应的正规方程，只能用数值优化方法来求解

其中一个最常见的是坐标下降法，它的核心思想是每次固定其他参数，只优化一个参数，交替进行直到收敛

优化目标：

\min_{\mathbf{a}} \frac{1}{2} \Vert \mathbf{y} - \mathbf{X} \mathbf{a} \Vert^2_2 + \lambda \Vert \mathbf{a} \Vert_1

TIP
$RSS(\mathbf{a})$ 变成 $\dfrac{1}{2} \Vert \mathbf{y} - \mathbf{X} \mathbf{a} \Vert^2_2$ 是因为在优化过程中这两个损失是等价的(都是二次函数)， $\dfrac{1}{2}$ 纯是为了求导完跟2消掉好看

计算部分残差：

\rho_j = \mathbf{x}_j^T(\mathbf{y} - \sum_{k \neq j} a_k \mathbf{x}_k)

其中 $\mathbf{x}_j$ 是 $\mathbf{X}$ 的第 $j$ 列， $a_k$ 是参数向量 $\mathbf{a}$ 的第 $k$ 个元素，左乘 $\mathbf{x}_j^T$ 是投影(内积)，部分残差的含义为「当第 $j$ 个参数不参与拟合时，剩余的误差与第 $j$ 个特征的相关程度」

软阈值算子：

a_j = \begin{cases} \dfrac{\rho_j}{\Vert \mathbf{x}_j \Vert^2} & \text{if } j = 0 \quad \text{(偏移量)} \\[1em] \dfrac{\rho_j - \lambda}{\Vert \mathbf{x}_j \Vert^2} & \text{if } \rho_j > \lambda \\[1em] \dfrac{\rho_j + \lambda}{\Vert \mathbf{x}_j \Vert^2} & \text{if } \rho_j < -\lambda \\[1em] 0 & \text{if } |\rho_j| \leq \lambda \end{cases}

原理是 $\rho_j$ 大于 $|\lambda|$ 就往里收一点，小于就直接变0，把偏移量单独拿出来是防止偏移量被限制活动范围，导致无法拟合偏离中心的数据，除以 $\Vert \mathbf{x}_j \Vert^2$ 是为了归一化特征的尺度

部分更新残差向量：

为了避免每次迭代都重新计算巨大的矩阵乘法，所以采用增量更新

\mathbf{r}_{\text{new}} = \mathbf{r}_{\text{old}} - \mathbf{x}_j(a_j - a_j^{\text{old}})

矩阵的条件数#

理论上可逆不代表工程上真的能用，因为数值计算中可能会遇到病态矩阵，即矩阵的条件数非常大，导致解不稳定

矩阵 $\mathbf{X}$ 的条件数定义为：

\kappa(\mathbf{X}) = \Vert \mathbf{X} \Vert \cdot \Vert \mathbf{X}^{-1} \Vert

其中 $\Vert \cdot \Vert$ 表示某种范数， $\kappa(\mathbf{X}) \in [1, +\infty)$

良态: $\kappa(\mathbf{X})$ 接近1，说明 $\mathbf{X}$ 的列向量之间强线性无关(高度正交)，解稳定
病态: $\kappa(\mathbf{X})$ 很大，说明 $\mathbf{X}$ 的列向量之间接近线性相关，解不稳定，对输入的微小变化可能导致输出的巨大变化

Hager算法#

Hager算法用于估算条件数公式中的 $\Vert \mathbf{A}^{-1} \Vert$ :

优化目标：

\max \Vert \mathbf{A^{-1}} \mathbf{v} \Vert_1 \quad \text{其中 } \Vert \mathbf{v} \Vert_1 = 1

在单位球面上找一个向量 $\mathbf{v}$ ，使得被 $\mathbf{A^{-1}}$ 变换后的向量最长

先随机选一个初始向量 $\mathbf{b}$
解方程组 $\mathbf{A} \mathbf{x} = \mathbf{b}$ (使用LU分解)，得到 $\mathbf{x} = \mathbf{A^{-1}} \mathbf{b}$
对 $f= \Vert \mathbf{x} \Vert_1$ 求梯度

\nabla_{\mathbf{b}} f = (\frac{\partial \mathbf{x}}{\partial \mathbf{b}})^T \nabla_{\mathbf{\mathbf{x}}} f = (A^{-1})^T \nabla_{\mathbf{x}} f = (A^T)^{-1} \nabla_{\mathbf{x}} f

因为绝对值函数的梯度是符号函数，所以 $\nabla_{\mathbf{x}} f = \text{sign}(\mathbf{x})$

得到:

\begin{aligned} \mathbf{z} &= (A^T)^{-1} \text{sign}(\mathbf{x}) \\[1em] A^T \mathbf{z} &= \text{sign}(\mathbf{x}) \end{aligned}

解方程组 $\mathbf{A}^T \mathbf{z} = \text{sign}(\mathbf{x})$ 得到 $\mathbf{z}$ (刚刚已经对 $A$ 进行过LU分解，此步很快)，得到绝对值最大分量的索引 $j=\arg \max_i |z_i|$ ，如果 $|z_j| > z^T\mathbf{b}$ ，说明在 $j$ 的方向上可以找到一个更大的值，所以更新 $\mathbf{b}$ 为单位向量 $\mathbf{e}_j$ 并重复以上过程，否则说明已经找到最大值了，此时 $\mathbf{b}$ 就是算法结果