圖文詳解牛頓迭代算法原理及Python實現(xiàn)

更新時間：2022年08月10日 14:17:15 作者：Mr.Winter`

牛頓迭代法又稱為牛頓-拉夫遜（拉弗森）方法，它是牛頓在17世紀(jì)提出的一種在實數(shù)域和復(fù)數(shù)域上近似求解方程的方法。本文將利用圖文詳解牛頓迭代算法原理及實現(xiàn)，需要的可以參考一下

1.引例

給定如圖所示的某個函數(shù)，如何計算函數(shù)零點x₀

在數(shù)學(xué)上我們?nèi)绾翁幚磉@個問題？

最簡單的辦法是解方程f(x)=0，在代數(shù)學(xué)上還有著名的零點判定定理

如果函數(shù)y=f(x)在區(qū)間[a,b]上的圖象是連續(xù)不斷的一條曲線，并且有f(a)⋅f(b)<0,那么函數(shù)y=f(x)在區(qū)間(a,b)內(nèi)有零點，即至少存在一個c∈(a,b),使得f(c)=0,這個c也就是方程f(x)=0的根。

然而，數(shù)學(xué)上的方法并不一定適合工程應(yīng)用，當(dāng)函數(shù)形式復(fù)雜，例如出現(xiàn)超越函數(shù)形式；非解析形式，例如遞推關(guān)系時，精確的方程解析一般難以進(jìn)行，因為代數(shù)上還沒發(fā)展出任意形式的求根公式。而零點判定定理求解效率也較低，需要不停試錯。

因此，引入今天的主題——牛頓迭代法，服務(wù)于工程數(shù)值計算。

2.牛頓迭代算法求根

記第k輪迭代后，自變量更新為x_k，令目標(biāo)函數(shù)f(x)在x=x_k泰勒展開：

f(x)=f(x_k?)+f′(x_k?)(x−x_k?)+o(x)

我們希望下一次迭代到根點，忽略泰勒余項，令f(x_k+1)=0，則

x_k+1?=x_k?−f(x_k?)/f'(x_k?)?

不斷重復(fù)運(yùn)算即可逼近根點。

在幾何上，上面過程實際上是在做f(x)在x=x_k處的切線，并求切線的零點，在工程上稱為局部線性化。如圖所示，若x_k在x₀的左側(cè)，那么下一次迭代方向向右。

若x_k在x₀的右側(cè)，那么下一次迭代方向向左。

3.牛頓迭代優(yōu)化

將優(yōu)化問題轉(zhuǎn)化為求目標(biāo)函數(shù)一階導(dǎo)數(shù)零點的問題，即可運(yùn)用上面說的牛頓迭代法。

具體地，記第k輪迭代后，自變量更新為x_k ，令目標(biāo)函數(shù)f(x)在x=x_k泰勒展開：

f(x)=f(x_k?)+f′(x_k?)(x−x_k?)+1/2?f′′(x_k?)(x−x_k?)²+o(x)

兩邊求導(dǎo)得

f′(x)=f′(x_k?)+f′′(x_k?)(x−x_k?)

令f′(x_k+1?)=f′(x_k?)+f′′(x_k?)(x_k+1?−x_k?)=0，從而得到

x_k+1?=x_k?−f′(x_k?)/f'′(x_k?)?

對于向量x=[x₁?? x_2????x_d??]^T，將上述迭代公式推廣為

x_k+1?=x_k?−[∇²f(x_k?)]⁻¹∇f(x_k?)

其中∇²f(x_k?)是Hessian矩陣，當(dāng)其正定時可以保證牛頓優(yōu)化算法往減小的方向迭代

牛頓法的特點如下：

① 以二階速率向最優(yōu)點收斂，迭代次數(shù)遠(yuǎn)小于梯度下降法，優(yōu)化速度快；

梯度下降法的解析參考圖文詳解梯度下降算法的原理及Python實現(xiàn)

②學(xué)習(xí)率為[∇²f(x_k?)]⁻¹ ，包含更多函數(shù)本身的信息，迭代步長可實現(xiàn)自動調(diào)整，可視為自適應(yīng)梯度下降算法；

③ 耗費CPU計算資源多，每次迭代需要計算一次Hessian矩陣，且無法保證Hessian矩陣可逆且正定，因而無法保證一定向最優(yōu)點收斂。

在實際應(yīng)用中，牛頓迭代法一般不能直接使用，會引入改進(jìn)來規(guī)避其缺陷，稱為擬牛頓算法簇，其中包含大量不同的算法變種，例如共軛梯度法、DFP算法等等，今后都會介紹到。

4 代碼實戰(zhàn)：Logistic回歸

import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt
import matplotlib as mpl
from Logit import Logit

'''
* @breif: 從CSV中加載指定數(shù)據(jù)
* @param[in]: file -> 文件名
* @param[in]: colName -> 要加載的列名
* @param[in]: mode -> 加載模式, set: 列名與該列數(shù)據(jù)組成的字典, df: df類型
* @retval: mode模式下的返回值
'''
def loadCsvData(file, colName, mode='df'):
    assert mode in ('set', 'df')
    df = pd.read_csv(file, encoding='utf-8-sig', usecols=colName)
    if mode == 'df':
        return df
    if mode == 'set':
        res = {}
        for col in colName:
            res[col] = df[col].values
        return res

if __name__ == '__main__':
    # ============================
    # 讀取CSV數(shù)據(jù)
    # ============================
    csvPath = os.path.abspath(os.path.join(__file__, "../../data/dataset3.0alpha.csv"))
    dataX = loadCsvData(csvPath, ["含糖率", "密度"], 'df')
    dataY = loadCsvData(csvPath, ["好瓜"], 'df')
    label = np.array([
        1 if i == "是" else 0
        for i in list(map(lambda s: s.strip(), list(dataY['好瓜'])))
    ])

    # ============================
    # 繪制樣本點
    # ============================
    line_x = np.array([np.min(dataX['密度']), np.max(dataX['密度'])])
    mpl.rcParams['font.sans-serif'] = [u'SimHei']
    plt.title('對數(shù)幾率回歸模擬\nLogistic Regression Simulation')
    plt.xlabel('density')
    plt.ylabel('sugarRate')
    plt.scatter(dataX['密度'][label==0],
                dataX['含糖率'][label==0],
                marker='^',
                color='k',
                s=100,
                label='壞瓜')
    plt.scatter(dataX['密度'][label==1],
                dataX['含糖率'][label==1],
                marker='^',
                color='r',
                s=100,
                label='好瓜')

    # ============================
    # 實例化對數(shù)幾率回歸模型
    # ============================
    logit = Logit(dataX, label)

    # 采用牛頓迭代法
    logit.logitRegression(logit.newtomMethod)
    line_y = -logit.w[0, 0] / logit.w[1, 0] * line_x - logit.w[2, 0] / logit.w[1, 0]
    plt.plot(line_x, line_y, 'g-', label="牛頓迭代法")

    # 繪圖
    plt.legend(loc='upper left')
    plt.show()

其中更新權(quán)重代碼為

    '''
    * @breif: 牛頓迭代法更新權(quán)重
    * @param[in]: None
    * @retval: 優(yōu)化參數(shù)的增量dw
    '''
    def newtomMethod(self):
        wTx = np.dot(self.w.T, self.X).reshape(-1, 1)
        p = Logit.sigmod(wTx)
        dw_1 = -self.X.dot(self.y - p)
        dw_2 = self.X.dot(np.diag((p * (1 - p)).reshape(self.N))).dot(self.X.T)
        dw = np.linalg.inv(dw_2).dot(dw_1)
        return dw