数值优化 —— 拟牛顿法 L-BFGS

主要来源和参考：

https://felix.blog.csdn.net/article/details/46389869

L-BFGS算法中的方向的计算方法: （这个伪代码形式的算法描述并不是很OK，不建议使用，建议跳过）

关键公式：

\(s_k = x_{k+1} - x_k\)

\(y_k = grad_{k+1} - grad_k\)

DeepSeek生成：

L-BFGS算法介绍与伪代码

思考过程：

好的，我现在需要帮助用户深入理解L-BFGS算法的具体细节。之前已经给出了介绍和伪代码，但用户现在要求“具体的描述”，可能需要更详细的数学推导、步骤解释或者实际应用中的注意事项。首先，我需要回顾之前的对话，确保不重复信息，同时补充用户可能需要的更深层次内容。

用户可能希望了解L-BFGS与标准BFGS的区别，特别是内存限制如何影响Hessian近似的构建。此外，双循环递归的具体实现可能还需要更详细的解释，包括每一步的数学依据。另外，初始Hessian的缩放因子gamma的选择也是一个关键点，需要更深入的解释。

还需要考虑用户的应用场景，比如大规模机器学习中的参数优化，这里可能需要讨论存储向量对的选择策略，如何平衡内存和收敛速度。此外，线搜索的具体实现方法，如回溯线搜索或满足Wolfe条件的方法，用户可能希望了解如何与L-BFGS结合使用。

另外，用户可能对L-BFGS的收敛性分析感兴趣，比如在凸和非凸问题下的表现，以及超线性收敛的条件。不过，这部分可能比较理论，是否需要深入取决于用户的需求层次。此外，伪代码中的参数m的选择建议，以及实际编程实现时的注意事项，比如数值稳定性问题，比如如何处理y_k和s_k的内积避免除以零，这些都需要提及。

需要确认是否有遗漏的关键步骤，比如如何处理初始迭代中向量对不足的情况，或者在存储队列已满时如何替换旧的向量对。此外，可能还需要解释为什么双循环递归能有效计算Hessian逆的近似，以及其背后的数学原理，比如通过递归应用BFGS的更新公式，而无需显式存储矩阵。

最后，可能需要给出一个具体的例子，比如在逻辑回归中的应用，或者对比其他优化算法如Adam的效果，但用户可能更关注算法本身的具体描述，所以例子可能不是必须的。重点还是放在算法步骤的详细分解和数学推导上，确保用户能够理解每一步的动机和实现方式。

（由于技术原因，联网搜索暂不可用）

L-BFGS 具体描述

L-BFGS（Limited-memory BFGS）是一种高效的拟牛顿优化算法，专为高维参数空间设计。以下是其核心组件和数学细节的完整解析：

L-BFGS伪代码详解（Python风格）

def L_BFGS(x0, grad_f, m=10, max_iter=100):
    x = x0.copy()
    g = grad_f(x)
    S, Y = [], []  # 存储s和y的队列
    
    for _ in range(max_iter):
        # 双循环计算搜索方向d = -H_k g
        q = g.copy()
        alphas = []
        for i in reversed(range(len(S))):
            rho_i = 1.0 / (Y[i].T @ S[i])  # 避免重复计算
            alpha_i = rho_i * (S[i].T @ q)
            alphas.append(alpha_i)
            q -= alpha_i * Y[i]
        
        # 计算初始Hessian尺度gamma
        if len(S) == 0:
            gamma = 1.0
        else:
            s_old, y_old = S[-1], Y[-1]
            gamma = (s_old.T @ y_old) / (y_old.T @ y_old)  # 确保正定性
        
        r = gamma * q
        
        for i in range(len(S)):
            beta_i = (Y[i].T @ r) * rho_i  # rho_i已预计算
            r += S[i] * (alphas[-(i+1)] - beta_i)  # 反向索引
        
        d = -r
        
        # 线搜索（需满足Wolfe条件）
        alpha = line_search_wolfe(f, grad_f, x, d)
        x_new = x + alpha * d
        g_new = grad_f(x_new)
        
        # 更新存储队列
        s_k = x_new - x
        y_k = g_new - g
        if y_k.T @ s_k > 1e-10:  # 避免数值问题
            S.append(s_k)
            Y.append(y_k)
            if len(S) > m:
                S.pop(0)
                Y.pop(0)
        
        x, g = x_new, g_new
    
    return x

