机器学习数学基础专题【左扬精讲】— 最小二乘的梯度下降：数值求解、NumPy/PyTorch 实现与大规模实战
机器学习数学基础专题【左扬精讲】— 最小二乘的梯度下降：数值求解、NumPy/PyTorch 实现与大规模实战
本文是传播网络基础算法系列的第三篇，将深入探讨如何用梯度下降算法来求解最小二乘问题。前两篇文章分别介绍了最小二乘法的数学原理和梯度下降算法的核心思想，本篇将这座桥梁：讲解为什么可以用梯度下降求解最小二乘、如何实现完整的解决方案，以及在大规模场景下的工程实践。
最小二乘 + 梯度下降 数值稳定性 NumPy 实现 PyTorch 实现 分布式求解
学习重点
必须掌握：
- 最小二乘问题的梯度形式推导
- NumPy 实现批量/随机/小批量梯度下降
- PyTorch 实现自动求导的最小二乘求解器
- 收敛判断标准与最优学习率选择
需要了解：
- 共轭梯度法与高斯-牛顿法
- 分布式最小二乘求解
- 稀疏矩阵的特殊处理
一、从正规方程到梯度下降
Why — 有了正规方程，为什么还需要梯度下降？
正规方程 β) = (X)^TX)^-1X)^Ty) 给出了最小二乘问题的解析解，理论上一步到位。但在实际工程中，正规方程面临几个根本性困难：
- 计算复杂度：矩阵求逆是 O(p^3)，当特征数 p 达到百万级时不可行
- 存储问题：需要存储 X^TX 矩阵，大小是 p × p
- 数值稳定性：当 X^TX 接近奇异时，求逆会放大误差
- 增量更新：有新数据时无法增量更新，需要重新计算
How — 梯度下降如何解决这些问题
梯度下降的核心优势：
- 计算量可控：每次迭代只需 O(np)，不需矩阵求逆
- 内存友好：只需存储当前参数向量 O(p)
- 在线学习：可以流式处理新数据
- 易于分布式：可以天然并行化
代价是：只能找到近似解（但对于机器学习场景，足够精确即可），且需要调参（学习率、迭代次数）。
两种方法的对比速查表
| 方面 | 正规方程 | 梯度下降 |
|---|---|---|
| 复杂度（时间） | O(p^3 + np^2) | O(k · np) |
| 复杂度（空间） | O(p^2) | O(p) |
| 解的类型 | 解析解（精确） | 近似解（迭代） |
| 增量更新 | 不支持 | 支持 |
| 大规模数据 | 困难 | 容易 |
| 调参难度 | 无需调参 | 需要调参 |
其中 n 是样本数，p 是特征数，k 是迭代次数。
二、最小二乘梯度推导
What — 最小二乘问题的梯度形式
最小二乘问题的目标函数是残差平方和：
J(β) = (1)/(2n)Σ_i=1^n(y_i - x_i^Tβ)^2 = (1)/(2n)||y) - X)/(β)||^2
这里的 (1)/(2) 是为了抵消求导时产生的 2，(1)/(n) 是为了归一化（与样本数量无关）。
How — 梯度推导过程
对 J(β) 求偏导：
(∂ J)/(∂ β) = (∂)/(∂ β)[(1)/(2n)Σ_i=1^n(y_i - x_i^Tβ)^2]
对第 i 个样本：
(∂)/(∂ β)[(1)/(2)(y_i - x_i^Tβ)^2] = -(y_i - x_i^Tβ_x_i
对全部样本求和：
(∂ J)/(∂ β) = -(1)/(n)Σ_i=1^n(y_i - x_i^Tβ_x_i = (1)/(n_X^T(X)/(β) - y)
直观理解：梯度向量指向误差增大的方向
梯度 ∇ J = (1)/(n_X^T(X)/(β) - y) 可以这样理解：
- (X)/(β) - y) 是残差向量（预测值减去真实值）
- X)^T 将每个样本的残差投影回特征空间
- 整个梯度是每个样本残差的加权平均
当我们沿着 -∇ J 方向更新时，实际上是在减少最"突出"的误差方向。
向量化实现的关键
单个样本的梯度：g_i = -(y_i - x_i^Tβ_x_i
批量样本的梯度（向量化）：
∇ J = (1)/(n_X^T(X)/(β) - y)
这就是为什么 NumPy 的矩阵运算比 Python 循环快几十倍的原因：一条语句就能完成全部 n 个样本的梯度计算。
三、NumPy 实现：从零构建
从零实现一个完整的最小二乘梯度下降求解器
本节将用纯 NumPy 实现一个功能完整的线性回归求解器，涵盖批量、随机、小批量三种模式，并支持动量和学习率衰减。
核心类设计
import numpy as np
from typing import Tuple, Optional, Literal
class LeastSquaresGD:
"""
基于梯度下降的最小二乘求解器
支持三种模式：
- 'batch': 批量梯度下降
- 'stochastic': 随机梯度下降
- 'mini-batch': 小批量梯度下降
"""
def __init__(self,
learning_rate: float = 0.01,
n_iterations: int = 1000,
batch_size: Optional[int] = None,
momentum: float = 0.0,
lr_decay: float = 0.0,
verbose: bool = True):
"""
参数:
learning_rate: 学习率 α
n_iterations: 最大迭代次数
batch_size: batch大小 (None=全部, 1=随机, >1=小批量)
momentum: 动量系数 β (0=不使用)
lr_decay: 学习率衰减率
verbose: 是否打印训练过程
"""
self.lr = learning_rate
self.n_iterations = n_iterations
self.batch_size = batch_size
self.momentum = momentum
self.lr_decay = lr_decay
self.verbose = verbose
self.theta = None # 参数向量
self.velocity = None # 动量速度
self.loss_history = [] # 损失历史
def _compute_gradient(self, X: np.ndarray, y: np.ndarray,
theta: np.ndarray) -> np.ndarray:
"""计算梯度：∇J = (1/n) * X^T * (X*θ - y)"""
n = len(y)
residuals = X @ theta - y
gradient = (1/n) * (X.T @ residuals)
return gradient
def _compute_loss(self, X: np.ndarray, y: np.ndarray,
theta: np.ndarray) -> float:
"""计算损失：J = (1/2n) * ||X*θ - y||^2"""
n = len(y)
residuals = X @ theta - y
return (1/(2*n)) * np.sum(residuals ** 2)
def _get_batches(self, X: np.ndarray, y: np.ndarray) -> Tuple:
"""生成批次"""
n_samples = len(y)
indices = np.random.permutation(n_samples)
if self.batch_size is None:
# 批量模式：返回全部数据
yield X, y
elif self.batch_size == 1:
# 随机模式：每次一个样本
for idx in indices:
yield X[idx:idx+1], y[idx:idx+1]
else:
# 小批量模式
n_batches = int(np.ceil(n_samples / self.batch_size))
for b in range(n_batches):
start = b * self.batch_size
end = min((b + 1) * self.batch_size, n_samples)
batch_idx = indices[start:end]
yield X[batch_idx], y[batch_idx]
def fit(self, X: np.ndarray, y: np.ndarray) -> 'LeastSquaresGD':
"""
训练模型
参数:
X: 特征矩阵 (n_samples, n_features)
y: 目标向量 (n_samples,)
"""
n_samples, n_features = X.shape
# 初始化参数（零初始化）
self.theta = np.zeros(n_features)
self.velocity = np.zeros(n_features)
self.loss_history = []
if self.verbose:
print(f"开始训练: {n_samples} 样本, {n_features} 特征")
print(f"模式: {'批量' if self.batch_size is None else ('随机' if self.batch_size == 1 else '小批量')}")
print("-" * 50)
for iteration in range(self.n_iterations):
# 当前学习率（带衰减）
current_lr = self.lr / (1 + self.lr_decay * iteration)
# 遍历每个batch
for X_batch, y_batch in self._get_batches(X, y):
# 计算梯度
gradient = self._compute_gradient(X_batch, y_batch, self.theta)
# 动量更新
if self.momentum > 0:
self.velocity = (self.momentum * self.velocity +
(1 - self.momentum) * gradient)
gradient = self.velocity
# 参数更新
self.theta = self.theta - current_lr * gradient
# 记录损失（使用全量数据评估）
loss = self._compute_loss(X, y, self.theta)
self.loss_history.append(loss)
# 打印进度
if self.verbose and iteration % 100 == 0:
print(f"Iter {iteration:4d}: Loss = {loss:.8f}, LR = {current_lr:.6f}")
return self
def predict(self, X: np.ndarray) -> np.ndarray:
"""预测"""
return X @ self.theta
def score(self, X: np.ndarray, y: np.ndarray) -> float:
"""计算 R² 分数"""
y_pred = self.predict(X)
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
return 1 - ss_res / ss_tot
使用示例
# 生成模拟数据
np.random.seed(42)
n_samples = 1000
n_features = 10
# 生成随机特征矩阵
X = np.random.randn(n_samples, n_features)
# 生成真实参数
true_theta = np.random.randn(n_features)
# 生成目标变量（带噪声）
y = X @ true_theta + np.random.randn(n_samples) * 0.5
# 划分训练集和测试集
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# 1. 批量梯度下降
print("=" * 50)
print("批量梯度下降 (BGD)")
print("=" * 50)
model_bgd = LeastSquaresGD(
learning_rate=0.1,
n_iterations=1000,
batch_size=None # 批量模式
).fit(X_train, y_train)
print(f"测试集 R²: {model_bgd.score(X_test, y_test):.4f}")
# 2. 随机梯度下降
print("\n" + "=" * 50)
print("随机梯度下降 (SGD)")
print("=" * 50)
model_sgd = LeastSquaresGD(
learning_rate=0.01,
n_iterations=1000,
batch_size=1 # 随机模式
).fit(X_train, y_train)
print(f"测试集 R²: {model_sgd.score(X_test, y_test):.4f}")
# 3. 小批量梯度下降（带动量）
print("\n" + "=" * 50)
print("小批量梯度下降 + 动量 (MBGD + Momentum)")
print("=" * 50)
model_mbgd = LeastSquaresGD(
learning_rate=0.1,
n_iterations=1000,
batch_size=32,
momentum=0.9, # 动量
lr_decay=0.001 # 学习率衰减
).fit(X_train, y_train)
print(f"测试集 R²: {model_mbgd.score(X_test, y_test):.4f}")
# 4. 与正规方程对比
print("\n" + "=" * 50)
print("正规方程解（ground truth）")
print("=" * 50)
theta_normal = np.linalg.lstsq(X_train, y_train, rcond=None)[0]
y_pred_normal = X_test @ theta_normal
ss_res = np.sum((y_test - y_pred_normal) ** 2)
ss_tot = np.sum((y_test - np.mean(y_test)) ** 2)
r2_normal = 1 - ss_res / ss_tot
print(f"正规方程 R²: {r2_normal:.4f}")
收敛可视化
import matplotlib.pyplot as plt
# 可视化收敛过程
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左图：损失曲线对比
ax1 = axes[0]
ax1.plot(model_bgd.loss_history, label='BGD', alpha=0.8)
ax1.plot(model_sgd.loss_history, label='SGD', alpha=0.8)
ax1.plot(model_mbgd.loss_history, label='MBGD+Momentum', alpha=0.8)
ax1.set_xlabel('Iteration')
ax1.set_ylabel('Loss')
ax1.set_title('Loss Convergence Comparison')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 右图：预测结果对比
ax2 = axes[1]
y_pred_bgd = model_bgd.predict(X_test)
y_pred_mbgd = model_mbgd.predict(X_test)
ax2.scatter(y_test, y_pred_mbgd, alpha=0.5, label='MBGD+Momentum', s=10)
ax2.scatter(y_test, y_pred_normal, alpha=0.5, label='Normal Equation', s=10)
ax2.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()],
'r--', label='Perfect Prediction')
ax2.set_xlabel('True Values')
ax2.set_ylabel('Predictions')
ax2.set_title('Prediction Comparison')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
四、PyTorch 实现：自动求导
PyTorch 的优势：自动微分
PyTorch 提供了自动微分（autograd）功能，我们只需定义前向传播，梯度计算由框架自动完成。这不仅简化了代码，还支持复杂的非线性模型。
方法一：使用 autograd 自动求导
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
class LinearRegressionTorch:
"""
PyTorch 实现的线性回归（使用 autograd）
"""
def __init__(self, input_dim: int, learning_rate: float = 0.01):
"""
参数:
input_dim: 特征维度
learning_rate: 学习率
"""
self.input_dim = input_dim
# 定义模型（线性层）
self.model = nn.Linear(input_dim, 1)
# 定义损失函数：MSE
self.criterion = nn.MSELoss(reduction='mean')
# 定义优化器：SGD
self.optimizer = optim.SGD(self.model.parameters(), lr=learning_rate)
# 记录训练历史
self.loss_history = []
def fit(self, X: np.ndarray, y: np.ndarray,
n_epochs: int = 1000, batch_size: int = 32,
verbose: bool = True) -> 'LinearRegressionTorch':
"""
训练模型
参数:
X: 特征矩阵 (n_samples, n_features)
y: 目标向量 (n_samples,)
n_epochs: 训练轮数
batch_size: batch大小
verbose: 是否打印进度
"""
# 转换为 PyTorch 张量
X_tensor = torch.tensor(X, dtype=torch.float32)
y_tensor = torch.tensor(y, dtype=torch.float32).unsqueeze(1)
# 创建 DataLoader
dataset = TensorDataset(X_tensor, y_tensor)
dataloader = DataLoader(dataset, batch_size=batch_size, shuffle=True)
self.model.train() # 设置为训练模式
for epoch in range(n_epochs):
epoch_loss = 0.0
for batch_X, batch_y in dataloader:
# 前向传播
predictions = self.model(batch_X)
# 计算损失
loss = self.criterion(predictions, batch_y)
# 反向传播
self.optimizer.zero_grad() # 清零梯度
loss.backward() # 自动计算梯度
# 更新参数
self.optimizer.step()
epoch_loss += loss.item() * len(batch_X)
# 记录平均损失
avg_loss = epoch_loss / len(dataset)
self.loss_history.append(avg_loss)
if verbose and epoch % 100 == 0:
print(f"Epoch {epoch:4d}: Loss = {avg_loss:.8f}")
return self
def predict(self, X: np.ndarray) -> np.ndarray:
"""预测"""
self.model.eval() # 设置为评估模式
with torch.no_grad():
X_tensor = torch.tensor(X, dtype=torch.float32)
predictions = self.model(X_tensor).numpy().squeeze()
return predictions
def get_weights(self) -> np.ndarray:
"""获取模型参数"""
return self.model.weight.data.numpy().squeeze()
def get_bias(self) -> float:
"""获取偏置"""
return self.model.bias.data.numpy().squeeze()
# 使用示例
print("=" * 50)
print("PyTorch 自动求导实现")
print("=" * 50)
# 数据准备
np.random.seed(42)
torch.manual_seed(42)
n_samples = 1000
n_features = 5
X = np.random.randn(n_samples, n_features)
true_theta = np.random.randn(n_features)
y = X @ true_theta + np.random.randn(n_samples) * 0.3
# 训练
model_torch = LinearRegressionTorch(
input_dim=n_features,
learning_rate=0.1
)
model_torch.fit(X, y, n_epochs=1000, batch_size=32)
# 评估
y_pred = model_torch.predict(X)
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r2 = 1 - ss_res / ss_tot
print(f"\nR² Score: {r2:.4f}")
print(f"学习到的权重: {model_torch.get_weights()}")
print(f"真实权重: {true_theta}")
方法二：手动实现 autograd 函数
对于需要更精细控制的场景，可以自定义 autograd 函数：
import torch
from torch.autograd import Function
class LeastSquaresLoss(Function):
"""
自定义最小二乘损失 autograd 函数
展示梯度计算的内部原理
"""
@staticmethod
def forward(ctx, X: torch.Tensor, y: torch.Tensor,
theta: torch.Tensor) -> torch.Tensor:
"""
前向传播：计算 J = (1/2n) * ||X*θ - y||^2
"""
n = y.shape[0]
residual = X @ theta - y
loss = (1/(2*n)) * torch.sum(residual ** 2)
# 保存反向传播需要的变量
ctx.save_for_backward(X, y, residual)
ctx.n = n
return loss
@staticmethod
def backward(ctx, grad_output: torch.Tensor):
"""
反向传播：计算梯度 ∂J/∂θ = (1/n) * X^T * (X*θ - y)
"""
X, y, residual = ctx.saved_tensors
n = ctx.n
# 梯度计算（与手动推导一致）
gradient = (1/n) * (X.T @ residual)
return None, None, gradient # X和y不需要梯度
def train_manual_autograd(X: np.ndarray, y: np.ndarray,
lr: float = 0.1,
n_epochs: int = 1000) -> np.ndarray:
"""
使用自定义 autograd 函数训练
"""
# 转换为 PyTorch 张量（需要梯度）
X_t = torch.tensor(X, dtype=torch.float32, requires_grad=False)
y_t = torch.tensor(y, dtype=torch.float32, requires_grad=False)
theta = torch.zeros(X.shape[1], dtype=torch.float32, requires_grad=True)
optimizer = torch.optim.SGD([theta], lr=lr)
for epoch in range(n_epochs):
# 前向传播（使用自定义损失）
loss = LeastSquaresLoss.apply(X_t, y_t, theta)
# 反向传播
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 200 == 0:
print(f"Epoch {epoch}: Loss = {loss.item():.8f}")
return theta.detach().numpy()
# 测试
print("=" * 50)
print("自定义 autograd 函数实现")
print("=" * 50)
theta_learned = train_manual_autograd(X, y, lr=0.1, n_epochs=1000)
print(f"\n学习到的参数: {theta_learned}")
print(f"真实参数: {true_theta}")
PyTorch vs NumPy 选择建议
- 学习目的：使用 NumPy，理解底层原理
- 简单线性模型：NumPy 足够
- 复杂模型（神经网络）：必须用 PyTorch
- 需要 GPU 加速：必须用 PyTorch
- 生产环境：考虑 scikit-learn 或专业库
五、收敛判断与最优学习率
收敛判断的三种标准
- 损失阈值：当损失低于某个阈值时停止
- 梯度阈值：当梯度范数足够小时停止
- 相对变化：当连续两次损失的相对变化足够小时停止
实现早停机制
class EarlyStoppingGD:
"""
带早停的梯度下降求解器
"""
def __init__(self,
patience: int = 10,
min_delta: float = 1e-6,
restore_best: bool = True):
"""
参数:
patience: 容忍多少个epoch没有改善
min_delta: 认为有改善的最小变化量
restore_best: 是否恢复到最佳参数
"""
self.patience = patience
self.min_delta = min_delta
self.restore_best = restore_best
self.best_loss = float('inf')
self.best_theta = None
self.counter = 0
self.should_stop = False
def check(self, loss: float, theta: np.ndarray) -> bool:
"""
检查是否应该停止训练
返回:
True 表示应该停止，False 表示继续训练
"""
if loss < self.best_loss - self.min_delta:
# 有改善
self.best_loss = loss
self.best_theta = theta.copy()
self.counter = 0
else:
# 没有改善
self.counter += 1
if self.counter >= self.patience:
self.should_stop = True
return True
return False
def get_best_params(self) -> np.ndarray:
"""获取最佳参数"""
if self.restore_best and self.best_theta is not None:
return self.best_theta
return None
def fit_with_early_stopping(X: np.ndarray, y: np.ndarray,
lr: float = 0.1,
max_epochs: int = 10000,
patience: int = 20) -> Tuple[np.ndarray, list]:
"""
带早停的梯度下降训练
"""
n_features = X.shape[1]
theta = np.zeros(n_features)
loss_history = []
early_stopping = EarlyStoppingGD(patience=patience)
for epoch in range(max_epochs):
# 计算梯度
gradient = (1/len(y)) * (X.T @ (X @ theta - y))
# 更新参数
theta = theta - lr * gradient
# 计算损失
loss = (1/(2*len(y))) * np.sum((X @ theta - y) ** 2)
loss_history.append(loss)
# 检查早停
if early_stopping.check(loss, theta):
print(f"早停触发！Epoch {epoch}, Best Loss: {early_stopping.best_loss:.8f}")
break
if epoch % 500 == 0:
print(f"Epoch {epoch}: Loss = {loss:.8f}")
best_theta = early_stopping.get_best_params()
return best_theta, loss_history
# 测试早停
print("=" * 50)
print("早停机制测试")
print("=" * 50)
theta_early, history_early = fit_with_early_stopping(
X, y, lr=0.1, patience=20
)
最优学习率的选择
使用 Learning Rate Finder 方法：
def find_best_learning_rate(X: np.ndarray, y: np.ndarray,
lr_start: float = 1e-8,
lr_end: float = 10,
n_steps: int = 100) -> Tuple[float, list]:
"""
学习率finder：找到最优学习率
方法：逐渐增大学习率，记录损失变化
"""
n_features = X.shape[1]
losses = []
learning_rates = []
# 指数级增长的学习率
factor = (lr_end / lr_start) ** (1 / n_steps)
theta = np.zeros(n_features)
for step in range(n_steps):
lr = lr_start * (factor ** step)
learning_rates.append(lr)
# 计算梯度并更新
gradient = (1/len(y)) * (X.T @ (X @ theta - y))
theta = theta - lr * gradient
# 计算损失
loss = (1/(2*len(y))) * np.sum((X @ theta - y) ** 2)
losses.append(loss)
# 如果损失爆炸，停止
if loss > 10 * losses[0] or np.isnan(loss):
break
# 找到损失下降最快的点
best_idx = np.argmin(losses)
best_lr = learning_rates[best_idx]
return best_lr, list(zip(learning_rates, losses))
# 使用学习率finder
print("=" * 50)
print("学习率 Finder")
print("=" * 50)
best_lr, lr_results = find_best_learning_rate(X, y)
print(f"推荐学习率: {best_lr:.6f}")
# 可视化
lrs, losses = zip(*lr_results)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(lrs, losses)
plt.xscale('log')
plt.xlabel('Learning Rate')
plt.ylabel('Loss')
plt.title('Learning Rate Finder')
plt.axvline(x=best_lr, color='r', linestyle='--', label=f'Best LR = {best_lr:.6f}')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
plt.plot(losses)
plt.xlabel('Steps')
plt.ylabel('Loss')
plt.title('Loss vs Steps')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
六、大规模场景下的工程实践
大规模数据的挑战
当数据量达到百万级甚至更大时，单机训练会遇到：
- 内存不够：无法一次性加载所有数据
- 计算太慢：单次迭代需要处理所有样本
- 精度问题：浮点数精度随数据量增大而降低
解决方案一：流式训练（Online Learning）
数据流式到达，每次只用最新的样本更新：
class OnlineLeastSquaresGD:
"""
在线最小二乘梯度下降
特点：每次只处理一个样本，适合流式数据
"""
def __init__(self, learning_rate: float = 0.01,
forgetting_factor: float = 0.999):
"""
参数:
learning_rate: 学习率
forgetting_factor: 遗忘因子（给历史数据较小的权重）
"""
self.lr = learning_rate
self.ff = forgetting_factor
self.theta = None
self.n_updates = 0
def partial_fit(self, x: np.ndarray, y: float) -> np.ndarray:
"""
部分拟合：用单个样本更新模型
参数:
x: 单个样本特征 (n_features,)
y: 单个样本标签
"""
x = x.reshape(-1, 1) # (n_features, 1)
# 初始化
if self.theta is None:
self.theta = np.zeros(x.shape[0])
# 计算预测和残差
y_pred = self.theta @ x.squeeze()
residual = y_pred - y
# 计算梯度：∇J = (y_pred - y) * x
gradient = residual * x.squeeze()
# 更新参数（带遗忘）
effective_lr = self.lr * (self.ff ** self.n_updates)
self.theta = self.theta - effective_lr * gradient
self.n_updates += 1
return self.theta.copy()
def predict(self, X: np.ndarray) -> np.ndarray:
"""预测"""
return X @ self.theta
# 模拟流式数据
print("=" * 50)
print("在线学习（流式数据）")
print("=" * 50)
online_model = OnlineLeastSquaresGD(learning_rate=0.01)
# 模拟流式到达的数据
np.random.seed(42)
for i in range(1000):
x_stream = np.random.randn(5)
y_stream = x_stream @ true_theta + np.random.randn() * 0.3
online_model.partial_fit(x_stream, y_stream)
if i % 200 == 0:
print(f"已处理 {i} 个样本...")
print(f"\n最终参数: {online_model.theta}")
print(f"真实参数: {true_theta}")
解决方案二：分布式训练
当单机无法处理时，可以使用数据并行：
def distributed_gradient_descent(X: np.ndarray, y: np.ndarray,
theta: np.ndarray,
lr: float,
n_workers: int = 4) -> np.ndarray:
"""
简化的分布式梯度下降
思路：将数据分片到多个worker，每个worker计算本地梯度，
然后汇总更新（类似Parameter Server架构）
参数:
X: 特征矩阵
y: 目标向量
theta: 当前参数
lr: 学习率
n_workers: worker数量
"""
n_samples = len(y)
samples_per_worker = n_samples // n_workers
# 收集所有worker的梯度
all_gradients = []
for i in range(n_workers):
start = i * samples_per_worker
end = start + samples_per_worker if i < n_workers - 1 else n_samples
X_chunk = X[start:end]
y_chunk = y[start:end]
# 计算该worker负责的梯度
gradient = (1/n_samples) * (X_chunk.T @ (X_chunk @ theta - y_chunk))
all_gradients.append(gradient)
# 聚合梯度（简单平均）
aggregated_gradient = np.mean(all_gradients, axis=0)
# 更新参数
theta = theta - lr * aggregated_gradient
return theta
# 使用分布式训练
print("=" * 50)
print("分布式梯度下降")
print("=" * 50)
n_samples = 10000
n_features = 10
X_dist = np.random.randn(n_samples, n_features)
y_dist = X_dist @ true_theta + np.random.randn(n_samples) * 0.5
theta_dist = np.zeros(n_features)
for epoch in range(1000):
theta_dist = distributed_gradient_descent(X_dist, y_dist, theta_dist,
lr=0.1, n_workers=4)
if epoch % 200 == 0:
loss = (1/(2*n_samples)) * np.sum((X_dist @ theta_dist - y_dist) ** 2)
print(f"Epoch {epoch}: Loss = {loss:.6f}")
print(f"\n分布式结果: {theta_dist}")
print(f"真实参数: {true_theta}")
实际生产环境的选择
- 中小规模数据（<10万样本）：直接用 sklearn 的 LinearRegression 或 SGDRegressor
- 大规模数据（10万-1亿样本）：使用 Spark MLlib 的 LinearRegressionWithSGD
- 超大规模（>1亿样本）：使用 Vowpal Wabbit 或专用分布式机器学习平台
七、完整实战案例
实战：用梯度下降求解实际问题
本节用一个完整的案例串联所有知识点：预测房价。
数据集说明
使用模拟的房价数据集，包含 5 个特征：
- 面积（平方米）
- 卧室数量
- 房龄（年）
- 距市中心距离（公里）
- 绿化率
import numpy as np
import matplotlib.pyplot as plt
# 生成模拟房价数据
np.random.seed(42)
n_samples = 500
# 特征：面积、卧室数、房龄、距市中心距离、绿化率
X = np.zeros((n_samples, 5))
X[:, 0] = np.random.uniform(50, 200, n_samples) # 面积
X[:, 1] = np.random.randint(1, 6, n_samples) # 卧室数
X[:, 2] = np.random.uniform(0, 30, n_samples) # 房龄
X[:, 3] = np.random.uniform(1, 20, n_samples) # 距市中心
X[:, 4] = np.random.uniform(0.1, 0.8, n_samples) # 绿化率
# 真实参数
true_theta = np.array([2.5, 5.0, -1.0, -3.0, 10.0]) # 每单位特征的价格影响
base_price = 50 # 基础价格（万元）
# 生成房价（带噪声）
y = (base_price + X @ true_theta +
np.random.randn(n_samples) * 5) # 万元
print("=" * 60)
print("房价预测实战")
print("=" * 60)
print(f"数据集大小: {n_samples} 样本, {X.shape[1]} 特征")
print(f"\n特征统计:")
feature_names = ['面积(m²)', '卧室数', '房龄(年)', '距市中心(km)', '绿化率']
for i, name in enumerate(feature_names):
print(f" {name}: 均值={X[:,i].mean():.2f}, 标准差={X[:,i].std():.2f}")
print(f"\n房价范围: {y.min():.1f} - {y.max():.1f} 万元")
print(f"房价均值: {y.mean():.1f} 万元")
# 特征标准化（非常重要！）
X_mean = X.mean(axis=0)
X_std = X.std(axis=0)
X_scaled = (X - X_mean) / X_std
# 使用我们之前定义的 LeastSquaresGD 类
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
X_scaled, y, test_size=0.2, random_state=42
)
# 训练多个模型对比
results = {}
# 1. 批量梯度下降
model1 = LeastSquaresGD(learning_rate=0.5, n_iterations=2000)
model1.fit(X_train, y_train)
results['BGD'] = {
'model': model1,
'train_r2': model1.score(X_train, y_train),
'test_r2': model1.score(X_test, y_test)
}
# 2. 小批量 + 动量
model2 = LeastSquaresGD(
learning_rate=0.5,
n_iterations=2000,
batch_size=32,
momentum=0.9
)
model2.fit(X_train, y_train)
results['MBGD+Momentum'] = {
'model': model2,
'train_r2': model2.score(X_train, y_train),
'test_r2': model2.score(X_test, y_test)
}
# 3. 小批量 + 动量 + 学习率衰减
model3 = LeastSquaresGD(
learning_rate=0.5,
n_iterations=2000,
batch_size=32,
momentum=0.9,
lr_decay=0.01
)
model3.fit(X_train, y_train)
results['MBGD+Momentum+Decay'] = {
'model': model3,
'train_r2': model3.score(X_train, y_train),
'test_r2': model3.score(X_test, y_test)
}
# 4. 正规方程对比
theta_normal = np.linalg.lstsq(X_train, y_train, rcond=None)[0]
y_pred_train = X_train @ theta_normal
y_pred_test = X_test @ theta_normal
ss_res_train = np.sum((y_train - y_pred_train) ** 2)
ss_tot_train = np.sum((y_train - np.mean(y_train)) ** 2)
ss_res_test = np.sum((y_test - y_pred_test) ** 2)
ss_tot_test = np.sum((y_test - np.mean(y_test)) ** 2)
results['Normal Eq.'] = {
'model': None,
'train_r2': 1 - ss_res_train / ss_tot_train,
'test_r2': 1 - ss_res_test / ss_tot_test
}
print("\n" + "=" * 60)
print("模型对比")
print("=" * 60)
print(f"{'模型':<20} {'训练R²':>12} {'测试R²':>12}")
print("-" * 44)
for name, res in results.items():
print(f"{name:<20} {res['train_r2']:>12.4f} {res['test_r2']:>12.4f}")
# 可视化结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 1. 损失曲线
ax1 = axes[0, 0]
for name, res in results.items():
if res['model'] is not None:
ax1.plot(res['model'].loss_history, label=name, alpha=0.8)
ax1.set_xlabel('Iteration')
ax1.set_ylabel('Loss')
ax1.set_title('Loss Convergence')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. 预测 vs 真实
ax2 = axes[0, 1]
best_model_name = 'MBGD+Momentum'
y_pred_best = results[best_model_name]['model'].predict(X_test)
ax2.scatter(y_test, y_pred_best, alpha=0.6, s=30)
ax2.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()],
'r--', lw=2, label='Perfect Prediction')
ax2.set_xlabel('True Price (万元)')
ax2.set_ylabel('Predicted Price (万元)')
ax2.set_title(f'Prediction vs Truth ({best_model_name})')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 3. 残差分布
ax3 = axes[1, 0]
residuals = y_test - y_pred_best
ax3.hist(residuals, bins=30, edgecolor='black', alpha=0.7)
ax3.axvline(x=0, color='r', linestyle='--')
ax3.set_xlabel('Residual (万元)')
ax3.set_ylabel('Frequency')
ax3.set_title('Residual Distribution')
ax3.grid(True, alpha=0.3)
# 4. 特征重要性
ax4 = axes[1, 1]
final_theta = results[best_model_name]['model'].theta
colors = ['green' if c > 0 else 'red' for c in final_theta]
ax4.barh(feature_names, final_theta, color=colors, alpha=0.7)
ax4.axvline(x=0, color='black', linestyle='-')
ax4.set_xlabel('Coefficient (标准化后的影响)')
ax4.set_title('Feature Importance (Standardized Coefficients)')
ax4.grid(True, alpha=0.3, axis='x')
plt.tight_layout()
plt.show()
# 输出预测示例
print("\n" + "=" * 60)
print("预测示例")
print("=" * 60)
print(f"{'实际房价':>12} {'预测房价':>12} {'误差':>12}")
print("-" * 36)
for i in range(10):
print(f"{y_test[i]:>12.2f} {y_pred_best[i]:>12.2f} {y_test[i]-y_pred_best[i]:>+12.2f}")
实战总结
- 特征标准化：必须步骤，否则不同量纲的特征会导致收敛困难
- 学习率选择：通过 Learning Rate Finder 或实验确定
- 早停机制：防止过拟合，节省训练时间
- 动量+衰减：加速收敛，提高稳定性
FAQ 常见问题解答
Q1. 最小二乘梯度下降和正规方程哪个更精确？
正规方程给出解析解（理论上精确），梯度下降给出近似解。但梯度下降的"不精确"在实际中通常不是问题，因为数据本身有噪声，解到一定程度就足够好了。
Q2. 为什么需要特征标准化？
两个原因：① 不同量纲的特征会让损失函数等高线是椭圆而非圆形，梯度下降路径曲折；② 学习率的选取会很困难，某些方向步伐太大而某些方向太小。
Q3. 什么时候用 L1 正则化而不是 L2？
当需要稀疏解（特征选择）时用 L1（Lasso）。L1 正则化会让不重要的特征系数变为零，实现自动特征选择。L2 正则化（Ridge）只是收缩系数，不会变零。
Q4. 如何判断模型是否过拟合？
训练集 R² 高、测试集 R² 低。如果训练集和测试集差距很大，说明过拟合。可以使用早停、正则化、增加数据量或减少特征来解决。
Q5. 为什么随机梯度下降比批量梯度下降收敛更快？
因为单位计算时间内，随机梯度下降可以迭代更多次。虽然单次方向不如批量准确，但"更多次迭代"弥补了方向的不准确性，最终反而更快。
Q6. 学习率衰减和动量可以同时用吗？
可以，而且通常效果更好。动量帮助跨越平坦区域和减少振荡，学习率衰减帮助在接近最优时精细调整。两者结合是深度学习的标准做法。
Q7. 残差分析有什么作用？
诊断模型假设是否成立。如果残差分布偏离正态、存在系统性模式或异常值，说明模型假设有问题，需要改进模型或数据处理。
Q8. 如何处理特征共线性？
三种方法：① 使用岭回归（L2 正则化）；② 删除高度相关的特征；③ 使用 PCA 降维。正规方程遇到共线性时矩阵接近奇异，梯度下降没有这个问题。
Q9. 梯度下降一定会收敛吗？
对于凸函数，使用合适的学习率，梯度下降一定收敛。对于非凸函数（如神经网络），只能保证收敛到局部最优或鞍点。
Q10. 为什么神经网络训练用随机初始化而不是零初始化？
零初始化会导致对称性问题。如果所有权重都是零，所有神经元的梯度相同，无法学习不同特征。随机初始化打破对称性，让不同神经元学习不同特征。
Q11. 什么是二阶优化？为什么不常用？
二阶优化使用 Hessian 矩阵（如牛顿法）。收敛速度比一阶快（理论上二次收敛），但计算 Hessian 及其逆是 O(p^2) 到 O(p^3)，内存和计算量都太大。拟牛顿法（如 L-BFGS）试图近似 Hessian，但在大规模深度学习中仍然不实用。
Q12. 在线学习和批量学习的区别是什么？
批量学习需要所有数据一次性加载，训练后模型固定；在线学习数据流式到达，每次用新样本更新模型，无需存储历史数据。在线学习适合数据持续产生的场景。
Q13. 什么是 federated learning？
联邦学习。数据分散在多个客户端（如手机），每个客户端本地训练并只上传梯度，服务器聚合梯度更新全局模型。隐私保护场景下很有用。
Q14. 梯度爆炸时如何诊断？
观察损失值。如果损失突然变成 NaN 或 inf，或者参数值突然变得非常大（超出 float 表示范围），说明梯度爆炸。可以使用梯度裁剪来缓解。
Q15. 为什么深度学习用小批量而不是大批量？
两个原因：①大批量内存消耗大；②大批量泛化能力通常不如小批量（近期研究）。小批量（32/64/128）是精度和效率的平衡点。
Q16. 如何选择正则化参数？
交叉验证。尝试不同的正则化强度（如 [10^-4, 10^-3, ..., 10^3]），选择验证集表现最好的。如果训练集和验证集差距随正则化增大而减小，说明正则化在起作用。
Q17. 什么是最小二乘的递归形式？
Recursive Least Squares (RLS)。每次新增一个样本时，不需要重新计算整个矩阵求逆，而是用递推公式更新参数。适合在线学习场景。
Q18. 共轭梯度法和梯度下降有什么区别？
共轭梯度法在 n 步内就能找到 n 维二次函数的精确最小值，而普通梯度下降可能需要无限步。共轭梯度法介于梯度下降和二阶方法之间。
Q19. 如何处理异方差问题？
使用加权最小二乘法。给方差大的样本分配较小的权重，让模型更关注方差小的样本。权重可以选择方差的倒数。
Q20. 什么情况下梯度下降比正规方程更好？
大规模数据、特征数量大、在线学习、需要正则化。当数据量达到百万级或特征数达到十万级时，正规方程的 O(p^3) 复杂度不可接受，梯度下降成为唯一选择。
三篇文章核心知识点回顾
- 最小二乘法：最小化残差平方和，正规方程给出解析解
- 梯度下降：迭代优化，沿负梯度方向更新参数
- 两者结合：用梯度下降求解最小二乘，适合大规模场景
- 关键技巧：特征标准化、学习率选择、早停机制、动量优化
Roadmap 后续预告
机器学习数学基础专题到此告一段落。我们从最小二乘法的理论基础出发，经过梯度下降算法的核心原理，最终掌握了用梯度下降求解最小二乘的完整方法。后续我们将探索更深入的主题：
- 神经网络基础：从感知机到多层神经网络
- 反向传播算法：链式法则的工程实现
- 卷积神经网络：图像处理的利器
- 循环神经网络：序列数据的处理
浙公网安备 33010602011771号