牛顿法与拟牛顿法

牛顿法和拟牛顿法是求解无约束最优化的常用方法,有收敛速度快的优点. 牛顿法属于迭代算法,每一步需要求解目标函数的海赛矩阵的逆矩阵,计算复杂. 拟牛顿法通过正定矩阵近似海赛矩阵的逆矩阵,简化了这个过程.

牛顿法

对于无约束优化 min⁡x∈Rnf(x) \min_{x\in R^n} f(x) x∈Rnmin​f(x) x∗x^*x∗是目标的极小值点.

假设f(x)f(x)f(x)有二阶连续偏导数,第k次迭代值为x(k)x^(k)x(k),将f(x)f(x)f(x)在x(k)x^(k)x(k)附近进行二阶展开: f(x)=f(x(k))+f′(x(k))Δx+12f′′(x(k)Δx2=f(x(k))+f′(x(k))(x−x(k))+12f′′(x(k))(x−x(k))2 \begin{aligned} f(x) &= f(x^{(k)}) + f'(x^{(k)})\Delta x + \frac{1}{2} f''(x^{(k)}\Delta x^2\\ & = f(x^{(k)}) + f'(x^{(k)})(x - x^{(k)}) + \frac{1}{2} f''(x^{(k)})(x - x^{(k)})^2 \end{aligned} f(x)​=f(x(k))+f′(x(k))Δx+21​f′′(x(k)Δx2=f(x(k))+f′(x(k))(x−x(k))+21​f′′(x(k))(x−x(k))2​ 当这里的xxx是高维时 f(x)=f(x(k))+gkT(x−x(k))+12(x−x(k))TH(x(k))(x−x(k)) f(x) = f(x^{(k)}) +g_k^T(x-x^{(k)}) + \frac{1}{2}(x-x^{(k)})^TH(x^{(k)})(x-x^{(k)}) f(x)=f(x(k))+gkT​(x−x(k))+21​(x−x(k))TH(x(k))(x−x(k)) 其中,gk=g(x(k))=▽f(x(k))g_k = g(x^{(k)}) = \bigtriangledown f(x^{(k)})gk​=g(x(k))=▽f(x(k))是f(x)f(x)f(x)的梯度向量在点x(k)x^{(k)}x(k)的值,H(x(k))H(x^{(k)})H(x(k))是f(x)f(x)f(x)的海赛矩阵在点x(k)x^{(k)}x(k)的值.

海赛矩阵式二阶导数矩阵 H(x)=[∂2f∂xi∂xj]n×n H(x) = \Big[\frac{\partial^2f}{\partial x_i\partial x_j}\Big]_{n\times n} H(x)=[∂xi​∂xj​∂2f​]n×n​

f(x)f(x)f(x)有极值的必要条件是在极值点处一阶导数为0,即梯度向量为0. 当H(x(k))H(x^{(k)})H(x(k))为正定矩阵时,f(x)f(x)f(x)的极值为极小值.

假设▽f(x(k+1))=0\bigtriangledown f(x^{(k+1)}) = 0▽f(x(k+1))=0,我们有 ▽f(x)=f′(x)+f′′(x)(x−x(k)) \bigtriangledown f(x) = f'(x) + f''(x)(x - x^{(k)}) ▽f(x)=f′(x)+f′′(x)(x−x(k)) (1.1)▽f(x)=gk+Hk(x−x(k)) \bigtriangledown f(x) = g_k + H_k(x-x^{(k)})\tag{1.1} ▽f(x)=gk​+Hk​(x−x(k))(1.1) 所以 gk+Hk(x−x(k))=0x(k+1)=x(k)−Hk−1gk g_k + H_k(x-x^{(k)}) = 0\\ x^{(k+1)} = x^{(k)} - H_k^{-1}g_k\\ gk​+Hk​(x−x(k))=0x(k+1)=x(k)−Hk−1​gk​ 我们假设 Hkpk=−gk H_kp_k = -g_k Hk​pk​=−gk​ 则有 x(k+1)=x(k)+pk x^{(k+1)} = x^{(k)} + p_k x(k+1)=x(k)+pk​

算法过程

输入:目标函数f(x)f(x)f(x),梯度g(x)=▽f(x)g(x) = \bigtriangledown f(x)g(x)=▽f(x),海赛矩阵H(x)H(x)H(x),精度要求ϵ\epsilonϵ.

输出:f(x)f(x)f(x)的极小值点x∗x^*x∗.

  1. 取初始点x(0)x^{(0)}x(0),令k=0k=0k=0.
  2. 计算gkg_kgk​.
  3. 若∣∣gk∣∣&lt;ϵ||g_k||&lt;\epsilon∣∣gk​∣∣<ϵ,则停止计算,得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k).
  4. 计算HkH_kHk​,并求pkp_kpk​
  5. x(k+1)=x(k)+pkx^{(k+1)} = x^{(k)} + p_kx(k+1)=x(k)+pk​
  6. k=k+1k=k+1k=k+1,转2

拟牛顿法

基本思想

海赛矩阵的逆矩阵计算困难,考虑用一个n阶矩阵Gk=G(x(k))G_k = G(x^{(k)})Gk​=G(x(k))来近似替代Hk−1H_k^{-1}Hk−1​.

在(1.1)(1.1)(1.1)中代入x=x(k+1)x = x^{(k+1)}x=x(k+1),即有如下: ▽f(x(k+1))=gk+Hk(x(k+1)−x(k))gk+1−gk=Hk(x(k+1)−x(k)) \bigtriangledown f(x^{(k+1)}) = g_k + H_k(x^{(k+1)}-x^{(k)})\\ g_{k+1} - g_k = H_k(x^{(k+1)} - x^{(k)}) ▽f(x(k+1))=gk​+Hk​(x(k+1)−x(k))gk+1​−gk​=Hk​(x(k+1)−x(k))

记yk=g(k+1)−gk,δk=x(x+1)−x(k)y_k = g_{(k+1)} - g_k, \delta_k = x^{(x+1)} - x^{(k)}yk​=g(k+1)​−gk​,δk​=x(x+1)−x(k),则有 (2.1)yk=H)kδkHk−1yk=δk y_k = H)k\delta_k\tag{2.1}\\ H_k^{-1}y_k = \delta_k yk​=H)kδk​Hk−1​yk​=δk​(2.1) 我们称(2.1)(2.1)(2.1)为拟牛顿条件.

