用Python+NumPy实战PGD优化算法：从零实现带约束的线性回归

在机器学习的世界里，优化算法就像探险家的指南针，指引我们穿越复杂的参数空间寻找最优解。Projected Gradient Descent（PGD）作为梯度下降家族的重要成员，因其独特的"投影"机制而成为处理带约束优化问题的利器。本文将带你用Python和NumPy从零实现PGD算法，通过一个具体的线性回归案例，让你不仅理解数学公式，更能亲手构建和调试这个强大的工具。

1. 为什么需要PGD：当优化遇上约束

想象你正在设计一个推荐系统，需要确保权重参数都落在0到1之间；或者训练一个物理模型，某些参数必须满足能量守恒定律。这些场景都提出了一个共同需求：在特定约束条件下进行优化。

普通梯度下降对此束手无策——它可能让参数跑到约束范围之外。PGD的巧妙之处在于每次参数更新后，通过 投影操作 将参数拉回可行域。这种"走一步，拉一步"的策略，既保持了梯度下降的简洁性，又满足了约束条件。

典型应用场景包括：

• 非负矩阵分解（NMF）中的非负约束
• 投资组合优化中的权重总和为1
• 物理模拟中的参数合法性保持

# 简单的投影函数示例：将参数限制在[0,1]区间
def project_to_unit_interval(x):
    return np.clip(x, 0, 1)

2. PGD核心原理拆解：三步理解投影梯度下降

PGD算法的每个迭代可以分解为三个关键步骤：

1. 梯度计算 ：与传统梯度下降相同，计算当前参数处的梯度
2. 参数更新 ：沿负梯度方向移动一定步长
3. 投影操作 ：将更新后的参数映射到约束集合中

数学表达上，第t次迭代的更新公式为：

x_{t+1} = Π_C(x_t - α∇f(x_t))

其中Π_C表示向约束集合C的投影操作。

注意：投影操作的定义依赖于约束集合的形状。对于简单的区间约束，投影就是截断；对于更复杂的约束（如球约束或仿射约束），投影可能需要求解一个子优化问题。

3. 实战准备：构建带约束的线性回归问题

让我们设计一个具体的优化问题：假设我们有一组房屋面积与价格的数据，需要拟合一个线性模型，但要求斜率参数在[-1,1]之间，截距项非负。

问题定义 ：

minimize ½||Xw - y||²
subject to w[0] ∈ [-1,1], w[1] ≥ 0

首先准备模拟数据：

import numpy as np
import matplotlib.pyplot as plt

# 生成模拟数据
np.random.seed(42)
X = 2 * np.random.rand(100, 1)  # 面积特征
y = 4 + 3 * X + np.random.randn(100, 1)  # 真实价格（带噪声）

# 添加偏置项
X_b = np.c_[np.ones((100, 1)), X]

4. 从零实现PGD算法

现在让我们完整实现PGD算法，包含梯度计算、参数更新和投影三个核心组件：

def pgd_linear_regression(X, y, constraints, alpha=0.1, max_iter=1000, tol=1e-6):
    """
    带约束的线性回归PGD实现
    
    参数：
    X: 特征矩阵 (m samples x n features)
    y: 目标值 (m x 1)
    constraints: 约束条件列表 [(min,max) for each param]
    alpha: 学习率
    max_iter: 最大迭代次数
    tol: 收敛阈值
    
    返回：
    w: 优化后的参数
    losses: 损失历史记录
    """
    m, n = X.shape
    w = np.random.randn(n, 1)  # 随机初始化参数
    
    losses = []
    for i in range(max_iter):
        # 计算梯度
        grad = X.T.dot(X.dot(w) - y) / m
        
        # 参数更新
        w_new = w - alpha * grad
        
        # 投影操作
        for j in range(n):
            low, high = constraints[j]
            w_new[j] = np.clip(w_new[j], low, high)
        
        # 检查收敛
        loss = 0.5 * np.mean((X.dot(w) - y)**2)
        losses.append(loss)
        
        if np.linalg.norm(w_new - w) < tol:
            break
            
        w = w_new
    
    return w, losses

# 定义约束：w0 ∈ [-1,1], w1 ≥ 0
constraints = [(-1, 1), (0, None)]

# 运行PGD
w_pgd, losses = pgd_linear_regression(X_b, y, constraints, alpha=0.1, max_iter=1000)

5. 算法可视化与调试技巧

理解优化算法最好的方式就是观察它的运行过程。让我们可视化损失下降和参数轨迹：

# 损失函数下降曲线
plt.figure(figsize=(12, 4))
plt.subplot(121)
plt.plot(losses)
plt.xlabel('Iteration')
plt.ylabel('Loss')
plt.title('Loss Convergence')

# 参数空间轨迹
w0_vals = np.linspace(-2, 2, 100)
w1_vals = np.linspace(-1, 5, 100)
W0, W1 = np.meshgrid(w0_vals, w1_vals)
loss_surface = np.array([0.5 * np.mean((X_b.dot([w0, w1]) - y)**2) 
                         for w0, w1 in zip(W0.ravel(), W1.ravel())]).reshape(W0.shape)

