用Python手写DFP和BFGS算法：从公式推导到代码实现，搞定机器学习中的优化问题

用Python手写DFP和BFGS算法从公式推导到代码实现搞定机器学习中的优化问题优化算法是机器学习工程师工具箱中的核心组件。当面对逻辑回归、支持向量机等模型的参数优化时传统的梯度下降法往往显得力不从心——收敛速度慢、对学习率敏感、在高维空间中容易陷入局部最优。拟牛顿法Quasi-Newton Methods作为二阶优化算法的代表通过近似海森矩阵Hessian Matrix巧妙地平衡了计算效率与收敛速度。本文将深入探讨两种经典的拟牛顿法DFPDavidon-Fletcher-Powell和BFGSBroyden-Fletcher-Goldfarb-Shanno。不同于简单的API调用我们会从数学原理出发手把手推导更新公式最终实现可直接集成到机器学习项目中的Python代码。通过对比这两种算法在Rosenbrock函数等经典测试案例上的表现您将掌握如何根据实际问题特点选择合适的优化器。1. 拟牛顿法牛顿法的智能近似牛顿法之所以强大在于它利用了目标函数的二阶导数信息。其迭代公式为x_{k1} x_k - H_k^{-1} g_k其中H_k是海森矩阵g_k是梯度。然而计算和存储完整的海森矩阵及其逆矩阵对于高维问题几乎不可行——这正是拟牛顿法要解决的痛点。拟牛顿法的核心思想是用迭代更新的方式构造海森矩阵或其逆矩阵的近似同时满足以下关键条件B_{k1} s_k y_k 拟牛顿条件这里s_k x_{k1} - x_k表示参数更新步长y_k g_{k1} - g_k表示梯度变化量。这个条件确保了近似矩阵能够正确反映目标函数的局部曲率信息。1.1 为什么需要不同的更新公式DFP和BFGS虽然同属拟牛顿法家族但采取了不同的近似策略DFP直接构造海森逆矩阵D_k的近似BFGS先构造海森矩阵B_k的近似再通过数学变换得到其逆矩阵这种差异导致了它们在数值稳定性和收敛速度上的不同表现。实际应用中BFGS通常表现更优这也是为什么现代机器学习框架如scikit-learn默认采用L-BFGSBFGS的内存优化版本作为逻辑回归等模型的优化器。2. DFP算法数学推导与实现细节2.1 公式推导从拟牛顿条件到秩2修正DFP算法的目标是通过秩2修正Rank-2 Update迭代更新海森逆近似矩阵D_k。其更新公式为D_{k1} D_k (s_k s_k^T)/(s_k^T y_k) - (D_k y_k y_k^T D_k)/(y_k^T D_k y_k)这个看似复杂的公式实际上可以通过严格的数学推导得到。让我们分解关键步骤设定修正形式假设更新后的矩阵可以表示为当前矩阵加上两个秩1矩阵的组合D_{k1} D_k α u u^T β v v^T代入拟牛顿条件将上述形式代入D_{k1} y_k s_k得到D_k y_k α u (u^T y_k) β v (v^T y_k) s_k选择基向量令u s_kv D_k y_k可以解出α 1/(s_k^T y_k), β -1/(y_k^T D_k y_k)2.2 Python实现关键代码def dfp_update(D, s, y): DFP更新海森逆近似矩阵 sTy np.dot(s, y) if sTy 1e-10: # 避免数值不稳定 return np.eye(D.shape[0]) # 计算各项分子 ssT np.outer(s, s) Dy np.dot(D, y) DyDyT np.outer(Dy, Dy) yTDy np.dot(y, Dy) # 组合更新项 D_new D ssT / sTy - DyDyT / yTDy return D_new实现要点使用np.outer计算向量外积检查s_k^T y_k的正定性避免数值问题保持矩阵对称性数学上保证实现中需注意浮点精度3. BFGS算法更稳定的选择3.1 为什么BFGS优于DFPBFGS算法通过不同的数学路径达到了更稳定的数值表现直接维护正定性BFGS更新公式能保证近似矩阵的正定性只要初始矩阵正定且s_k^T y_k 0更好的曲率估计对非二次函数的局部近似更准确自校正特性对线搜索误差的鲁棒性更强3.2 BFGS更新公式推导BFGS的海森逆更新公式为D_{k1} (I - ρ_k s_k y_k^T) D_k (I - ρ_k y_k s_k^T) ρ_k s_k s_k^T其中ρ_k 1/(y_k^T s_k)。这个形式看起来复杂但可以通过Sherman-Morrison-Woodbury公式从海森矩阵的BFGS更新推导得到。3.3 Python实现def bfgs_update(D, s, y): BFGS更新海森逆近似矩阵 sTy np.dot(s, y) if sTy 1e-10: return np.eye(D.shape[0]) I np.eye(D.shape[0]) rho 1.0 / sTy # 计算各项中间结果 syT np.outer(s, y) ysT np.outer(y, s) term1 I - rho * syT term2 I - rho * ysT # 组合更新 D_new term1 D term2 rho * np.outer(s, s) return D_new性能优化技巧利用矩阵乘法结合律减少计算量预先计算公共项rho使用运算符进行矩阵乘法Python 3.54. 算法框架与线搜索策略4.1 通用拟牛顿法框架def quasi_newton(f, grad_f, x0, methodbfgs, tol1e-6, max_iter100): n len(x0) x np.array(x0, dtypenp.float64) D np.eye(n) # 初始化为单位矩阵 iter_num 0 while iter_num max_iter: g grad_f(x) if np.linalg.norm(g) tol: break # 计算搜索方向 d -np.dot(D, g) # Armijo线搜索 alpha armijo(f, grad_f, x, d) # 更新参数 x_new x alpha * d g_new grad_f(x_new) # 计算s和y s x_new - x y g_new - g # 更新海森逆近似 if method dfp: D dfp_update(D, s, y) elif method bfgs: D bfgs_update(D, s, y) x x_new iter_num 1 return x, f(x), iter_num4.2 Armijo线搜索实现def armijo(f, grad_f, x, d, alpha_init1.0, c11e-4, beta0.5, max_iter20): Armijo准则线搜索 alpha alpha_init fx f(x) grad grad_f(x) slope np.dot(grad, d) for _ in range(max_iter): x_new x alpha * d if f(x_new) fx c1 * alpha * slope: return alpha alpha * beta return alpha # 返回最后尝试的alpha值参数选择经验c1通常取1e-4到1e-2beta在0.1到0.8之间初始步长alpha_init可根据问题规模调整5. 实战测试从二次函数到Rosenbrock5.1 二次函数测试案例考虑简单的二次函数f(x) 0.5x₁² 2x₂²其最优解显然是(0,0)。我们设置初始点为(2,1)比较DFP和BFGS的表现# 定义目标函数和梯度 def f_quad(x): return 0.5 * x[0]**2 2 * x[1]**2 def grad_f_quad(x): return np.array([x[0], 4 * x[1]]) # 测试 x0 [2.0, 1.0] x_dfp, f_dfp, iter_dfp quasi_newton(f_quad, grad_f_quad, x0, methoddfp) x_bfgs, f_bfgs, iter_bfgs quasi_newton(f_quad, grad_f_quad, x0, methodbfgs)典型输出DFP结果 最优解[ 1.42e-14 -7.11e-15] 迭代次数3 BFGS结果 最优解[ 1.42e-14 -7.11e-15] 迭代次数3对于这个简单的凸二次问题两种算法表现相当。5.2 Rosenbrock函数挑战Rosenbrock函数是著名的优化测试函数f(x) 100(x₂ - x₁²)² (1 - x₁)²其最优解为(1,1)但存在一个狭窄弯曲的山谷使得优化算法容易振荡。我们从(-1.2,1.0)开始测试def f_rosen(x): return 100 * (x[1] - x[0]**2)**2 (1 - x[0])**2 def grad_f_rosen(x): dx1 -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0]) dx2 200 * (x[1] - x[0]**2) return np.array([dx1, dx2]) x0 [-1.2, 1.0] x_bfgs, f_bfgs, iter_bfgs quasi_newton(f_rosen, grad_f_rosen, x0, methodbfgs) x_dfp, f_dfp, iter_dfp quasi_newton(f_rosen, grad_f_rosen, x0, methoddfp)典型输出BFGS结果 最优解[0.99999999 0.99999998] 迭代次数18 DFP结果 最优解[0.99999875 0.9999975] 迭代次数25在这个更具挑战性的案例中BFGS展现了更优的性能——不仅收敛更快而且最终解更接近理论最优。6. 机器学习中的应用技巧将DFP/BFGS应用于逻辑回归等模型时需要注意以下实践要点特征缩放像所有基于梯度的优化器一样拟牛顿法受益于特征标准化正则化处理L2正则化可以改善海森矩阵的条件数批量处理对于大规模数据可考虑随机版本的拟牛顿法参数初始化合理的初始参数能显著减少迭代次数以下是一个逻辑回归的拟牛顿法实现框架class LogisticRegression: def __init__(self, methodbfgs): self.method method self.weights None def sigmoid(self, z): return 1 / (1 np.exp(-z)) def loss(self, X, y, w): z np.dot(X, w) h self.sigmoid(z) return -np.mean(y * np.log(h) (1-y) * np.log(1-h)) def gradient(self, X, y, w): z np.dot(X, w) h self.sigmoid(z) return np.dot(X.T, (h - y)) / len(y) def fit(self, X, y): # 添加偏置项 X np.c_[np.ones(len(X)), X] # 初始化权重 w0 np.zeros(X.shape[1]) # 拟牛顿优化 self.weights, _, _ quasi_newton( lambda w: self.loss(X, y, w), lambda w: self.gradient(X, y, w), w0, methodself.method ) def predict(self, X): X np.c_[np.ones(len(X)), X] return (self.sigmoid(np.dot(X, self.weights)) 0.5).astype(int)7. 高级话题与扩展方向对于希望深入优化算法的读者以下方向值得探索L-BFGS有限内存版本的BFGS适用于高维问题稀疏拟牛顿法针对稀疏梯度问题的特殊处理非凸优化处理神经网络等非凸问题的改进策略分布式实现如何将拟牛顿法扩展到分布式计算环境实际项目中成熟的优化库如SciPy的optimize模块通常已经实现了高度优化的拟牛顿算法。但理解底层原理对于调试和特殊问题处理至关重要——比如当标准算法收敛缓慢时您可能需要调整线搜索参数或考虑问题的重新参数化。