使用BFGS公式实现阻尼牛顿法, 并实现L-BFGS:有限内存BFGS

本文侧重于数学推导, 相信读者已经对Logistic Regression已经有了充分的认识.

文章主要参考目录如下:

Broyden类算法：BFGS算法的迭代公式推导

12-lect-QN.pdf (pku.edu.cn)

代码可见仓库

mr-wuliu/MECHINELEARNING

符号约定

$$\begin{aligned} &s^k = x^{k+1}-x^k\\ &y^k = \nabla f(x^{k+1}) - \nabla f(x^k)\\ &B^k = \nabla^2f(x^k)\\ &H^k = (B^k)^{-1} \end{aligned}$$

$\nabla$为梯度

基本框架:

拟牛顿法框架

$$\begin{aligned} &初始化x_{n\times 1},H^k_{n\times n}, k = 0 \\ &\text{while }\|\nabla f(x_k) \|> \varepsilon \text{ do}\\ &\ \ \ \ d^k = - H^k\nabla f(x_k)\\ &\ \ \ \ wolfe搜索产生步长a_k, x^{k+1} = x^k +a_kd^k\\ &\ \ \ \ \text{update}H^{k+1} \\ &\ \ \ \ k\leftarrow k+1\\ &\text{end while} \end{aligned}$$

计算梯度

$$\nabla\mathscr{J}(\theta)=\frac1m X^T(\bar y - y)$$

如果加正则化项, 损失函数的形式为:

$$\mathscr{J}(\theta) = -\sum^m_{i=1}(y^{(i)}\ln p_1 + (1-y^{(i)}\ln (1-p_1) ) +\frac{\lambda }{2m}\|\theta\|^2_2$$

其梯度为:

$$\nabla\mathscr{J}(\theta)=\frac1m X^T(\bar y - y) +\frac{\lambda}{m}\theta$$

计算Hessian矩阵

设$f(x)$二阶连续可微. $\nabla f(x)$在$x^{k+1}$处y一阶泰勒近似为:

$$\nabla f(x) = \nabla f(x^{k+1}) +\nabla^2f(x^{k+1})(x-x^{k+1{}})\\ +O(\|x-x^{k+1}\|_2^2)$$

如果我们忽略高阶项, 可以得到:

$$B^{k+1}s^k = y^k \tag{1}$$

或者

$$H^{k+1} y^k = s^k\tag{2}$$

公式(1)和公式(2)称为割线方程.

使用BFGS公式:

$$H^{k+1}=\left(I-\frac{s^k(y^k)^T}{(s^k)^Ty^k}\right)^TH^k \left( I-\frac{s^k(y^k)^T}{(s^k)^Ty^k} \right) +\frac{s^k(s^k)^T}{(s^k)^Ty^k)}$$

其核心思想是对$H^k$进行秩二更新.

问题在于, $H^{k+1}$是根据割线方程近似得出的, 迭代开始时它就不是一个准确的二阶导逆矩阵(初始值为单位矩阵). 只不过随着迭代进程, 会快速收敛到正确的Hessian矩阵. 割线方程如果是基于添加正则化项的公式得出的, 我认为$H^{k+1}$就不需要添加正则化项了.

秩二更新

对于拟牛顿矩阵$B^k\in\mathbb{R}^{n\times n}$, 设$0\neq u,v\in\mathbb{R}^{n}$, 且$a,b\in\mathbb{R}$待定. 则秩二更新形式为:

$$B^{k+1} = B^k + auu^T +bvv^T$$

代入割线方程可得:

$$(au^Ts^k)u + (bv^Ts^k)v = y^k - B^ks^k$$

取合适的$a,b,u,v$, 得到:

$$B^{k+1} = B^k +\frac{y^k(y^k)^T}{(s^k)^Ty^k} - \frac{B^ks^k(B^ks^k)^T}{(s^k)^TB^ks^k}$$

对其使用两次两次Sherman-Morrison公式,:

$$(A+\frac{uv^T}{t})^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{t+ v^TA^{-1}u}$$

经过繁琐的化简, 最终就可以得到基于$H^k$更新的BFGS公式.

推导也可以使用SMW公式, 但是我对该公式不太熟悉.

Wolfe准则搜索步长

拟牛顿矩阵正定的充分条件可以是:

1. $H^k$正定
2. 满足曲率条件$(s^k)^Ty^k>0, \forall k \in \mathbb{N}$

使用wolfe准则搜索步长, 可以满足曲率条件, 保证拟牛顿矩阵的正定性.

