• 在上一篇文章《機器學習算法推導&實現(xiàn)——邏輯斯蒂回歸》中禾酱，我們分別使用了梯度下降法和牛頓法來求解對數(shù)似然函數(shù)的最優(yōu)化問題微酬，從實例中我們也能發(fā)現(xiàn)牛頓法的收斂速度遠快于梯度下降法，但是在上文中颤陶，直接使用了牛頓法的結(jié)果佩伤，并沒有進行相應推導蕊温，故本文一方面是補上牛頓法的推導，另一方面是展開討論下擬牛頓法。

牛頓法推導

• 牛頓法的推導得先從泰勒公式說起：

  
  
  0.PNG

• 假設現(xiàn)在函數(shù)f(x)迭代了k次的值為X(k)操灿，則在X(k)上進行二階泰勒展開可近似得到以下公式：

  
  
  1.PNG

• 我們要求得f(x)的極小值，則必要條件是f(x)在極值點處的一階導數(shù)為0盟蚣，即：

  
  
  2.PNG

• 因為我們把每輪迭代求得的滿足目標函數(shù)極小值的x作為下一輪迭代的值刺覆，因此我們可以假設第k+1輪的值就是最優(yōu)解：

  
  
  3.PNG
• 代入二階泰勒展開并求導可得：
  
  4.PNG
  
  令：
  
  5.PNG
  6.PNG
  
  可得最終的優(yōu)化公式為：
  
  7.PNG
• 雖然根據(jù)推導，我們已經(jīng)知道了牛頓法的迭代方法姨涡，但是在實際應用過程中衩藤，我們會發(fā)現(xiàn)海塞矩陣的逆矩陣往往計算比較復雜，于是又有了擬牛頓法來簡化這一過程涛漂。

擬牛頓法的原理

• 在擬牛頓法中赏表，我們考慮優(yōu)化出一個n階矩陣D來代替海塞矩陣的逆矩陣，首先我們得先來看看替代矩陣D要滿足什么條件：

• 首先根據(jù)二階泰勒展開匈仗，我們令：

  
  
  8.PNG

• 代入f(x)瓢剿，并求導，可得：

  
  
  9.PNG

• 再令：

  
  
  5.PNG
  6.PNG

• 最終可得矩陣D需滿足的條件：

  
  
  10.PNG

• 另外在每次迭代中可以更新矩陣D悠轩，故可以假設以下更新條件：

  
  
  11.PNG
  
  

  由此我們可以發(fā)現(xiàn)海塞矩陣逆矩陣的近似矩陣D(x)的選擇條件比較靈活间狂，可以有多種具體的實現(xiàn)方法。

1火架、DFP算法（Davidon-Fletcher-Powell）推導

• 在DFP算法中鉴象，我們假設：
  
  12.PNG
  
  *在這里我們知道凸函數(shù)的二階導數(shù)海塞矩陣(Hessian)是對稱矩陣，因此我們假設想要替代的D, P, Q矩陣也是對稱矩陣何鸡，即矩陣的轉(zhuǎn)置等于矩陣本身炼列。
• 兩邊乘以可得：
  
  13.PNG

• 為了滿足上式要求，我們不妨令：

  
  
  14.PNG
• 我們先來求P矩陣音比，最簡單的就是令：
  
  16.PNG
• 接著兩邊乘以可得：
  
  18.PNG
  
  求得俭尖，并代入P矩陣，可求得：
  
  19.PNG
• 同樣的方法我們求解下Q矩陣：
  
  20.PNG
• 最終，DFP算法的迭代公式為：
  
  21.PNG

python實現(xiàn)DFP算法

• 我們繼承上一篇文章邏輯回歸算法的類稽犁，并且增加擬牛頓法的優(yōu)化算法焰望。
• 我們先看下算法的主函數(shù)部分买决，在繼承邏輯回歸類的基礎上主要有2點改動：一是把"method"參數(shù)提前齐苛，在實例化的時候就要求定位好優(yōu)化算法；二是重寫了train方法承二，衍生出其他擬牛頓算法虑椎。

from ML_LogisticRegression import LogisticRegressionSelf as logit
import numpy as np
?
class allLogitMethod(logit):
    def __init__(self, method):
        """init class"""
        super().__init__()
        self.method = method
        self.graList = []               #記錄梯度模長的列表
