数值计算之 牛顿法与函数极值

  • 前言
  • 最速下降法
  • 牛顿法
  • 牛顿法分析
  • 代码示例
  • 后记

前言

本篇继续优化理论的算法学习,牛顿法。

最速下降法

首先回顾上次提到的梯度下降法(其实就是最速下降法):通过求取多元函数在某个点处的梯度,沿着梯度的反方向前进,直到达到迭代停止条件。

对于多元函数(实值向量函数)f(x)f(\bf x)f(x),其在x0\bf x_0x0​处的泰勒展开可表示为:
f(x)=f(x0)+∇f(x0)T(x−x0)+12(x−x0)TH(x0)(x−x0)f({\bf x}) = f({\bf x_0})+{\bf \nabla} f({\bf x_0})^T({\bf x-x_0}) + \frac {1}{2} ({\bf x-x_0})^T{\bf H(x_0)}({\bf x-x_0}) f(x)=f(x0​)+∇f(x0​)T(x−x0​)+21​(x−x0​)TH(x0​)(x−x0​)
其中,∇,H\nabla, H∇,H分别是函数梯度,海森矩阵。

梯度下降法直接取一阶梯度的反方向作为优化方向,因此称为最速下降法(每次迭代的方向都是下降最快的方向)。

牛顿法

回到多元泰勒展开:
f(x)=f(x0)+∇f(x0)(x−x0)+12(x−x0)TH(x0)(x−x0)f({\bf x}) = f({\bf x_0})+{\bf \nabla} f({\bf x_0})({\bf x-x_0}) + \frac {1}{2} ({\bf x-x_0})^T{\bf H(x_0)}({\bf x-x_0}) f(x)=f(x0​)+∇f(x0​)(x−x0​)+21​(x−x0​)TH(x0​)(x−x0​)
f(x)f(\bf x)f(x)的极小值处的导数应当等于0\bf 00。现在取初始点x0\bf x_0x0​,将f(x)f(\bf x)f(x)在x0\bf x_0x0​处的展开式对x\bf xx进行求导:
∇f(x0)T+(x−x0)TH(x0)=0∇f(x0)=−H(x0)(x−x0)x=x0−H(x0)−1∇f(x0){\bf \nabla} f({\bf x_0})^T+{(\bf x-x_0)}^T{\bf H(x_0)}={\bf 0} \\ {\bf \nabla}f({\bf x_0})=-{\bf H(x_0)}({\bf x-x_0}) \\ \quad \\ {\bf x} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}} ∇f(x0​)T+(x−x0​)TH(x0​)=0∇f(x0​)=−H(x0​)(x−x0​)x=x0​−H(x0​)−1∇f(x0​)
这样就得到了新的x\bf xx,这就是牛顿法。