如果HkH_kHk​是正定的,那么可以保证牛顿法搜索方向pkp_kpk​是下降方向:

因为搜索方向是pk=−λgkp_k = -\lambda g_kpk​=−λgk​ x=x(k)+λpk=x(k)−λHk(−1)gk x = x^{(k)} + \lambda p_k = x^{(k)} - \lambda H_k^{(-1)}g_k x=x(k)+λpk​=x(k)−λHk(−1)​gk​

所以f(x)f(x)f(x)在x(k)x^{(k)}x(k)的泰勒展开式可以近似写成: f(x)=f(x(k))−λgkTHk−1gk f(x) = f(x^{(k)}) - \lambda g_k^TH_k^{-1}g_k f(x)=f(x(k))−λgkT​Hk−1​gk​ 由于Hk−1H_k^{-1}Hk−1​正定,所以有gkTHk−1gk&gt;0g_k^TH_k^{-1}g_k&gt;0gkT​Hk−1​gk​>0. 当λ\lambdaλ为一个充分小的正数时,总有f(x)&lt;f(x(k))f(x)&lt;f(x^{(k)})f(x)<f(x(k)),也就是说pkp_kpk​是下降方向.

拟牛顿法将GkG_kGk​作为Hk−1H_k^{-1}Hk−1​的近似,要求矩阵GkG_kGk​满足同样的条件,每次迭代矩阵GkG_kGk​都是正定的,且GkG_kGk​要满足拟牛顿条件: Gk1yk=δkG_{k_1}y_k = \delta_kGk1​​yk​=δk​

按照拟牛顿条件选择GkG_kGk​作为Hk−1H_k^{-1}Hk−1​的近似或选择BkB_kBk​作为HkH_kHk​的近似的算法称为拟牛顿法.

按照拟牛顿条件,每次迭代都可以更新矩阵: Gk+1=Gk+ΔGkG_{k+1} = G_k +\Delta G_kGk+1​=Gk​+ΔGk​

有多种具体实现方法

DFP(Davidon-Fletcher-Powell)算法

DFP选择Gk+1G_{k+1}Gk+1​的方法是,假设每一次迭代Gk+1G_{k+1}Gk+1​都是由GkG_kGk​加上两个附加项构成,我们假设这两个附加项分别为Qk、PkQ_k、P_kQk​、Pk​,则有 Gk+1=Gk+Qk+PkGk+1yk=Gkyk+Qkyk+Pkyk G_{k+1} = G_k + Q_k + P_k\\ G_{k+1}y_k = G_ky_k + Q_ky_k + P_ky_k Gk+1​=Gk​+Qk​+Pk​Gk+1​yk​=Gk​yk​+Qk​yk​+Pk​yk​

使Gk+1G_{k+1}Gk+1​满足拟牛顿条件,可以使 Pkyk=δkQkyk=−Gkyk P_ky_k = \delta_k\\ Q_ky_k = - G_ky_k Pk​yk​=δk​Qk​yk​=−Gk​yk​