注意！！！

由于上面的DeepSeek生成的伪代码实在是有些不好理解（可能还是有一些错误存在），因此我们又使用DeepSeek生成了一个可以运行的LBFG的代码（有一丢丢的小错误，人工修改后可以成功运行），给出可以运行的LBFG代码：

import numpy as np
from scipy.optimize import line_search

def rosenbrock(x):
    """Rosenbrock测试函数"""
    return 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2

def rosenbrock_grad(x):
    """Rosenbrock梯度"""
    return np.array([
        -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0]),
        200 * (x[1] - x[0]**2)
    ])

def lbfgs(f, grad_f, x0, m=10, max_iter=100, tol=1e-6):
    """
    L-BFGS优化算法实现
    
    参数:
        f - 目标函数
        grad_f - 梯度函数
        x0 - 初始点
        m - 存储的向量对数
        max_iter - 最大迭代次数
        tol - 梯度收敛阈值
        
    返回:
        x - 优化后的参数
    """
    x = x0.copy().astype(float)
    g = grad_f(x)
    s_list = []  # 存储s向量 (x的变化量)
    y_list = []  # 存储y向量 (梯度的变化量)
    
    for k in range(max_iter):
        # ------------------- 计算搜索方向 d = -H_k g -------------------
        q = g.copy()
        alpha_list = []
        
        # 反向循环 (从最新到最旧)
        for i in reversed(range(len(s_list))):
            s = s_list[i]
            y = y_list[i]
            rho = 1.0 / (y.T @ s)
            alpha = rho * (s.T @ q)
            alpha_list.append(alpha)
            q -= alpha * y
        
        # 计算初始Hessian尺度因子gamma
        if len(s_list) == 0:
            gamma = 1.0
        else:
            s_prev, y_prev = s_list[-1], y_list[-1]
            gamma = (s_prev.T @ y_prev) / (y_prev.T @ y_prev)
        
        r = gamma * q
        
        # 正向循环 (从最旧到最新)
        for i in range(len(s_list)):
            s = s_list[i]
            y = y_list[i]
            rho = 1.0 / (y.T @ s)
            beta = rho * (y.T @ r)
            r += s * (alpha_list[-(i+1)] - beta)  # 注意反向索引
        
        d = -r
        
        # ------------------- 线搜索 -------------------
        # 使用SciPy的强Wolfe线搜索
        alpha, _, _, _, _, _ = line_search(f, grad_f, x, d, gfk=g, maxiter=100)
        if alpha is None:  # 线搜索失败时使用保守步长
            alpha = 1e-4
        
        # ------------------- 更新参数 -------------------
        print(x, d, alpha)
        x_new = x + alpha * d
        g_new = grad_f(x_new)
        s = x_new - x
        y = g_new - g
        
        # ------------------- 存储管理 -------------------
        # 只保留正定条件 y.T@s > 0 的向量对
        if y.T @ s > 1e-10:
            s_list.append(s)
            y_list.append(y)
            if len(s_list) > m:
                s_list.pop(0)
                y_list.pop(0)
        
        # 准备下一次迭代
        x = x_new
        g = g_new
        
        # 检查收敛条件
        if np.linalg.norm(g) < tol:
            print(f"在 {k+1} 次迭代后收敛")
            break
    
    return x

# 测试优化器 -------------------------------------------------
if __name__ == "__main__":
    # 初始点 (经典测试点 [-1.2, 1.0])
    x0 = np.array([-1.2, 1.0])
    
    # 运行L-BFGS
    x_opt = lbfgs(rosenbrock, rosenbrock_grad, x0, m=5, max_iter=100)
    
    print("\n优化结果:")
    print(f"最小值点: {x_opt}")
    print(f"函数值: {rosenbrock(x_opt):.2e}")
    print(f"真实最小值位置: [1.0, 1.0]")

运行效果：