《機器學習_06_優化_擬牛頓法實現(DFP,BFGS)》

一.簡介

經過前面幾節的介紹，你們能夠直觀的感覺到：對於大部分機器學習模型，咱們一般會將其轉化爲一個優化問題，因爲模型一般較爲複雜，難以直接計算其解析解，咱們會採用迭代式的優化手段，用數學語言描述以下：

\[\min_{v^k} f(x^k+v^k) \]

這裏目標函數爲\(f(x)\)，當前優化變量爲\(v^k\)，目標便是找到一個\(v^k\)對當前的\(x^k\)進行更新，使得函數值儘量的下降，若是目標函數一階可微，對其做一階泰勒展開便可獲得以下梯度降低的更新公式：

二.梯度降低

對目標函數做一階泰勒展開：

\[f(x^k+v^k)=f(x^k)+\triangledown f(x^k)^Tv^k \]

因此要使得\(f(x^k+v^k)<f(x^k)\)只須要使\(\triangledown f(x^k)^Tv^k<0\)便可，而若是取：

\[v^k=-\lambda_k \triangledown f(x^k),\lambda_k>0 \]

則必定能使\(\triangledown f(x^k)^Tv^k<0\)，因此，咱們就獲得了梯度降低的更新公式：

\[x^{k+1}=x^k+v^k=x^k-\lambda_k \triangledown f(x^k) \]

這裏\(\lambda_k\)通常能夠設置一個較小的定值，或者初始設置一個定值並使其隨着迭代次數增長而遞減；另外更好的作法是利用一維搜索：\(min_{\lambda_k}f(x^k-\lambda_k\triangledown f(x^k))\)求一個最優的\(\lambda_k\)，接下來咱們想一下下面兩個問題：

（1）梯度降低法必定能使得函數值降低嗎？

（2）若它能使函數值降低，則它是最優的方向嗎？

對於第一個問題，泰勒展開實際上是有一個條件的，那就是\(v^k\rightarrow 0\)，再結合上面的更新公式，若是\(\lambda_k\)取得過大時，咱們是不能省略泰勒展開後面的項的，並且後面項的取值也不必定能保證小於0，因此有時咱們設置的學習率\(\lambda_k\)較大時，函數值反而會上升。

因此，當\(v^k\)的取值大到不能忽略後面的項時，泰勒展開的二階項取值就必需要考慮其中了，因此這時梯度降低法未必時最優的方向，接下來咱們看下二階展開的狀況：

三.牛頓法

對其做二階泰勒展開：

\[f(x^k+v^k)=f(x^k)+\triangledown f(x^k)^Tv^k+\frac{1}{2}{v^k}^T\triangledown^2 f(x^k)v^k\\ =f(x^k)+g_k^Tv^k+\frac{1}{2}{v^k}^TH_kv^k \]

這裏爲了方便，記\(g_k=g(x^k)=\triangledown f(x^k),H_k=H(x^k)=\triangledown^2 f(x^k)\)，\(H_k\)表示Hessian矩陣在\(x^k\)處的取值，Hessian矩陣的定義：

\[H(x)=[\frac{\partial f}{\partial x_i\partial x_j}]_{n\times n} \]

對於大部分機器學習模型，一般目標函數是凸的，因此\(H(x)\)半正定，即對於\(\forall v^k\)，都有\({v^k}^TH_kv^k\geq0\)，此時，\(f(x^k)+g_k^Tv^k+\frac{1}{2}{v^k}^TH_kv^k\)是關於\(v^k\)的凸二次規劃問題，因此最優的\(v^k\)在其梯度爲0處取得：

\[\frac{\partial f(x^k+v^k)}{\partial v^k}=g^k+H^kv^k=0\\ \Rightarrow v^k=H_k^{-1}(-g_k) \]

能夠發現牛頓法對比梯度降低法，其實牛頓法是對梯度法的方向進行了一個改變\(H_k^{-1}\)，因此，咱們能夠獲得牛頓法的更新公式：

\[x^{k+1}=x^k-\lambda_kH_k^{-1}g_k\\ =x^k+\lambda_k p_k \]

這裏記\(p_k=-H_k^-1g_k\)；

能夠發現牛頓法的複雜有點高，由於要求解\(H_k^{-1}\)，那麼有沒有方便一點的方法呢？好比構建一個矩陣去近似\(H_k\)或者\(H_k^{-1}\)，這即是擬牛頓法的基本思想

四.擬牛頓條件

上面說到了利用一個矩陣去近似Hessian矩陣或者Hessian矩陣的逆，那麼這個近似矩陣須要知足怎樣的條件呢？咱們仍是從二階泰勒展開出發，稍微變換一下：

