最小二乘法 Python 3.12 实现：从公式推导到 10 行代码拟合直线

在数据分析与机器学习领域，线性回归是最基础也最重要的算法之一。而最小二乘法作为实现线性回归的核心数学工具，其价值不仅在于理论上的优雅，更在于实践中的高效。本文将带你从数学公式推导开始，最终用 Python 3.12 实现一个仅需 10 行代码的最小二乘法直线拟合。

1. 最小二乘法的数学本质

当我们面对一组散点数据时，最小二乘法的目标是找到一条直线，使得所有数据点到这条直线的垂直距离（残差）的平方和最小。这种"平方和最小"的思想，在数学上体现为优化问题。

对于一元线性回归模型 y = ax + b，最小二乘法的核心公式推导如下：

1. 定义残差平方和（SSE）：

   SSE = Σ(y_i - (ax_i + b))²
   

2. 对参数 a 和 b 分别求偏导并令其为零：

   ∂SSE/∂a = -2Σx_i(y_i - ax_i - b) = 0
   ∂SSE/∂b = -2Σ(y_i - ax_i - b) = 0
   

3. 解这个方程组得到闭式解：

   a = (nΣx_i y_i - Σx_i Σy_i) / (nΣx_i² - (Σx_i)²)
   b = (Σy_i - aΣx_i) / n
   

这个推导过程展示了最小二乘法如何将几何直觉转化为可计算的数学表达式。值得注意的是，这种方法得到的解是全局最优的，这是最小二乘法相比梯度下降等迭代方法的优势所在。

2. Python 实现前的准备工作

在开始编码前，我们需要明确几个关键点：

• 数据表示 ：使用 NumPy 数组存储特征 x 和标签 y
• 计算效率 ：利用向量化运算替代循环
• 数值稳定性 ：处理除零等边界情况

以下是实现所需的基础计算步骤：

计算项	数学表达式	Python实现
样本量	n	len(x)
x求和	Σx	np.sum(x)
y求和	Σy	np.sum(y)
xy乘积和	Σxy	np.dot(x, y)
x平方和	Σx²	np.dot(x, x)

3. 10 行核心实现代码

基于上述数学原理，以下是完整的 Python 实现：

import numpy as np

def least_squares_fit(x, y):
    n = len(x)
    sum_x, sum_y = np.sum(x), np.sum(y)
    sum_xy = np.dot(x, y)
    sum_xx = np.dot(x, x)
    
    denominator = n * sum_xx - sum_x**2
    a = (n * sum_xy - sum_x * sum_y) / denominator
    b = (sum_y * sum_xx - sum_x * sum_xy) / denominator
    
    return a, b

这段代码的精妙之处在于：

1. 完全基于 NumPy 的向量化运算，避免了低效的 Python 循环
2. 严格遵循数学推导公式，保证了计算精度
3. 代码行数虽少，但包含了所有关键计算步骤

4. 代码验证与可视化

为了验证我们的实现是否正确，让我们用实际数据测试：

# 测试数据
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 4, 5, 4, 5])

# 计算回归系数
a, b = least_squares_fit(x, y)
print(f"斜率a: {a:.4f}, 截距b: {b:.4f}")

# 可视化
import matplotlib.pyplot as plt
plt.scatter(x, y, color='blue', label='原始数据')
plt.plot(x, a*x + b, color='red', label='拟合直线')
plt.legend()
plt.show()

执行结果应该显示：

斜率a: 0.6000, 截距b: 2.2000

为了进一步验证，我们可以与 scikit-learn 的结果对比：

from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(x.reshape(-1,1), y)
print(f"Sklearn结果 - 斜率: {model.coef_[0]:.4f}, 截距: {model.intercept_:.4f}")

两者结果应该完全一致，这验证了我们手动实现的正确性。

5. 性能优化与边界处理

虽然上述实现已经相当简洁，但在实际应用中我们还需要考虑一些优化和健壮性处理：

1. 数值稳定性 ：当分母接近零时的处理
2. 内存效率 ：对于超大数组的优化
3. 类型检查 ：确保输入数据的有效性

改进后的版本如下：

def robust_least_squares(x, y, eps=1e-10):
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    
    n = len(x)
    if len(y) != n:
        raise ValueError("x和y长度必须相同")
    
    sum_x, sum_y = np.sum(x), np.sum(y)
    sum_xy = np.dot(x, y)
    sum_xx = np.dot(x, x)
    
    denominator = n * sum_xx - sum_x**2
    if abs(denominator) < eps:
        raise ValueError("数据共线性过高，无法计算")
    
    a = (n * sum_xy - sum_x * sum_y) / denominator
    b = (sum_y - a * sum_x) / n
    
    return a, b

这个增强版增加了输入验证、数值稳定性检查和浮点精度控制，更适合生产环境使用。

6. 数学原理与代码的对应关系

为了更深入理解代码背后的数学，让我们分解关键步骤：

1. 向量点积与求和 ：

   np.dot(x, y)  # 对应数学中的 Σx_i y_i
   

   这种向量化运算不仅代码简洁，而且底层使用优化过的BLAS库，速度比手动循环快数十倍

2. 分母计算 ：

   n * sum_xx - sum_x**2  # 对应 nΣx_i² - (Σx_i)²
   

   这个值实际上就是 x 的方差乘以 n²，分母过小意味着 x 值几乎无变化

3. 系数求解 ：

   a = (n * sum_xy - sum_x * sum_y) / denominator
   

   这正是协方差与方差的比值，体现了 x 和 y 的线性关系强度

7. 实际应用案例

最小二乘法在现实中有广泛应用，例如：

• 经济学 ：预测GDP增长与失业率的关系
• 工程学 ：传感器数据的校准曲线拟合
• 生物学 ：药物剂量与反应的关系建模

以下是一个房价预测的示例：

# 房屋面积(平米)与价格(万元)
areas = np.array([50, 70, 90, 110, 130])
prices = np.array([180, 230, 280, 320, 400])

a, b = least_squares_fit(areas, prices)
print(f"每平米价格: {a:.2f}万元, 基础价格: {b:.2f}万元")

# 预测120平米房屋价格
predicted = a * 120 + b
print(f"预测价格: {predicted:.2f}万元")

8. 扩展思考

虽然我们的实现已经相当完善，但仍有值得探讨的扩展方向：

1. 多元线性回归 ：如何扩展到多个特征的情况
2. 正则化 ：加入L1/L2正则项防止过拟合
3. 加权最小二乘 ：对不同数据点赋予不同权重
4. 非线性扩展 ：通过基函数变换处理非线性关系

这些高级话题都建立在本文介绍的最小二乘基础之上，理解了这个核心算法，后续学习会更加顺畅。

最小二乘法的价值不仅在于它本身，更在于它体现的优化思想。从高斯发现它用于天体轨道计算，到如今成为机器学习的基石，这个诞生于18世纪的算法依然闪耀着智慧的光芒。