可以给出更具体的牛顿法求极值过程:

  1. 确定实值向量函数f(x)f(\bf x)f(x)的初始点x0\bf x_0x0​,梯度∇f{\bf \nabla} f∇f的表达式,海森矩阵H(x)\bf H(x)H(x)的表达式,终止条件δ,ϵ\delta, \epsilonδ,ϵ
  2. 计算梯度∇f(x0){\bf \nabla}f({\bf x_0})∇f(x0​),如果∣∣∇f(x0)∣∣<δ||{\bf \nabla}f({\bf x_0})||<\delta∣∣∇f(x0​)∣∣<δ,终止迭代,得到解x=x0\bf x=x_0x=x0​;否则进入第3步
  3. 计算海森矩阵,并迭代x′=x0−H(x0)−1∇f(x0){\bf x'} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}}x′=x0​−H(x0​)−1∇f(x0​)
  4. 计算Δf=f(x′)−f(x0)\Delta f= f({\bf x'}) - f(\bf x_0)Δf=f(x′)−f(x0​),如果Δf<ϵ\Delta f < \epsilonΔf<ϵ,终止迭代,得到解x=x′\bf x=x'x=x′;否则x0=x′\bf x_0=x'x0​=x′,回到第2步

牛顿法分析

现在考虑牛顿法的优化迭代表达式:
x′=x0−H(x0)−1∇f(x0){\bf x'} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}} x′=x0​−H(x0​)−1∇f(x0​)

将牛顿法的迭代方向与梯度方向做内积:
−H(x0)−1∇f(x0)⋅∇f(x0)=−∇f(x0)TH(x0)−T∇f(x0)=−∇f(x0)TH(x0)−1∇f(x0)-{{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}} \cdot {{\bf \nabla} f{\bf(x_0)}} \\ = -{{\bf \nabla} f{\bf(x_0)}}^T {{\bf H(x_0)}}^{-T}{{\bf \nabla} f{\bf(x_0)}} \\ = -{{\bf \nabla} f{\bf(x_0)}}^T {{\bf H(x_0)}}^{-1}{{\bf \nabla} f{\bf(x_0)}} \\ −H(x0​)−1∇f(x0​)⋅∇f(x0​)=−∇f(x0​)TH(x0​)−T∇f(x0​)=−∇f(x0​)TH(x0​)−1∇f(x0​)
当海森矩阵是正定矩阵时,迭代方向与梯度方向的内积为负,即函数值在不断减小。这就涉及到了优化理论的重要问题:什么条件能够满足海森矩阵的正定性?

牛顿法的优点:优化收敛速度比梯度下降法更快。原因是牛顿法不仅考虑到函数的一阶梯度,也考虑到了海森矩阵(二阶梯度),在每个优化点处附近都能一次求得最佳迭代值。

牛顿法的缺点:每次迭代的速度较慢。这是因为迭代值的计算中涉及到了海森矩阵的逆,使得单步迭代效率比较低。

代码示例

和上次一样,我写了一个同牛顿法寻找二元函数极小值的代码,二元五次函数,如下所示:

import numpy as npdef partial_derivate_xy(x, y, F):dx = (F(x + 0.001, y) - F(x, y))/0.001dy = (F(x, y + 0.001) - F(x, y))/0.001return dx, dydef partial_partial_derivate_xy(x, y, F):dx, dy = partial_derivate_xy(x, y, F)dxx = (partial_derivate_xy(x + 0.001, y, F)[0] - dx) / 0.001dyy = (partial_derivate_xy(x, y + 0.001, F)[1] - dy) / 0.001dxy = (partial_derivate_xy(x, y + 0.001, F)[0] - dx) / 0.001dyx = (partial_derivate_xy(x + 0.001, y, F)[1] - dy) / 0.001return dxx, dyy, dxy, dyxdef non_linear_func(x, y):fxy = 0.5 * (x ** 2 + y ** 2)return fxydef non_linear_func_2(x, y):fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2return fxydef non_linear_func_3(x, y):fxy = 0.5 * (x ** 2 - y ** 2)return fxydef non_linear_func_4(x, y):fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5return fxydef newton_optim(x, y, F):dx, dy = partial_derivate_xy(x, y, F)dxx, dyy, dxy, dyx = partial_partial_derivate_xy(x, y, F)grad = np.array([[dx], [dy]])hessian = np.array([[dxx, dxy], [dyx, dyy]])hessian_inv = np.linalg.inv(hessian)vec_delta = np.matmul(hessian_inv, grad)vec_opt = np.array([[x], [y]]) - vec_deltax_opt = vec_opt[0][0]y_opt = vec_opt[1][0]return x_opt, y_opt, vec_deltadef optimizer(x0, y0, F, th=0.01):x = x0y = y0counter = 0while True:x_opt, y_opt, vec_delta = newton_optim(x, y, F)if np.linalg.norm(vec_delta) < th:breakx = x_opty = y_optcounter = counter + 1print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y))return x, yif __name__ == '__main__':x0 = 2.y0 = 2.result_x, result_y = optimizer(x0, y0, non_linear_func_4)print(result_x, result_y)

牛顿法迭代11次后得到优化解。

梯度下降法在lr=0.01时迭代149次得到优化解。

运行1000次时,牛顿法用时0.586s,梯度下降法用时0.285s,可见虽然牛顿法迭代次数少,但速度慢。

后记

下篇将进入高斯牛顿法,后续还会有LM算法、拟牛顿法等。

数值计算之 牛顿法与函数极值相关推荐

  1. 数值计算之 梯度下降法与函数极值

    数值计算之 梯度下降法与函数极值 前言 微积分基础 一元函数的极值,导数与泰勒展开 多元函数的泰勒展开 梯度下降法 梯度方向 终止条件 代码举例 后记 前言 本篇将开始介绍优化算法.首先是梯度下降法, ...

  2. MATLAB求解二元(多元)函数极值

    matlab求解二元函数极值 依然是机房中的R2010a版本 命令: 1.x=fminsearch(fun,x0)或x=fminunc(fun,x0)求极小值点x,初值选为x0 2.[x,fmin]= ...

  3. Python实现多元线性回归方程梯度下降法与求函数极值

    梯度下降法 梯度下降法的基本思想可以类比为一个下山的过程. 假设这样一个场景:一个人被困在山上,需要从山上下来(找到山的最低点,也就是山谷).但此时山上的浓雾很大,导致可视度很低:因此,下山的路径就无 ...

  4. matlab解方程最值点,MATLAB解方程与函数极值

    1.线性方程数值求解 主要是用到了计算方法里的LU分解等不过是加快了求解速度而已相对于inv(A)*b或者A\b 2.非线性方程数值求解 1 单变量非线性方程求解 在MATLAB中提供了一个fzero ...

  5. 基于Matlab的神经网络结合遗传算法在非线性函数极值寻优中的应用

    本微信图文利用神经网络进行非线性函数数据的拟合并通过遗传算法对训练后的神经网络进行非线性函数极值寻优.

  6. 如何利用神经网络结合遗传算法进行非线性函数极值寻优(2)

    如何利用神经网络结合遗传算法进行非线性函数极值寻优

  7. python求函数极值_python 遗传算法求函数极值的实现代码

    废话不多说,大家直接看代码吧! """遗传算法实现求函数极大值-Zjh""" import numpy as np import rando ...

  8. MATLAB学习笔记(七)——MATLAB解方程与函数极值

    (一)线性方程组求解 包含n个未知数,由n个方程构成的线性方程组为: 其矩阵表示形式为: 其中 一.直接求解法 1.左除法 x=A\b; 如果A是奇异的,或者接近奇异的.MATLAB会发出警告信息的. ...

  9. 运用遗传算法求解函数极值(fortran)

    运用遗传算法求解函数极值(fortran) 写在前面 遗传算法的前世今生 算法步骤简介 遗传算法的主体结构 开始求解: 结果显示: 最后再来说一些需要注意的地方 写在前面 这篇文章适合一些应急学习最优 ...

最新文章

  1. POJ - 1330 Nearest Common Ancestors tanjan_LCA
  2. 主营无线部件 高通与TDK创立合资公司
  3. mac 上的环境变量配置
  4. 【Cocosd2d实例教程八】Cocos2d实现碰撞检测(含实例)
  5. Eclipse中将Java项目转换成Web项目的方法
  6. oracle报V27的错误解决办法,oracle11g ora-27154 past/wait 错误解决方法
  7. mysql组件化_MySql笔记
  8. android led灯框架_LED面板灯的特点:应用领域、产品结构与产品分类
  9. Python编程基础13:文件读写操作
  10. java m.find()_正则有关问题 m.groupCount() 和 m.find()
  11. Android Multimedia框架总结(四)MediaPlayer中从Java层到C++层类关系及prepare及之后其他过程
  12. react引入html2canvas和jspdf生成PDF打印及下载
  13. isbn书号查询php代码,php根据isbn书号查询amazon网站上的图书信息的示例_PHP
  14. android 释放摄像头,android – 为什么Camera需要在onPause()而不是onstop()方法中释放?...
  15. express 配置ip
  16. mysql netbeans_使用Netbeans操作MySQL数据库
  17. 线程和进程的区别 线程和进程有什么不同
  18. RAID磁盘阵列的几种模式
  19. Excel的简单编程
  20. idea中的Diagram功能,查看类图

热门文章

  1. python爬b站评论_学习笔记(1):写了个python爬取B站视频评论的程序
  2. 一款DCDC稳压开关电源转换器HM1509
  3. 嵌入式Linux开发调优之一:系统与内核
  4. 安卓安装kali教程
  5. Java模拟生意参谋登录_用Excel实现生意参谋爬虫,伪装登陆状态
  6. 百度翻译参数逆向过程
  7. 【数据结构】哈希表的概念及应用
  8. MSN Media协议分析
  9. 保时捷纯电自动驾驶跑车上路,翻版特斯拉?
  10. 2006世界杯32强人体彩绘队服样式(捷克)