\[f(x^{k+1})=f(x^k)+g_k^T(x^{k+1}-x^k)+\frac{1}{2}(x^{k+1}-x^k)^TH_k(x^{k+1}-x^k) \]

兩邊對\(x^{k+1}\)求偏導可得：

\[g_{k+1}=g_k+H_k(x^{k+1}-x^k) \]

這即是擬牛頓條件，爲了方便，記\(y_k=g_{k+1}-g_k,\delta_k=x^{k+1}-x^k\)，因此：

\[y_k=H_k\delta_k \]

因此，擬牛頓法也要知足和\(H_k\)同樣的性質：

（1）正定性；

（2）知足擬牛頓條件

接下來，簡單證實一下若是知足性質（1）：正定性，更新時能夠知足函數值降低，假設\(G_k\)去近似\(H_k^{-1}\)，因此：\(G_k\succ 0\)，那麼迭代公式爲：

\[x^{k+1}=x^k-\lambda_kG_kg_k,\lambda_k>0 \]

將其帶入二階泰勒展開式中：

\[f(x^{k+1})=f(x^k)-\lambda_kg_k^TG_kg_k+\frac{1}{2}\lambda_k^2g_k^TG_k^TH_kG_kg_k \]

一般\(\lambda_k^2<<\lambda_k\)，因此能夠省略第三項，而第二項因爲\(G_k\succ 0\)，因此\(-\lambda_kg_k^TG_kg_k< 0\)，因此\(f(x^{k+1})<f(x^k)\)

五.DFP算法

DFP算法即是利用\(G_k\)近似\(H_k^{-1}\)的一種算法，它的構造很tricky，它假設每一步迭代中矩陣\(G_{k+1}\)是由\(G_k\)加上兩個附加項構成的，即：

\[G_{k+1}=G_k+P_k+Q_k \]

這裏\(P_k,Q_k\)是待定項，因爲須要知足擬牛頓條件，因此：

\[G_{k+1}y_k=\delta_k=G_ky_k+P_ky_k+Q_ky_k \]

這裏作一個tricky的假設：

\[P_ky_k=\delta_k\\ Q_ky_k=-G_ky_k \]

這樣的\(P_k,Q_k\)不難找到：

\[P_k=\frac{\delta_k\delta_k^T}{\delta_k^Ty_k}\\ Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} \]

因此，矩陣\(G_{k+1}\)的迭代公式：

\[G_{k+1}=G_k+\frac{\delta_k\delta_k^T}{\delta_k^Ty_k}-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} \]

能夠證實，只要初始矩陣\(G_0\)正定對稱，則迭代過程當中的每一個矩陣\(G_k\)均正定對稱，接下來對其進行代碼實現：

import numpy as np
"""
DPF擬牛頓法,封裝到ml_models.optimization模塊,與梯度降低法配合使用
"""


class DFP(object):
    def __init__(self, x0, g0):
        """

        :param x0: 初始的x
        :param g0: 初始x對應的梯度
        """
        self.x0 = x0
        self.g0 = g0
        # 初始化G0
        self.G0 = np.eye(len(x0))

    def update_quasi_newton_matrix(self, x1, g1):
        """
        更新擬牛頓矩陣
        :param x1:
        :param g1:
        :return:
        """
        # 進行一步更新
        y0 = g1 - self.g0
        delta0 = x1 - self.x0
        self.G0 = self.G0 + delta0.dot(delta0.T) / delta0.T.dot(y0)[0][0] - self.G0.dot(y0).dot(y0.T).dot(self.G0) / y0.T.dot(
            self.G0).dot(y0)[0][0]

    def adjust_gradient(self, gradient):
        """
        對原始的梯度作調整
        :param gradient:
        :return:
        """
        return self.G0.dot(gradient)

應用到LogisticRegression

咱們試一試將DFP算法應用到LogisticRegression，修改的地方以下：

fit函數追加以下的一個判斷：

elif self.solver == 'dfp':
    self.dfp = None
    self._fit_sgd(x, y)

_fit_sgd函數中，在梯度更新前作以下調整：

if self.solver == 'dfp':
    if self.dfp is None:
        self.dfp = optimization.DFP(x0=self.w, g0=dw)
    else:
        # 更新一次擬牛頓矩陣
        self.dfp.update_quasi_newton_matrix(self.w, dw)
    # 調整梯度方向
    dw = self.dfp.adjust_gradient(dw)

"""
梯度降低和DFP作一下對比
"""
from sklearn.datasets import make_classification
import matplotlib.pyplot as plt
%matplotlib inline

data, target = make_classification(n_samples=200, n_features=2, n_classes=2, n_informative=1, n_redundant=0,
                                   n_repeated=0, n_clusters_per_class=1)
