线性回归中的最小二乘法：直接法与梯度下降的比较

最小二乘法（Least Squares Method）学习博客

0. 直观理解最小二乘法

最小二乘法（Least Squares Method）是一种用于数据拟合的方法，它的核心思想是“最小化误差的平方和”。

假设我们有一组数据点，并希望找到一条最优的直线或曲线来尽可能贴合这些点。由于数据通常存在噪声或者误差，无法完美拟合所有点，因此最小二乘法的目标是找到一个最佳拟合，使得所有点到拟合曲线的垂直距离的平方和最小。

通俗地说，就像是在一堆散落的数据点中找一根“最合理”的线，使得数据点到这条线的总体偏差最小。

1. 最小二乘法解决的问题

在数据拟合或机器学习任务中，我们经常需要找到一个最优的模型，使其能够最好地描述数据的规律。例如，给定一组数据点 $x_i, y_i)$ ，希望找到一个线性模型：

$y = w x + b$

或者更一般的线性回归模型：

$y = w_1 x_1 + w_2 x_2 + ... + w_n x_n + b$

最小二乘法的目标是找到参数 $\mathbf{w}$ 和 $b$ ，使得模型预测值与真实值之间的误差平方和最小，即：

$\min_{\mathbf{w}, b} \sum_{i=1}^{m} (y_i - (w_1 x_{i1} + w_2 x_{i2} + ... + w_n x_{in} + b))^2$

这里， $m$ 是样本数， $n$ 是特征数。

2. 最小二乘法的求解方法

2.1 直接法（正规方程）

假设我们用矩阵形式表示线性回归模型：

$\mathbf{y} = X \mathbf{w} + \varepsilon$

其中，

$\mathbf{y} \in \mathbb{R}^{m \times 1}$ 是目标值向量。
$\in \mathbb{R}^{m \times (n+1)}$ 是输入数据矩阵（包括一个全 1 的列表示偏置项）。
$\mathbf{w} \in \mathbb{R}^{(n+1) \times 1}$ 是待求参数向量。
$\varepsilon$ 是误差项。

最小二乘法的目标是最小化误差平方和，即最小化：

$J(\mathbf{w}) = ||X\mathbf{w} - \mathbf{y}||^2$

通过求导并设导数为 0，可以推导出最优解的封闭形式：

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

2.2 迭代法（梯度下降）

当数据规模很大时，直接求解 $\mathbf{w} = (X^T X)^{-1} X^T \mathbf{y}$ 可能会遇到计算问题：

计算 $X^T X$ 需要 $O(n^2m)$ 的时间复杂度，求逆 $O(n^3)$ ，当 $n$ 很大时计算成本极高。
$X^T X$ 可能是病态矩阵（接近奇异矩阵），导致数值计算不稳定。

梯度下降是一种优化方法，通过迭代更新参数来逐步找到最优解。其基本思想是计算损失函数 $J(\mathbf{w})$ 对 $\mathbf{w}$ 的梯度，并沿着梯度的负方向更新参数：

$\mathbf{w} := \mathbf{w} - \alpha \nabla J(\mathbf{w})$

其中，学习率 $\alpha$ 控制每次更新的步长。

梯度计算如下：

$\nabla J(\mathbf{w}) = \frac{1}{m} X^T (X \mathbf{w} - \mathbf{y})$

因此，更新规则为：

$\mathbf{w} := \mathbf{w} - \frac{\alpha}{m} X^T (X \mathbf{w} - \mathbf{y})$

3. 直接法 vs 迭代法

方法	计算复杂度	适用情况
直接法（正规方程）	$O(n^2 m + n^3)$	适用于特征维度较小的数据集（如 $n < 10^3$ ）
迭代法（梯度下降）	每次迭代 $O (nm)$	适用于大规模数据集（如 $n > 10^4$ ）

在 $n$ 很大时，直接求逆是不可行的，因此我们使用梯度下降等迭代方法来近似求解。

4. 代码实现

4.1 直接法（正规方程）

import numpy as np

def least_squares_normal_eq(X, y):
    return np.linalg.inv(X.T @ X) @ X.T @ y

# 示例数据
X = np.array([[1, 1], [1, 2], [1, 3]])  # 增加一列1表示偏置项
y = np.array([[1], [2], [3]])

w = least_squares_normal_eq(X, y)
print("最优参数:", w)

4.2 迭代法（梯度下降）

import numpy as np

def gradient_descent(X, y, lr=0.01, epochs=1000):
    m, n = X.shape
    w = np.zeros((n, 1))  # 初始化权重
    for _ in range(epochs):
        gradient = (1/m) * X.T @ (X @ w - y)
        w -= lr * gradient
    return w

# 示例数据
X = np.array([[1, 1], [1, 2], [1, 3]])
y = np.array([[1], [2], [3]])

w = gradient_descent(X, y)
print("最优参数:", w)