lbfgs.py

import numpy as np
from strong_wolfe import strong_wolfe
floatX = np.longdouble
def lbfgs(fhandle, x0, m=10, max_iter=200000, ftol=1e-8, step_tol=1e-8,gtol=1e-8):
    x_k = np.asarray(x0, dtype=floatX).ravel().copy()
    g_k = fhandle(x_k, 'grad')#梯度
    f_k = fhandle(x_k, 'f')
    p = -g_k.copy()#初始前进方向
    S   = []
    Y   = []
    Rho = []
    for k in range(1, max_iter + 1):
        mu = strong_wolfe(fhandle, x_k, p)#强wolfe条件线搜索
        x_kp1 = x_k + mu * p#位置更新
        # ---------- 本次位置更新完成，计算下一次位置更新的前进方向 ----------
        g_kp1 = fhandle(x_kp1, 'grad')
        f_kp1 = fhandle(x_kp1, 'f')
        y_k = g_kp1 - g_k
        s_k = x_kp1 - x_k
        # ---------- 收敛判断 ----------
        if (abs(f_kp1 - f_k) < ftol and
                np.linalg.norm(s_k) < step_tol and
                np.linalg.norm(g_kp1)<gtol):
            return x_kp1, f_kp1, k
        yTs = np.dot(y_k, s_k)
        if yTs < 1e-20:#曲率条件修正
            p = -g_kp1.copy()  # 本步使用最速下降
            x_k = x_kp1.copy()
            g_k = g_kp1.copy()
            f_k = f_kp1.copy()
            continue
        rho_k = 1.0 / np.dot(y_k, s_k)
        if len(S) == m+1:# 最多保留m+1步，丢弃最老的一对
            S.pop(0)
            Y.pop(0)
            Rho.pop(0)
        S.append(s_k.copy())
        Y.append(y_k.copy())
        Rho.append(rho_k.copy())
        q = g_kp1.copy()
        alphas = []
        # 第一循环：从最新到最旧（i = k, k-1, ..., k-m）
        for j in range(len(S)-1, -1, -1):
            alpha = Rho[j] * np.dot(S[j], q)
            q -= alpha * Y[j]
            alphas.append(alpha)
        if len(Y) >= 2:
            gamma = np.dot(Y[-2],S[-2]) / np.dot(Y[-2],Y[-2])
        else:
            gamma = np.dot(Y[-1],S[-1]) / np.dot(Y[-1],Y[-1])
        r = gamma * q
        # 第二循环：从最旧到最新（i = k-m, ..., k）
        for j in range(len(S)):
            beta = Rho[j] * np.dot(Y[j],r)
            alpha_j = alphas.pop()
            r += S[j] * (alpha_j - beta)
        p = -r
        g_k = g_kp1.copy()
        f_k = f_kp1.copy()
        x_k = x_kp1.copy()
    return x_k, f_k, max_iter