?
    #重寫下train方法
    def train(self, X, y, n_iters=1000, learning_rate=0.01):
        """fit model"""
        if X.ndim < 2:
            raise ValueError("X must be 2D array-like!")
        self.trainSet = X
        self.label = y
        if self.method.lower() == "gradient":
            self._LogisticRegressionSelf__train_gradient(n_iters, learning_rate)
        elif self.method.lower() == "newton":
            self._LogisticRegressionSelf__train_newton(n_iters)
        elif self.method.lower() == "dfp":
            self.__train_dfp(n_iters, learning_rate)
        else:
            raise ValueError("method value not found!")
        return  

• 在寫算法主體前震鹉，我們先看下DFP的簡易優(yōu)化過程：

1. 初始化參數(shù)W0，初始化替代矩陣Dk0捆姜，計算初始梯度Gk0传趾，迭代次數(shù)k;
2. While ( k < n_iters ) and ( ||Gk0|| > e ):
   1. 計算出W0需要優(yōu)化的方向Pk = Dk*Gk0
   2. 更新參數(shù)W1如下：
      for n in N:
      • 嘗試用一維搜索方法改變Pk的學習率，
      • 求得最小似然值泥技，
      • 求出W1和deltaW（W1 - W0）
   3. 更新梯度Gk1浆兰，并求deltaG（Gk1 - Gk0）
   4. 更新DK1，利用上述推導的DFP的迭代公式
   5. 重新賦值W0珊豹，Dk0簸呈，Gk0
   6. k += 1
3. end.

• 我們可以大致看到DFP算法的主體部分有外循環(huán)和內(nèi)循環(huán)兩個部分組成，所以我們也是將內(nèi)外循環(huán)分開來寫店茶。
• 我們先看內(nèi)循環(huán)部分蜕便。內(nèi)循環(huán)是已經(jīng)知道了參數(shù)W的優(yōu)化方向的基礎上，想要找到最合適的學習率贩幻，使得更新后的W求得的似然值最小玩裙。實際有很多方法，在這里為了簡單段直，筆者只是對學習率進行了i次方的操作。具體函數(shù)如下：

    #一維搜索法求出最優(yōu)lambdak溶诞，更新W后鸯檬，使得似然值最小
    def __updateW(self, X, Y, lambdak, W0, Pk):
        """此處對lambdak的處理僅簡單用1~i次方來逐步變小，以選取到最小似然值的lambdak"""
        min_LLvalue = np.inf
        W1 = np.zeros(W0.shape)
        for i in range(10):
            Wi = W0 - (lambdak**i)*Pk
            Ypreprob, LLvalue = self.PVandLLV(X, Y, Wi)
            if LLvalue < min_LLvalue:
                min_LLvalue = LLvalue
                W1 = np.copy(Wi)
                deltaW = - (lambdak**i)*Pk
                bestYpreprob = Ypreprob
        return W1, deltaW, min_LLvalue, bestYpreprob

• 再看外循環(huán)部分螺垢。

    #新增擬牛頓法-DFP優(yōu)化算法
    def __train_dfp(self, n_iters, learning_rate):
        """Quasi-Newton Method DFP（Davidon-Fletcher-Powell）"""
        n_samples, n_features = self.trainSet.shape
        X = self.trainSet
        y = self.label
        #合并w和b喧务，在X尾部添加一列全是1的特征
        X2 = np.hstack((X, np.ones((n_samples, 1))))
        #將y轉(zhuǎn)置變?yōu)?n_samples,1)的矩陣
        Y = np.expand_dims(y, axis=1)
        #初始化特征系數(shù)W，初始化替代對稱矩陣
        W = np.zeros((1, n_features+1))
        Dk0 = np.eye(n_features+1)
        #計算初始的預測值枉圃、似然值功茴，并記錄似然值
        Ypreprob, LL0 = self.PVandLLV(X2, Y, W)
        self.llList.append(LL0)
        #根據(jù)初始的預測值計算初始梯度，并記錄梯度的模長
        Gk0 = self._LogisticRegressionSelf__calGradient(X2, Y, Ypreprob)
        graLength = np.linalg.norm(Gk0)
        self.graList.append(graLength)
        #初始化迭代次數(shù)
        k = 0
        while (k<n_iters) and (graLength>self.tol):
            #計算優(yōu)化方向的值Pk=Gk0.Dk0
            Pk = np.dot(Gk0, Dk0)
            #一維搜索更新參數(shù)孽亲，并保存求得的最小似然值
            W, deltaW, min_LLvalue, Ypreprob = self.__updateW(X2, Y, learning_rate, W, Pk)
            self.llList.append(min_LLvalue)
            #更新梯度Gk和deltaG坎穿，同時求得梯度的模長和更新前后的模長差值
            Gk1 = self._LogisticRegressionSelf__calGradient(X2, Y, Ypreprob)
            graLength = np.linalg.norm(Gk1)
            self.graList.append(graLength)
            deltaG = Gk1 - Gk0
            Gk0 = Gk1
            #更新替代矩陣Dk
            Dk1 = Dk0 + np.dot(deltaW.T, deltaW)/np.dot(deltaW, deltaG.T) - \
            np.dot(np.dot(np.dot(Dk0, deltaG.T), deltaG), Dk0)/np.dot(np.dot(deltaG, Dk0), deltaG.T)
            Dk0 = Dk1
            k += 1
        self.n_iters = k
        self.w = W.flatten()[:-1]
        self.b = W.flatten()[-1]
        Ypre = np.argmax(np.column_stack((1-Ypreprob,Ypreprob)), axis=1)
        self.accurancy = sum(Ypre==y)/n_samples
        print("第{}次停止迭代，梯度模長為{}，似然值為{}玲昧，準確率為{}".format(self.n_iters, self.graList[-1], self.llList[-1], self.accurancy))
        print("w:{};\nb:{}".format(self.w, self.b))
        return

