图文详解牛顿迭代法,牛顿不止力学三定律

举报
Mr.Winter 发表于 2022/03/22 19:17:48 2022/03/22
【摘要】 图文详解牛顿迭代算法原理+手推公式,附Python实战代码加深理解

image.png


1 引例

给定如图所示的某个函数,如何计算函数零点 x 0 x_0


在这里插入图片描述

在数学上我们如何处理这个问题?

最简单的办法是解方程 f ( x ) = 0 f(x)=0 ,在代数学上还有著名的零点判定定理

如果函数 y = f ( x ) y= f(x) 在区间 [ a , b ] [a,b] 上的图象是连续不断的一条曲线,并且有 f ( a ) f ( b ) < 0 f(a)·f(b)<0 ,那么函数 y = f ( x ) y= f(x) 在区间 ( a , b ) (a,b) 内有零点,即至少存在一个 c ( a , b ) c∈(a,b) ,使得 f ( c ) = 0 f(c)=0 ,这个 c c 也就是方程 f ( x ) = 0 f(x)= 0 的根。

然而,数学上的方法并不一定适合工程应用,当函数形式复杂,例如出现超越函数形式;非解析形式,例如递推关系时,精确的方程解析一般难以进行,因为代数上还没发展出任意形式的求根公式。而零点判定定理求解效率也较低,需要不停试错。

因此,引入今天的主题——牛顿迭代法,服务于工程数值计算。

2 牛顿迭代算法求根

记第 k k 轮迭代后,自变量更新为 x k x_k ,令目标函数 f ( x ) f(x) x = x k x=x_k 泰勒展开:

f ( x ) = f ( x k ) + f ( x k ) ( x x k ) + o ( x ) f\left( x \right) =f\left( x_k \right) +f'\left( x_k \right) \left( x-x_k \right) +o(x)

我们希望下一次迭代到根点,忽略泰勒余项,令 f ( x k + 1 ) = 0 f(x_{k+1})=0 ,则

x k + 1 = x k f ( x k ) f ( x k ) x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}

不断重复运算即可逼近根点。

在几何上,上面过程实际上是在做 f ( x ) f(x) x = x k x=x_k 处的切线,并求切线的零点,在工程上称为局部线性化。如图所示,若 x k x_k x 0 x_0 的左侧,那么下一次迭代方向向右。

在这里插入图片描述

x k x_k x 0 x_0 的右侧,那么下一次迭代方向向左。

在这里插入图片描述

3 牛顿迭代优化

将优化问题转化为求目标函数一阶导数零点的问题,即可运用上面说的牛顿迭代法。

具体地,记第 k k 轮迭代后,自变量更新为 x k x_k ,令目标函数 f ( x ) f(x) x = x k x=x_k 泰勒展开:

f ( x ) = f ( x k ) + f ( x k ) ( x x k ) + 1 2 f ( x k ) ( x x k ) 2 + o ( x ) f\left( x \right) =f\left( x_k \right) +f'\left( x_k \right) \left( x-x_k \right) +\frac{1}{2}f''\left( x_k \right) \left( x-x_k \right) ^2+o(x)

两边求导得

f ( x ) = f ( x k ) + f ( x k ) ( x x k ) f'\left( x \right) =f'\left( x_k \right) +f''\left( x_k \right) \left( x-x_k \right)

f ( x k + 1 ) = f ( x k ) + f ( x k ) ( x k + 1 x k ) = 0 f'\left( x_{k+1} \right) =f'\left( x_k \right) +f''\left( x_k \right) \left( x_{k+1}-x_k \right) =0 ,从而得到

x k + 1 = x k f ( x k ) f ( x k ) x_{k+1}=x_k-\frac{f'\left( x_k \right)}{f''\left( x_k \right)}

对于向量 x = [ x 1 x 2 x d ] T \boldsymbol{x}=\left[ \begin{matrix} x_1& x_2& \cdots& x_d\\\end{matrix} \right] ^T ,将上述迭代公式推广为

x k + 1 = x k [ 2 f ( x k ) ] 1 f ( x k ) {\boldsymbol{x}_{k+1}=\boldsymbol{x}_k-\left[ \nabla ^2f\left( \boldsymbol{x}_k \right) \right] ^{-1}\nabla f\left( \boldsymbol{x}_k \right) }

其中 2 f ( x k ) \nabla ^2f\left( \boldsymbol{x}_k \right) Hessian矩阵,当其正定时可以保证牛顿优化算法往 减小的方向迭代

牛顿法的特点如下:

① 以二阶速率向最优点收敛,迭代次数远小于梯度下降法,优化速度快;

梯度下降法的解析参考图文详解神秘的梯度下降算法原理(附Python代码)

②学习率为 [ 2 f ( x k ) ] 1 \left[ \nabla ^2f\left( \boldsymbol{x}_k \right) \right] ^{-1} ,包含更多函数本身的信息,迭代步长可实现自动调整,可视为自适应梯度下降算法;

③ 耗费CPU计算资源多,每次迭代需要计算一次Hessian矩阵,且无法保证Hessian矩阵可逆且正定,因而无法保证一定向最优点收敛。

在实际应用中,牛顿迭代法一般不能直接使用,会引入改进来规避其缺陷,称为拟牛顿算法簇,其中包含大量不同的算法变种,例如共轭梯度法、DFP算法等等,今后都会介绍到。

4 代码实战: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中加载指定数据
* @param[in]: file -> 文件名
* @param[in]: colName -> 要加载的列名
* @param[in]: mode -> 加载模式, set: 列名与该列数据组成的字典, 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数据
    # ============================
    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('对数几率回归模拟\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='好瓜')

    # ============================
    # 实例化对数几率回归模型
    # ============================
    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()

其中更新权重代码为

    '''
    * @breif: 牛顿迭代法更新权重
    * @param[in]: None
    * @retval: 优化参数的增量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

在这里插入图片描述


🔥 更多精彩专栏

🏠 欢迎加入社区和更多志同道合的朋友交流AI 技术社

👇配套代码 · 优质体验 · 系统知识 请关注👇


【版权声明】本文为华为云社区用户原创内容,未经允许不得转载,如需转载请自行联系原作者进行授权。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

举报
请填写举报理由
0/200