target = target.reshape(200, 1)

import os
os.chdir('../')
from ml_models.linear_model import LogisticRegression

sgd_model = LogisticRegression(epochs=50)
sgd_model.fit(data, target)

dfp_model = LogisticRegression(solver='dfp',epochs=50)
dfp_model.fit(data,target)

#損失函數對比
plt.plot(range(0, len(sgd_model.losses)), sgd_model.losses,'b')
plt.plot(range(0, len(dfp_model.losses)), dfp_model.losses,'r')

[<matplotlib.lines.Line2D at 0x169fb529588>]

能夠發現，大部分狀況下DFP比SGD收斂的更快，且收斂效果更好

#分類效果對比
sgd_model.plot_decision_boundary(data,target)
dfp_model.plot_decision_boundary(data,target)

六.BFGS算法

BFGS算法是用一個矩陣\(B_k\)去模擬海瑟矩陣\(H_k\)，它的更新公式一樣假設有兩個附加項：

\[B_{k+1}=B_k+P_k+Q_k \]

固然，它須要知足擬牛頓條件：

\[B_{k+1}\delta_k=y_k \]

因此：

\[B_{k+1}\delta_k=B_k\delta_k+P_k\delta_k+Q_k\delta_k=y_k \]

考慮，使\(P_k\)和\(Q_k\)知足下面兩個條件：

\[P_k\delta_k=y_k\\ Q_k\delta_k=-B_k\delta_k \]

能夠獲得知足條件的解：

\[P_k=\frac{y_ky_k^T}{y_k^T\delta_k}\\ Q_k=-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} \]

因此更新公式：

\[B_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} \]

一樣能夠證實，若是\(B_0\)正定對稱，那麼迭代過程當中的每一個矩陣\(B_k\)都是正定對稱的，因爲這裏是對\(H_k\)的近似，因此每次更新梯度時，還須要對\(B_k\)作求逆計算，咱們可使用兩次以下的Sherman-Morrison公式：

\[(A+uu^T)^{-1}=A^{-1}-\frac{A^{-1}uu^TA^{-1}}{1+u^TA^{-1}u} \]

獲得BFGS算法關於\(G_k\)的迭代公式：

\[G_{k+1}=(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})G_k(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})^T+\frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \]

接下來，進行代碼實現：

"""
BFGS擬牛頓法,封裝到ml_models.optimization模塊,與梯度降低法配合使用
"""


class BFGS(object):
    def __init__(self, x0, g0):
        """

        :param x0: 初始的x
        :param g0: 初始x對應的梯度
        """
        self.x0 = x0
        self.g0 = g0
        # 初始化B0
        self.B0 = np.eye(len(x0))

    def update_quasi_newton_matrix(self, x1, g1):
        """
        更新擬牛頓矩陣
        :param x1:
        :param g1:
        :return:
        """
        # 進行一步更新
        y0 = g1 - self.g0
        delta0 = x1 - self.x0
        self.B0 = self.B0 + y0.dot(y0.T) / y0.T.dot(delta0)[0][0] - self.B0.dot(delta0).dot(delta0.T).dot(self.B0) / \
                                                                    delta0.T.dot(self.B0).dot(delta0)[0][0]

    def adjust_gradient(self, gradient):
        """
        對原始的梯度作調整
        :param gradient:
        :return:
        """
        return np.linalg.pinv(self.B0).dot(gradient)

應用到LogisticRegression

fit函數追加以下的一個判斷：

elif self.solver == 'bfgs':
    self.bfgs = None
    self._fit_sgd(x, y)

_fit_sgd函數中，在梯度更新前作以下調整：

if self.solver == 'bfgs':
    if self.bfgs is None:
        self.bfgs = optimization.BFGS(x0=self.w, g0=dw)
    else:
        # 更新一次擬牛頓矩陣
        self.bfgs.update_quasi_newton_matrix(self.w, dw)
    # 調整梯度方向
    dw = self.bfgs.adjust_gradient(dw)

#訓練模型
bfgs_model = LogisticRegression(solver='bfgs',epochs=50)
bfgs_model.fit(data,target)

#損失函數對比
plt.plot(range(0, len(sgd_model.losses)), sgd_model.losses,'b')
plt.plot(range(0, len(dfp_model.losses)), dfp_model.losses,'r')
plt.plot(range(0, len(bfgs_model.losses)), bfgs_model.losses,'y')

[<matplotlib.lines.Line2D at 0x169fd6e7b38>]

能夠發現大部分狀況下BFGS會比DFS收斂更快

#查看效果
bfgs_model.plot_decision_boundary(data,target)