• 在這里還有一點小小的提示：筆者在實測過程中栖茉，發(fā)現(xiàn)DFP算法經(jīng)常會導致求解似然值和P(y=1|X)的時候報值溢出的錯誤。主要是python的numpy.exp(wx)中的wx過大導致的孵延，比如看下圖的報錯：
  
  22.PNG
• 所以在這里我們對求解似然值和P(y=1|X)的公式進行了一些調(diào)整吕漂。當wx為負數(shù)時，使用原公式尘应；當wx為正數(shù)時惶凝，使用轉(zhuǎn)換公式替換。其實質(zhì)是一致的犬钢。
  
  23.PNG
• 最后我們用擬牛頓法-DFP算法和牛頓法實測比較下苍鲜。
  1. 牛頓法：迭代了7次，迭代時長0.18s娜饵，似然值176.81坡贺，準確率0.94，模型參數(shù)如下圖所示：

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    from sklearn.datasets import make_classification
    import time
    X, y = make_classification(n_samples=1000, n_features=4)
    #1箱舞、自編的牛頓法進行擬合
    time_nt = time.time()
    logit_nt = allLogitMethod("newton")
    logit_nt.train(X, y, n_iters=100)
    print("迭代時長：", time.time()-time_nt)
    plt.plot(range(logit_nt.n_iters+1), logit_nt.llList)
    plt.show()

24.PNG

    #3、自編的擬牛頓法-DFP算法進行擬合
    time_dfp = time.time()
    logit_dfp = allLogitMethod("DFP")
    logit_dfp.train(X, y, n_iters=20, learning_rate=0.5)
    print("迭代時長：", time.time()-time_dfp)
    fig = plt.figure(figsize=(10,6)) 
    ax1 = fig.add_subplot(1,2,1)
    ax1.plot(range(logit_dfp.n_iters+1), logit_dfp.llList)
    ax2 = fig.add_subplot(1,2,2)
    ax2.plot(range(logit_dfp.n_iters+1), logit_dfp.graList)
    plt.show()

25.PNG

• 經(jīng)過兩種方法的對比电湘，我們發(fā)現(xiàn)兩者最終訓練出來的模型參數(shù)基本一致隔节，似然值也下降到同一水平；DFP雖然迭代18次比牛頓法多寂呛，但是訓練總時間只有0.038s怎诫，遠小于牛頓法的訓練總時間。
• 可見贷痪，DFP算法中的Dk矩陣完全可以逼近于牛頓法中的海塞矩陣的逆矩陣幻妓，而且DFP算法中的Dk矩陣的訓練也比二階導數(shù)矩陣的訓練來得方便與快速。

以上就是本次的全部內(nèi)容劫拢，謝謝閱讀肉津。（今天先寫到這，寫的好累舱沧，擬牛頓法的其他算法我們下次再接著推導與實現(xiàn)妹沙。）

全部代碼可前往以下地址下載：
https://github.com/shoucangjia1qu/ML_gzh/tree/master/LogisticRegression

?往期回顧：
機器學習算法推導&實現(xiàn)——邏輯斯蒂回歸
機器學習算法推導&實現(xiàn)——線性回歸