拟牛顿法BFGS算法与实现代码

import numpy as np
from scipy.optimize import minimize  # 用于对比验证

def bfgs_optimize(f, grad_f, x0, max_iter=1000, tol=1e-6):
    """
    BFGS 优化算法实现
    
    参数:
        f: 目标函数
        grad_f: 目标函数的梯度
        x0: 初始点
        max_iter: 最大迭代次数
        tol: 收敛容忍度（当梯度范数小于该值时停止）
    
    返回:
        x: 最优解
        history: 迭代历史记录（可选）
    """
    n = len(x0) 
    x = x0.copy()  # 浅拷贝  
    Hk = np.eye(n)  # 初始 Hessian 逆(即（B^(-1)_{k+1}）的近似 （即公式 1.5 中 的 Hk 取单位矩阵）; 
    history = [x.copy()]  # 记录迭代历史 第一个值; 
    
    for k in range(max_iter):
    
        gk = grad_f(x)  # 计算梯度
        
        if np.linalg.norm(gk) < tol:  # 检查是否收敛
            break
        
        dk = -np.dot(Bk, gk)  # 计算搜索方向  
        
        # 线搜索（这里采用简单的 Armijo 线搜索） 
        alpha = 1.0
        
        c = 0.25  # Armijo 条件参数   
        
        rho = 0.5  # 步长缩减因子  
        
        while f(x + alpha * dk) > f(x) + c * alpha * np.dot(gk, dk):
            alpha *= rho
        
        pk = alpha * dk  # 计算步长   n * 1 ; 
        x_new = x + sk   # sk 为 迭代变量差   dk | n* 1 
        
        # 计算新梯度
        gk_new = grad_f(x_new)
        qk = gk_new - gk   # qk 为 梯度差   q 
        
        # 更新 Hessian 逆的近似 H_{k+1} 公式 1.5 
        
        Hk_qk = np.dot(Hk, qk)           #  n * n , n * 1 ;  对于二维矩阵 计算矩阵积;
        qk_T_Hk_qk = np.dot(qk, Hk_qk)   #  对于一维矩阵计算内积   1 * n , n * 1; 
        pk_pkT = np.outer(pk, pk)        #  计算外积，outer里的矩阵都将展开为一维矩阵 得到 二维的n*n 矩阵;
        # Hk_pk_pkT_Hk = np.outer(Hk_qk, Hk_qk) / qk_T_Hk_qk  # 
        
        # BFGS 更新公式
        Hk = Hk + (1 + qk_T_Hk_qk / np.dot(qk, pk)) * (pk_pkT / np.dot(qk, pk)) - ( np.outer(pk, Hk_qk)+ np.outer(Hk_qk, pk)) / np.dot(qk, pk)
        # 根据公式1.5 
        x = x_new 
        
        history.append(x.copy())
    
    return x, history

# ===== 测试案例：Rosenbrock 函数优化 =====
def rosenbrock(x):
    """Rosenbrock 函数（经典优化测试函数）"""
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

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

# ===== 测试案例：Rosenbrock 函数优化 =====
def rosenbrock(x):
    """Rosenbrock 函数（经典优化测试函数）"""
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

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



# 初始点
x0 = np.array([-3.5, 3.5])

# 手动 BFGS 优化
x_opt, history = bfgs_optimize(rosenbrock, rosenbrock_grad, x0)
print("BFGS 手动优化结果:", x_opt)

# 用 SciPy 的 BFGS 验证。
res = minimize(rosenbrock, x0, method='BFGS', jac=rosenbrock_grad)
print("SciPy BFGS 优化结果:", res.x)

x1= np.array([-1,1])

结果：
BFGS 手动优化结果: [1. 1.]
SciPy BFGS 优化结果: [0.99999998 0.99999995]