plt.subplot(122)
plt.contourf(W0, W1, np.log(loss_surface + 1), levels=20, alpha=0.8)
plt.colorbar()
plt.plot(w_pgd[0], w_pgd[1], 'ro', markersize=10)
plt.xlim([-2, 2])
plt.ylim([-1, 5])
plt.xlabel('w0 (slope)')
plt.ylabel('w1 (intercept)')
plt.title('Parameter Space with Constraint')
plt.axvline(-1, color='k', linestyle='--')  # 约束边界
plt.axvline(1, color='k', linestyle='--')
plt.axhline(0, color='k', linestyle='--')
plt.tight_layout()
plt.show()

调试PGD的实用技巧：

• 学习率选择 ：从0.01开始尝试，观察损失是否稳定下降
• 投影验证 ：确保每次迭代后参数确实满足约束
• 梯度检查 ：用数值梯度验证解析梯度的正确性
• 约束敏感性 ：尝试不同约束条件，观察解的变化

6. PGD的变体与进阶应用

掌握了基础PGD后，我们可以探索它的几种实用变体：

1. 随机PGD (SPGD)
每次迭代随机选取一个样本计算梯度，适合大规模数据集：

def stochastic_pgd(X, y, constraints, alpha=0.01, n_epochs=50, batch_size=10):
    m, n = X.shape
    w = np.random.randn(n, 1)
    losses = []
    
    for epoch in range(n_epochs):
        indices = np.random.permutation(m)
        X_shuffled = X[indices]
        y_shuffled = y[indices]
        
        for i in range(0, m, batch_size):
            X_batch = X_shuffled[i:i+batch_size]
            y_batch = y_shuffled[i:i+batch_size]
            
            grad = X_batch.T.dot(X_batch.dot(w) - y_batch) / batch_size
            w = np.clip(w - alpha * grad, 
                        [c[0] for c in constraints], 
                        [c[1] for c in constraints])
            
        loss = 0.5 * np.mean((X.dot(w) - y)**2)
        losses.append(loss)
    
    return w, losses

2. 加速PGD
引入动量项加速收敛：

def accelerated_pgd(X, y, constraints, alpha=0.1, beta=0.9, max_iter=1000):
    m, n = X.shape
    w = np.random.randn(n, 1)
    v = np.zeros_like(w)
    losses = []
    
    for i in range(max_iter):
        grad = X.T.dot(X.dot(w) - y) / m
        v = beta * v + (1 - beta) * grad
        w = np.clip(w - alpha * v, 
                   [c[0] for c in constraints], 
                   [c[1] for c in constraints])
        
        loss = 0.5 * np.mean((X.dot(w) - y)**2)
        losses.append(loss)
    
    return w, losses

3. 自适应步长PGD
根据梯度变化自动调整步长：

def adaptive_pgd(X, y, constraints, alpha=0.1, rho=0.9, max_iter=1000):
    m, n = X.shape
    w = np.random.randn(n, 1)
    grad_sq = np.zeros_like(w)
    losses = []
    
    for i in range(max_iter):
        grad = X.T.dot(X.dot(w) - y) / m
        grad_sq = rho * grad_sq + (1 - rho) * grad**2
        adjusted_alpha = alpha / (np.sqrt(grad_sq) + 1e-8)
        w = np.clip(w - adjusted_alpha * grad, 
                   [c[0] for c in constraints], 
                   [c[1] for c in constraints])
        
        loss = 0.5 * np.mean((X.dot(w) - y)**2)
        losses.append(loss)
    
    return w, losses

7. 工程实践中的挑战与解决方案

在实际项目中应用PGD时，你可能会遇到以下几个典型挑战：

挑战1：投影计算复杂度高
对于复杂约束（如半正定约束），投影操作本身可能就是一个优化问题。解决方案：

• 寻找约束的特殊结构，设计高效投影算法
• 使用近似投影或代理约束
• 考虑交替方向乘子法（ADMM）等其他约束优化方法

挑战2：非凸问题的局部最优
PGD只能保证收敛到局部最优。应对策略：

• 多组随机初始化
• 结合模拟退火等随机搜索技术
• 使用更复杂的优化器如Adam的约束版本

挑战3：约束冲突或不可行
当约束过于严格时，可能没有可行解。诊断方法：

• 检查约束的相容性
• 逐步放松约束，观察解的变化
• 引入松弛变量处理硬约束

# 处理冲突约束的示例：松弛变量法
def pgd_with_slack(X, y, constraints, lambda_=0.1, alpha=0.1, max_iter=1000):
    m, n = X.shape
    w = np.random.randn(n, 1)
    s = np.zeros_like(w)  # 松弛变量
    losses = []
    
    for i in range(max_iter):
        grad_w = X.T.dot(X.dot(w) - y) / m + lambda_ * s
        grad_s = lambda_ * (w - s)
        
        w -= alpha * grad_w
        s += alpha * grad_s
        
        # 投影
        for j in range(n):
            low, high = constraints[j]
            if low is not None and w[j] < low:
                s[j] += w[j] - low
                w[j] = low
            if high is not None and w[j] > high:
                s[j] += w[j] - high
                w[j] = high
        
        loss = 0.5 * np.mean((X.dot(w) - y)**2) + 0.5 * lambda_ * np.sum(s**2)
        losses.append(loss)
    
    return w, losses

在真实项目中，我发现PGD对学习率的选择特别敏感。一个实用的技巧是先用较大的学习率快速接近解，然后逐步减小学习率进行精细调整。另一个常见陷阱是忘记检查投影后的参数是否真的满足了约束——特别是在处理复杂约束时，建议在代码中添加断言检查。