取 Pk=δkδkTδTykQ=−GkykykTGkykTGkyk P_k = \frac{\delta_k\delta_k^T}{\delta^Ty_k}\\ Q = - \frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} Pk​=δTyk​δk​δkT​​Q=−ykT​Gk​yk​Gk​yk​ykT​Gk​​

所以可以得到Gk+1G_{k+1}Gk+1​的迭代公式 (2.2)Gk+1=Gk−GkykykTGkykTGkyk+δkδkTδTyk G_{k+1} = G_k - \frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} + \frac{\delta_k\delta_k^T}{\delta^Ty_k}\tag{2.2} Gk+1​=Gk​−ykT​Gk​yk​Gk​yk​ykT​Gk​​+δTyk​δk​δkT​​(2.2)

算法过程

输入:f(x)f(x)f(x),梯度g(x)=▽f(x)g(x) = \bigtriangledown f(x)g(x)=▽f(x),精度ϵ\epsilonϵ.

输出:f(x)f(x)f(x)的极小点x∗x^*x∗

  1. 选取初始点x(0)x^{(0)}x(0),取G0G_0G0​为正定对称矩阵,令k=0k=0k=0
  2. 计算gkg_kgk​,若∣∣gk∣∣&lt;ϵ||g_k||&lt;\epsilon∣∣gk​∣∣<ϵ,则停止计算,得到近似解x∗=x(k)x^* = x^{(k)}x∗=x(k)
  3. 令pk=−Gkgkp_k = - G_kg_kpk​=−Gk​gk​
  4. 一维搜索: 求λk\lambda_kλk​使得f(x(k)+λkpk)=min⁡λ≥0f(x(k)+λpk)f(x^{(k)} + \lambda_kp_k) = \min_{\lambda\geq 0}f(x^{(k)}+\lambda p_k)f(x(k)+λk​pk​)=minλ≥0​f(x(k)+λpk​)
  5. 令x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λk​pk​
  6. 计算gk+1g_{k+1}gk+1​,若∣∣gk+1∣∣&lt;ϵ||g_{k+1}||&lt;\epsilon∣∣gk+1​∣∣<ϵ,则停止计算,得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k),否则按式(2.2)(2.2)(2.2)计算Gk+1G_{k+1}Gk+1​.
  7. k=k+1k=k+1k=k+1

BFGS(Boyden-Fletcher-Goldfarb-Shanno)算法

BFGS是最流行得拟牛顿算法.

基本思想: 考虑用BkB_kBk​逼近海赛矩阵HHH,则拟牛顿条件为:

Bk+1δk=ykB_{k+1}\delta_k = y_kBk+1​δk​=yk​ 令 Bk+1=Bk+Pk+QkBk+1δk=Bkδk+Pkδk+Qkδk B_{k+1} = B_k + P_k + Q_k\\ B_{k+1}\delta_k = B_k\delta_k + P_k\delta_k + Q_k\delta_k Bk+1​=Bk​+Pk​+Qk​Bk+1​δk​=Bk​δk​+Pk​δk​+Qk​δk​

使PkP_kPk​和QkQ_kQk​满足: Pkδk=ykQkδk=−Bkδk P_k\delta_k = y_k\\ Q_k\delta_k = -B_k\delta_k Pk​δk​=yk​Qk​δk​=−Bk​δk​ 取 Pk=ykykTykTδkQk=−BkδkδkTBkδkTBkδk P_k = \frac{y_ky_k^T}{y_k^T\delta_k}\\ Q_k = - \frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} Pk​=ykT​δk​yk​ykT​​Qk​=−δkT​Bk​δk​Bk​δk​δkT​Bk​​ 我们得到了Bk+1B_{k+1}Bk+1​的迭代公式: (2.3)Bk+1=Bk−BkδkδkTBkδkTBkδk+ykykTykTδk B_{k+1} = B_k - \frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} + \frac{y_ky_k^T}{y_k^T\delta_k}\tag{2.3} Bk+1​=Bk​−δkT​Bk​δk​Bk​δk​δkT​Bk​​+ykT​δk​yk​ykT​​(2.3) 可以证明,如果初始矩阵B0B_0B0​是正定的,那么迭代过程中的每个矩阵BkB_kBk​都是正定的.

算法过程

输入:f(x)f(x)f(x),梯度g(x)=▽f(x)g(x) = \bigtriangledown f(x)g(x)=▽f(x),精度ϵ\epsilonϵ.