定义:

设$d^k$是点$x^k$的下降方向, 若:

$$f(x^k + \alpha d^k)\leq f(x^k) + c_1\alpha \nabla f(x^k)^Td^k\\ \nabla f(x^k+\alpha d^k)^Td^k \geq c_2\nabla f(x^k)^Td^k$$

称步长$\alpha$满足wolfe准则, 其中$c_1 < c_2$且$c_1,c_2 \in (0,1 )$

可以使用二分法快速得到合适的$\alpha$

L-BFGS方法

BFGS方法速度很快, 但是仍然有不足之处. 对于特征数量很多的Hessian矩阵, 其计算和存储都是一个难以解决的问题. 很多时候我们可以通过加内存的方式来解决, 但是单纯通过硬件的提升是有上限的. 我们需要设计更好的算法来解决问题.

我们注意到, $H^{k+1}$的计算都是基于一个初始$H^0$和中间过程$s^k,y^k$得出的. 而且我们最终需要的只是一个牛顿方向$d^k$, 并不需要一个完整的Hessian矩阵. 如果存在一个算法能够直接得出牛顿方向$d^k$, 就能够大大节省内存. 甚至进一步提高了计算速度(大型矩阵运算的开销很大).

事实上确实存在这样一个神奇的算法, 能够基于$s^k,y^k$等信息, 直接计算$d^k$ .而$H^{k}$也不再显式存储, 而是通过保存有限个$s^k,y^k$.

"计算"Hessian矩阵

方便起见, 将BFGS中计算$H^k$的矩阵改写为:

$$H^{k+1} = (V^k)^TH^kV^k + \rho_ks^k(s^k)^T$$

其中

$$V^k =I-\frac{s^k(y^k)^T}{(s^k)^Ty^k};\rho_k =\frac1{(s^k)^Ty^k)}$$

如果对$H^{k+1}$进行有限的m次的展开, 可以得到:

$$\begin{aligned} H^{k+1}&=(\prod_{i=k-m}^kV^i)^TH^{k-m}(\prod_{i=k-m}^kV^i)\\ &+\rho_{k-m}(\prod_{i=k-m+1}^kV^i)^Ts^{k-m}(s^{k-m})^T(\prod_{i=k-m+1}^kV^i)\\ &+\cdots+\\ &+ \rho_{k-1}s^{k-1}(s^{k-1})^T \end{aligned}$$

看上去展开后依然还要存储$H^{k-m}$, 不过存在一个神奇算法, 能让我们不需要存储$H^{k-m}$, 也不存储$H^{k}$, 直接计算得到$d^k$.

双循环递归算法, 它采用迭代式的结构, 尽可能地节省内存开销. 并且避免了矩阵乘法, 大大加快了运算速度.

如果等式两边分别乘上$\nabla f(\theta^k)$ 可以得到第一:

