前言
本系列是科大25年秋课程《人工智能的数理科学基础》同步笔记。本文为第一章线性代数、微积分与优化的第三部分：最小二乘法相关内容。
系列文章
1、最小二乘法 (Least Squares Method)
最小二乘法是一种优化技术，常用于解决线性回归问题，旨在找到一个函数（通常是线性函数）的最佳拟合参数，使得理论预测值与实际观测值之间的残差平方和最小。
1.1 线性回归问题
考虑一个线性回归问题，我们有 $m$ 组数据 $(\boldsymbol x_i^{(1)}, \dots, \boldsymbol x_i^{(n)}, y_i)$，其中 $\boldsymbol x_i$ 是 $n$ 维特征向量， $y_i$ 是对应的观测值。我们猜测 $y_i$ 与 $\boldsymbol x_i$ 之间存在线性关系：
$$ y_i \approx w_1x_i^{(1)} + \dots + w_nx_i^{(n)} + b $$
我们的目标是找出系数 $w_1, \dots, w_n$ 和 $b$。对于第 $i$ 组数据，上述近似的残差（或误差）定义为：
$$ r_i = y_i - (w_1x_i^{(1)} + \dots + w_nx_i^{(n)} + b) $$
最小二乘法的核心思想是通过极小化残差平方和来求解这些系数：
$$ R(w_1, \dots, w_n, b) = \sum_{i=1}^m r_i^2 $$
1.2 矩阵向量形式
为了更简洁地表示和求解，我们将上述问题写成矩阵向量形式。记特征矩阵 $\boldsymbol A \in \mathbb R^{m \times (n+1)}$，参数向量 $\boldsymbol x \in \mathbb R^{n+1}$，以及观测值向量 $\boldsymbol y \in \mathbb R^m$：
$$ \boldsymbol A = \begin{bmatrix} x_1^{(1)} & \dots & x_1^{(n)} & 1 \\ \vdots & \ddots & \vdots & \vdots \\ x_m^{(1)} & \dots & x_m^{(n)} & 1 \end{bmatrix}, \quad \boldsymbol x = \begin{bmatrix} w_1 \\ \vdots \\ w_n \\ b \end{bmatrix}, \quad \boldsymbol y = \begin{bmatrix} y_1 \\ \vdots \\ y_m \end{bmatrix} $$
则残差平方和可以写为向量 2-范数的平方：
$$ R(\boldsymbol x) = \|\boldsymbol{Ax} - \boldsymbol y\|_2^2 $$
极小化 $R(\boldsymbol x)$ 的问题被称为最小二乘问题 (Least Squares Problem)。需要注意的是，方程 $R(\boldsymbol x) = 0$ 的解不一定存在，也不一定唯一。
两个线性最小二乘问题也可以合并为一个：
$$ \|\boldsymbol{A_1 x} - \boldsymbol y_1\|_2^2 + \|\boldsymbol{A_2 x} - \boldsymbol y_2\|_2^2 = \|\begin{bmatrix} \boldsymbol A_1 \ \boldsymbol A_2 \end{bmatrix} \boldsymbol x - \begin{bmatrix} \boldsymbol y_1 \ \boldsymbol y_2 \end{bmatrix}\|_2^2 = \|\boldsymbol{Ax} - \boldsymbol y\|_2^2 $$
其中
$$ \boldsymbol y = \begin{bmatrix} \boldsymbol y_1 \\ \boldsymbol y_2 \end{bmatrix},\boldsymbol A = \begin{bmatrix} \boldsymbol A_1 \\ \boldsymbol A_2 \end{bmatrix} $$
1.3 线性最小二乘问题的解法
为了找到 $R(\boldsymbol x)$ 的最小值，我们对其关于 $\boldsymbol x$ 求导，并令导数等于零。雅可比矩阵（梯度）为：
$$ \frac{\partial R}{\partial \boldsymbol x} = -2\boldsymbol A^\top (\boldsymbol{Ax} - \boldsymbol y) $$
令梯度为零，得到正规方程 (Normal Equation)：
$$ \boldsymbol A^\top \boldsymbol{Ax} = \boldsymbol A^\top \boldsymbol y $$
如果矩阵 $\boldsymbol A$ 是列满秩的，并且样本量 $m$ 足够大，那么正规方程有唯一解：
$$ \boldsymbol x = (\boldsymbol A^\top \boldsymbol A)^{-1} \boldsymbol A^\top \boldsymbol y $$
否则，如果 $\boldsymbol A^\top \boldsymbol A$ 不可逆，可以使用 Moore-Penrose 广义逆来计算，例如 Python 的 numpy.linalg.pinv 函数。
1.3.1 基于 QR 分解的解法
通过矩阵 $\boldsymbol A$ 的 QR 分解 ($\boldsymbol A = \boldsymbol{QR}$)，可以将最小二乘问题转化为：
$$ R(\boldsymbol x) = \|\boldsymbol{Ax} - \boldsymbol y\|_2^2 = \|\boldsymbol{QRx} - \boldsymbol y\|_2^2 = \|\boldsymbol{Rx} - \boldsymbol Q^\top \boldsymbol y\|_2^2 $$
如果 $\boldsymbol A$ 的秩为 $r \le n+1$，那么 $\boldsymbol R$ 和 $\boldsymbol{Rx}$ 只有前 $r$ 行非零。因此，可以进一步分解：
$$ R(\boldsymbol x) = \|(\boldsymbol{Rx})_{1:r} - (\boldsymbol Q^\top \boldsymbol y)_{1:r}\|_2^2 + \|-(\boldsymbol Q^\top \boldsymbol y)_{r+1:m}\|_2^2 $$
其中第二项是与 $\boldsymbol x$ 无关的常量。要使 $R(\boldsymbol x)$ 最小，我们只需要使第一项为零，即：
$$ (\boldsymbol{Rx})_{1:r} = (\boldsymbol Q^\top \boldsymbol y)_{1:r} $$
由于 $\boldsymbol R$ 是上三角矩阵，这个方程组可以通过回代法 (back substitution) 快速求解。
1.3.2 基于 SVD 分解的解法
通过矩阵 $\boldsymbol A$ 的奇异值分解 ($\boldsymbol A = \boldsymbol{U\Sigma V^\top}$)，最小二乘问题可以表示为：
$$ R(\boldsymbol x) = \|\boldsymbol{U\Sigma V^\top x} - \boldsymbol y\|_2^2 = \|\boldsymbol{\Sigma}(\boldsymbol V^\top \boldsymbol x) - \boldsymbol U^\top \boldsymbol y\|_2^2 $$
令 $\boldsymbol z = \boldsymbol V^\top \boldsymbol x$ 和 $\boldsymbol b = \boldsymbol U^\top \boldsymbol y$，问题变为 $\|\boldsymbol{\Sigma z} - \boldsymbol b\|_2^2$。由于 $\boldsymbol \Sigma$ 是对角矩阵（或准对角），这个系统可以很容易地求解 $\boldsymbol z$。然后通过 $\boldsymbol x = \boldsymbol V \boldsymbol z$ 得到 $\boldsymbol x$。
数值稳定性考虑：在计算 $\boldsymbol z = \boldsymbol V^\top \boldsymbol x$ 时，如果矩阵 $\boldsymbol A$ 存在很小的奇异值 $\sigma_r$，用它作分母可能会导致数值不稳定。因此，一些算法实现中会对小奇异值进行截断，以提高稳定性。
2、随机特征方法 (Random Feature Method, RFM)
随机特征方法是一种用于处理大规模非线性回归问题的技术，它通过将数据映射到高维随机特征空间，将非线性问题转化为线性问题，进而可以使用最小二乘法求解。
2.1 浅层神经网络的视角
考虑一个具有一个隐藏层的浅层神经网络的输出 $u(\boldsymbol x)$：
$$ u(\boldsymbol x) = \sum_{i=1}^m c_i \sigma(\boldsymbol w_i^\top \boldsymbol x + b_i) $$
这个表达式可以逐步拆解：
- 输入层：线性变换将输入向量 $\boldsymbol x$ 映射至 $\boldsymbol{Wx} + \boldsymbol b$。
- 隐藏层：逐项激活，其中 $\sigma(\cdot)$ 是非线性激活函数。
- 输出层：线性变换至输出 $u(\boldsymbol x)$。
在传统的深度学习方法中：
- 我们将 $u(\boldsymbol x)$ 的泛函设为目标函数，将网络参数 $c_i, \boldsymbol w_i, b_i$ 设为决策变量，建立优化问题。
- 通常使用随机梯度下降等随机优化算法求解网络参数。
2.2 随机特征方法的思想
随机特征方法的核心思想是：
- 将输入层的参数 $\boldsymbol w_i$ 和 $b_i$ 随机初始化并固定。这样，只有 $c_i$ 成为决策变量。
- 原非线性问题转化为一个关于 $c_i$ 的线性问题，可以使用确定性算法（如最小二乘法）求解。
这种方法将复杂的非线性优化问题转化为更简单的线性最小二乘问题。
2.3 深度学习方法与随机特征方法的对比
在处理某些问题时，深度学习方法和随机特征方法各有优劣：
- 深度学习方法：随机优化的每步迭代计算代价为 $O(Nm)$，其中 $N$ 是配点数，$m$ 是网络宽度。然而，当波数 $k$ 较大时，该问题难以通过随机优化方法准确求解，甚至难以收敛，这与频率原则 (frequency principle) 或谱偏置 (spectral bias) 有关。
- 随机特征方法：优化问题可以转化为关于 $c_i$ 的线性最小二乘问题，并使用矩阵分解方法求解，时间复杂度为 $O(Nm^2)$。此时，该方法对较大的波数 $k$ 并不敏感。
举例来说，以配点 $x_n$ 处的残差为例，记 $\sigma_{i,n} = \sigma(\boldsymbol w_i^\top x_n + b_i)$，则残差项为：
$$ u''(x_n) + k^2u(x_n) - g(x_n) = \sum_{i=1}^m (k^2 - 2\sigma_{i,n}(1 - \sigma_{i,n})) \sigma_{i,n}c_i - g(x_n) $$
当 $k$ 较大时，求解系数 $c_i$ 的困难程度与求解拟合问题几乎相同，随机特征方法在这种情况下表现出更好的稳定性。
3、Homework
- 设 $g(x) = \cos kx$，对于 $k = 1, 10, 100$，编写代码，应用随机特征方法求解方程 $u''(x) + k^2u(x) = g(x)$ (eq. (13))。参数 $w_i, b_i$ 的分布可以选取为 $[-1, 1]$ 上的均匀分布。
