一.簡介
通過前面幾節的介紹，大家可以直觀的感受到：對於大部分機器學習模型，我們通常會將其轉化為一個優化問題，由於模型通常較為復雜，難以直接計算其解析解，我們會采用迭代式的優化手段，用數學語言描述如下：
這里目標函數為\(f(x)\)，當前優化變量為\(v^k\)，目標即是找到一個\(v^k\)對當前的\(x^k\)進行更新，使得函數值盡可能的降低，如果目標函數一階可微，對其作一階泰勒展開即可得到如下梯度下降的更新公式：
二.梯度下降
對目標函數作一階泰勒展開：
所以要使得\(f(x^k+v^k)<f(x^k)\)只需要使\(\triangledown f(x^k)^Tv^k<0\)即可，而如果取：
則一定能使\(\triangledown f(x^k)^Tv^k<0\)，所以，我們就得到了梯度下降的更新公式：
這里\(\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\)的取值大到不能忽略后面的項時，泰勒展開的二階項取值就必須要考慮其中了，所以這時梯度下降法未必時最優的方向，接下來我們看下二階展開的情況：
三.牛頓法
對其作二階泰勒展開：
這里為了方便，記\(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)\)半正定，即對於\(\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處取得：
可以發現牛頓法對比梯度下降法，其實牛頓法是對梯度法的方向進行了一個改變\(H_k^{-1}\)，所以，我們可以得到牛頓法的更新公式：
這里記\(p_k=-H_k^-1g_k\)；
可以發現牛頓法的復雜有點高，因為要求解\(H_k^{-1}\)，那么有沒有方便一點的方法呢？比如構建一個矩陣去近似\(H_k\)或者\(H_k^{-1}\)，這便是擬牛頓法的基本思想
四.擬牛頓條件
上面說到了利用一個矩陣去近似Hessian矩陣或者Hessian矩陣的逆，那么這個近似矩陣需要滿足怎樣的條件呢？我們還是從二階泰勒展開出發，稍微變換一下：
兩邊對\(x^{k+1}\)求偏導可得：
這便是擬牛頓條件，為了方便，記\(y_k=g_{k+1}-g_k,\delta_k=x^{k+1}-x^k\)，所以：
所以，擬牛頓法也要滿足和\(H_k\)一樣的性質：
（1）正定性；
（2）滿足擬牛頓條件
接下來，簡單證明一下如果滿足性質（1）：正定性，更新時可以滿足函數值下降，假設\(G_k\)去近似\(H_k^{-1}\)，所以：\(G_k\succ 0\)，那么迭代公式為：
將其帶入二階泰勒展開式中：
通常\(\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\)加上兩個附加項構成的，即：
這里\(P_k,Q_k\)是待定項，由於需要滿足擬牛頓條件，所以：
這里做一個tricky的假設：
這樣的\(P_k,Q_k\)不難找到：
所以，矩陣\(G_{k+1}\)的迭代公式：
可以證明，只要初始矩陣\(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\)，它的更新公式同樣假設有兩個附加項：
當然，它需要滿足擬牛頓條件：
所以：
考慮，使\(P_k\)和\(Q_k\)滿足下面兩個條件：
可以得到滿足條件的解：
所以更新公式：
同樣可以證明，如果\(B_0\)正定對稱，那么迭代過程中的每個矩陣\(B_k\)都是正定對稱的，由於這里是對\(H_k\)的近似，所以每次更新梯度時，還需要對\(B_k\)做求逆計算，我們可以使用兩次如下的Sherman-Morrison公式：
得到BFGS算法關於\(G_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)