$$\begin{aligned} -d^k&=(\prod_{i=k-m}^kV^i)^TH^{k-m}(\prod_{i=k-m}^kV^i\nabla f(\theta^k)\\ &+\rho_{k-m}(\prod_{i=k-m+1}^kV^i)^Ts^{k-m}(s^{k-m})^T(\prod_{i=k-m+1}^kV^i)\nabla f(\theta^k)\\ &+\cdots+\\ &+ \rho_{k-1}s^{k-1}(s^{k-1})^T\nabla f(\theta^k) \end{aligned}$$

其中的$(\prod_{i=k-m+1}^kV^i)\nabla f(\theta^k)$过程可以存储起来, 用迭代的方式得到, 节省了计算量.这就是第一个循环.

如果记:

$$q =(\prod_{i=k-m}^{k-1}V^i)\nabla f(\theta^k)\\ \alpha_i=\rho_{k-j}s^{k-j}(\prod_{i=k-l+1}^{k-1}V^i)\nabla f(\theta^k)$$

那么有:

$$-d^k = (\prod_{i=k-m}^{k-1}V^i)H^{k-m}q + (\prod_{i=k-m+1}^{k-1}V^i)s^{k-m}\alpha_{k-m}+\cdots+s^{k-1}\alpha_{k-1}$$

另一个循换为了解决$d^k$方向的计算.

先计算$H^{k-m}q$. 展开$V^i$和后项合并可以得到

$$H^{k-m} - s^{k-m}\alpha_{k-m} + s^{k-m}\rho_{k-m}(y^{k})^TH^{k-m}\\ =H^{k-m} - s^{k-m}(\alpha_{k-m}-\rho_{k-m}(y^{k-m})^TH^0q)$$

依次递归就可以得到最终的结果$-d^k$

如此, L-BFGS算法只需要改动$H^k$的计算方式, 并保存有限$m$次迭代信息, 就能保证算法运行.

实现代码

将上述过程编写为代码即可. 实验得出, bfgs方法迭代次数远多于牛顿迭代法, 还是需要一定的次数的. 不过它利用了二阶矩阵的信息, 即使是近似的矩阵, 其收敛速度仍然是超线性的.

首先为L-BFGS和BFGS建立一个公共的可调用函数库.

import pandas as pd
import numpy as np

"""
Clementine Model
模块用于数据处理. 也部分作为公共库, 例如损失函数, 梯度等
Clementine 意为小橘子, 谐音局子, 表示写这个玩意儿像坐牢一样.
至于为什么不叫tangerine, 因为已经有这个库的名字了, 防止重名.
"""

'''
公共函数, 包含:
sigmoid() 激活函数
cost() 计算损失
nabla() 用解析的方式计算梯度
scrubbing() 数据加常数项, 并可选归一化
wolfe() 搜索合适的alpha值
predict() 结果预测
'''


def sigmoid(x):
    return 1.0 / (1 + np.exp(-x))


def cost(X, theta, Y, c=False, Lambda=0):
    """
    损失函数, 采用负对数似然函数作为衡量标准
    正则化项采用L2正则化, 默认正则化项为零
    :param X: 样本集
    :param theta: 待定系数
    :param Y: 标签集
    :param c: 函数系数
    :param Lambda: 拟合程度
    :return:
    """
    #h = sigmoid(X @ theta)
    if c:
        m, n = X.shape
        a = np.multiply(X @ theta, Y)
        b = np.log(1 + np.exp(X @ theta))
        loss = (1.0 / m) * np.sum(-a + b) \
               + Lambda * theta.T @ theta / 2
    else:
        a = np.multiply(X @  theta, Y)
        b = np.log(1 + np.exp(X @ theta))
        loss = np.sum( -a  + b )\
                +Lambda * theta.T @ theta / 2
    return loss


def nabla(X, theta, Y, Lambda=0, c=False):
    """
    实现nabla算子
    :param X:
    :param theta:
    :param Y:
    :param Lambda: 正则化程度
    :param c: 损失函数系数
    :return:
    """
    grad = X.T @ (sigmoid(X @ theta) - Y) + Lambda * theta
    if c:
        m, n = X.shape
        return (1.0 / m) * grad
    return grad


def scrubbing(data, stander=True):
    """
    对数据做归一化处理并且在第一列添加常数项
    :param data:
    :param stander:
    :return:
    """
    row, columns = data.shape
    ones = np.ones((row, 1))
    if stander:
        x = data
        data = (x - np.min(x, axis=0)) / (np.max(x, axis=0) - np.min(x, axis=0))
    data = np.append(ones, data, axis=1)
    return data


def wolfe(x, theta, y, d, c, Lambda, log=False):
    c1 = 0.1
    c2 = 0.9
    alpha = 25
    a = 0
    b = 50
    max_epoch = 40
    epoch = 0
    while a < b and epoch < max_epoch:
        if cost(x, theta + alpha * d, y, c=c, Lambda=Lambda) <= \
                cost(x, theta, y, c=c, Lambda=Lambda) + \
                c1 * alpha * nabla(x, theta, y, c=c, Lambda=Lambda).T @ d:
            if nabla(x, theta + alpha * d, y, c=c, Lambda=Lambda).T @ d >= \
                    c2 * nabla(x, theta, y, c=c, Lambda=Lambda).T @ d:
                break
            else:
                a = alpha
        else:
            b = alpha
        alpha = (a + b) / 2
        epoch += 1
    if log and epoch == max_epoch:
        alpha = 0.05
        print("wolfe搜索失败")
    return alpha

def predict(X, theta, y):
    m, n = X.shape
    pre = np.zeros((m, 1))

    for i in range(m):
        activate = sigmoid(np.dot(X[i, :], theta))
        if activate > 0.5:
            pre[i] = 1
    count = 0
    for i in range(m):
        if pre[i] == y[i]:
            count += 1
    return count / m

BFGS算法实现

def bfgs_method(x_train, y_label, max_iter=5000,
                epsilon=1e-4, Lambda=0,
                c=False, log=False):
    """
    实现阻尼牛顿法
    :param x_train:
    :param y_label:
    :param max_iter:
    :param epsilon: 精度
    :param Lambda: 正则化项
    :param c: 损失函数系数
    :param log: 日志选项
    :return:
    """

    m, n = x_train.shape
    theta_k = np.ones((n, 1))
    I = np.eye(n)
    Hk = I
    grad_k = ct.nabla(x_train, theta_k, y_label,
                      c=c, Lambda=Lambda)
    k = 0
    while np.linalg.norm(grad_k) > epsilon and k < max_iter:
        # 牛顿方向
        d = - Hk @ grad_k
        # 搜索步长
        alpha_k = ct.wolfe(x_train, theta_k, y_label, d,
                        c, Lambda)
        # 迭代
        theta_k1 = theta_k + alpha_k * d
        # 更新Hessian
        s_k = theta_k1 - theta_k
        grad_k1 = ct.nabla(x_train, theta_k1, y_label,
                           c=c, Lambda=Lambda)
        y_k = grad_k1 - grad_k
        if s_k.T @ y_k <= 0:
            print("矩阵非正定,求解失败, 退出")
            print(s_k.T @ y_k)
            break
        # BFGS公式
        Hk = ((I - (s_k @ y_k.T) / (s_k.T @ y_k))).T @ Hk @ ((I - (s_k @ y_k.T) / (s_k.T @ y_k))) + ((s_k @ s_k.T) / (s_k.T @ y_k))
        grad_k = grad_k1
        theta_k = theta_k1
        k += 1
        if log and  k % 200 == 0:
            print(k)
    if log and max_iter == k:
        print("迭代超范围,|grad|=",np.linalg.norm(grad_k))
    if log:
        print("损失:", ct.cost(x_train, theta_k, y_label, c=c, Lambda=Lambda))
    return theta_k

L-BFGS算法实现

import numpy as np
import clementine as ct
import time


def two_loop(s, y, rho, gk,col):
    """
    通过序列s , y rho 序列计算出d_k方向
    :param gk: 当然梯度
    :param s: 字典, 点差集合
    :param y: 字典, 梯度差集合
    :param rho: 字典, 标量集合
    :return: d_k, 牛顿方向
    """
    n = len(s)  # 计算列表长度

    # 选择近似矩阵

    a = np.ones((n, 1))  # a[i] 同样为标量
    q = gk.copy()  # n * 1 矩阵

    # 第一次循环
    for i in range(n - 1, -1, -1):
        a[i] = rho[i] @ s[i].T @ q
        q -= a[i] * y[i]

    # 选择近似矩阵
    h0 = np.eye(col)
    if n >= 1:
        h0 = (s[-1].T @ y[-1]) / (y[-1].T @ y[-1]) * h0
    # 第二次循环
    d = h0 @ q

    for i in range(n):
        b = rho[i] @ y[i].T @ d
        d += s[i] @ (a[i] - b)
    return d


def l_bfgs_method(x_train, y_label, max_iter=5000,
                  epsilon=1e-4, Lambda=0, m=5,
                  c=False, log=False):
    # 参数初始化
    row, col = x_train.shape
    theta_k = np.zeros((col, 1))  # 初始位置
    grad_k = ct.nabla(x_train, theta_k, y_label,
                      Lambda=Lambda, c=c)
    k = 0
    s, y, rho = [], [], []

    while np.linalg.norm(grad_k) > epsilon and k < max_iter:
        d_k = -two_loop(s, y, rho, grad_k,col)

        alpha = ct.wolfe(x_train, theta_k, y_label, d_k,
                      c=c, Lambda=Lambda)
        theta_k1 = theta_k + alpha * d_k

        s_k = theta_k1 - theta_k
        grad_k1 = ct.nabla(x_train, theta_k1, y_label,
                           Lambda=Lambda, c=c)
        y_k = grad_k1 - grad_k

        if s_k.T @ y_k <= 0:
            print("矩阵非正定,求解失败, 退出")
            print(s_k.T @ y_k)
            break
        rho.append(1 / (s_k.T @ y_k))
        s.append(s_k)
        y.append(y_k)
        if (len(rho) > m):
            rho.pop(0)
            s.pop(0)
            y.pop(0)
        # 更新
        k += 1
        if log and k % 10 == 0:
            print("迭代次数",k,"||grad||",np.linalg.norm(grad_k))
            print("loss", ct.cost(x_train, theta_k1, y_label,
                                  Lambda=Lambda, c=c))
        grad_k = grad_k1
        theta_k = theta_k1
    # end while

    return theta_k