输出:f(x)f(x)f(x)的极小点x∗x^*x∗

  1. 选取初始点x(0)x^{(0)}x(0),取B0B_0B0​为正定对称矩阵,令k=0k=0k=0
  2. 计算gkg_kgk​,若∣∣gk∣∣&lt;ϵ||g_k||&lt;\epsilon∣∣gk​∣∣<ϵ,则停止计算,得到近似解x∗=x(k)x^* = x^{(k)}x∗=x(k)
  3. 令pk=−Bkgkp_k = - B_kg_kpk​=−Bk​gk​
  4. 一维搜索: 求λk\lambda_kλk​使得f(x(k)+λkpk)=min⁡λ≥0f(x(k)+λpk)f(x^{(k)} + \lambda_kp_k) = \min_{\lambda\geq 0}f(x^{(k)}+\lambda p_k)f(x(k)+λk​pk​)=minλ≥0​f(x(k)+λpk​)
  5. 令x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λk​pk​
  6. 计算gk+1g_{k+1}gk+1​,若∣∣gk+1∣∣&lt;ϵ||g_{k+1}||&lt;\epsilon∣∣gk+1​∣∣<ϵ,则停止计算,得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k),否则按式(2.3)(2.3)(2.3)计算Bk+1B_{k+1}Bk+1​.
  7. k=k+1k=k+1k=k+1

Broden类算法

可以从BFGS算法的BkB_kBk​矩阵得到关于GkG_kGk​的迭代公式. 若记Gk=Bk−1、Gk+1=Bk+1−1G_k = B_k^{-1}、G_{k+1} = B_{k+1}^{-1}Gk​=Bk−1​、Gk+1​=Bk+1−1​,对(2.3)(2.3)(2.3)用两次Sherman-Morrison公式可得

Gk+1=(I−δkykTδkTyk)Gk(I−δkykTδkTyk)T+δkykTδkTyk G_{k+1} = \Big(I - \frac{\delta_ky_k^T}{\delta^T_ky_k}\Big)G_k\Big(I - \frac{\delta_ky_k^T}{\delta^T_ky_k}\Big)^T + \frac{\delta_ky_k^T}{\delta^T_ky_k} Gk+1​=(I−δkT​yk​δk​ykT​​)Gk​(I−δkT​yk​δk​ykT​​)T+δkT​yk​δk​ykT​​

参考文献

李航-统计学习方法

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数据结构与算法 - 时间复杂度与空间复杂度

    时间复杂度:时间复杂度的计算并不是计算程序具体运行的时间,而是算法执行语句的最大次数。 空间复杂度:类似于时间复杂度的讨论,一个算法的空间复杂度为该算法所耗费...

    進无尽
  • Core官方DI解析(3)-ServiceCallSite.md

    上一篇说过在整个DI框架中IServiceProviderEngine是核心,但是如果直接看IServiceProviderEngine派生类其实看不出也没什么...

    莫问今朝
  • SSH免密远程登录的配置与实现

    操作系统:CentOS Linux release 7.4.1708 (Core)

    耕耘实录
  • 多达 95% 的 HTTPS 链接能被黑客劫持

    HSTS 是一个被当前大多数 Web 浏览器所支持的 Web 安全策略,它可以帮助 Web 管理员保护他们的服务器和用户避免遭到 HTTPS 降级、中间人攻击和...

    Debian社区
  • 用JWT技术解决IM系统Socket长连接的身份认证痛点1、引言2、原作者3、系列文章5、完全搞懂什么是JWT技术6、我们是怎样使用JWT技术的?7、JWT技术的缺点8、点评附录:更多即时通讯方面的文

    本文引用了封宇《JWT技术解决IM系统的认证痛点》一文的部分内容,即时通讯网重新整理、增补和修订,感谢原作者的无私分享。

    JackJiang
  • 电脑小白学习软件开发(9)-C#基础数组最大值,最小值及排序

    上次说了枚举字符串以及数组的一部分知识点,其实这些东西枯燥的很。小编在以前学习的时候也是如此。虽然枯燥,但这是做所有项目的基础。今天主要讲解点数组的基础知识,这...

    做全栈攻城狮
  • Debian升级内核开启TCP_BBR 实现网络单边加速

    自从锐速发布以来,这款牛逼的单边加速神器的确为一些线路不太优秀的服务器带来了更优秀的体验。但是呢,过高的价格和不再低端售卖。导致了我们并无法实现一个免费好用的单...

    Debian社区
  • React Native APP签名打包release版本APK

    首先React Native开发的APP是无法通过Android Studio进行打包的,因为AS打包的APK,也是和debug版本一样,需要进行依托local...

    做全栈攻城狮
  • Google 宣布新拥堵控制算法 TCP BBR

    Google 宣布了 新拥堵控制算法 TCP BBR。Google 官方博客称新算法将 google.com 和 YouTube 的全球网络吞吐量平均改进了 4...

    Debian社区
  • 巧用OpenSSL完成md2、md4、md5、rmd160、sha、sha1等的验证

    版权声明:本文为耕耘实录原创文章,各大自媒体平台同步更新。欢迎转载,转载请注明出处,谢谢

    耕耘实录

扫码关注云+社区

领取腾讯云代金券