前言

无论是牛顿法，高斯牛顿法，还是LM算法，都需要计算海森矩阵或其近似，然后求海森矩阵的逆，来获得迭代增量，因此存在迭代不收敛问题（海森矩阵的正定性），并且计算效率较低。

拟牛顿法更进一步，将海森矩阵也作为一个迭代估计量，因而避免了海森矩阵的求逆运算，并提升了迭代稳定性。

拟牛顿条件

仍然是从优化函数的泰勒展开讲起：
f(x)=f(xk)+∇f(xk)T(x−xk)+12(x−xk)TH(xk)(x−xk)∇f(x)T=∇f(xk)T+H(xk)(x−xk)f({\bf x})=f({\bf x_k})+\nabla f({\bf x_k})^T({\bf x - x_k}) + \frac{1}{2}({\bf x - x_k})^T H({\bf x_k})({\bf x - x_k}) \\ \nabla f ({\bf x})^T = \nabla f({\bf x_k})^T + H({\bf x_k})({\bf x-x_k}) f(x)=f(xk)+∇f(xk)T(x−xk)+21(x−xk)TH(xk)(x−xk)∇f(x)T=∇f(xk)T+H(xk)(x−xk)
简写为：
J(x)−J(xk)=H(xk)(x−xk)J({\bf x}) - J({\bf x_k}) = H({\bf x_k})({\bf x-x_k}) J(x)−J(xk)=H(xk)(x−xk)
假设我们通过J(x)=0J(\bf x)=0J(x)=0求出了本次迭代增量Δx\Delta \bf xΔx，现在将x=xk+1=xk+Δx{\bf x=x_{k+1}=x_k}+\Delta {\bf x}x=xk+1=xk+Δx代入：
J(xk+1)−J(xk)=H(xk)(xk+1−xk)(1−1)J({\bf x_{k+1}})-J({\bf x_k}) = H({\bf x_{k}}) ({\bf x_{k+1} - x_k}) \quad (1-1) J(xk+1)−J(xk)=H(xk)(xk+1−xk)(1−1)
如果我们要通过迭代Gk+1G_{k+1}Gk+1近似H(xk)−1H({x_k})^{-1}H(xk)−1，那么GkG_kGk也要满足以上方程。这就是拟牛顿条件：
Gk+1(Jk+1−Jk)=(xk+1−xk)setJk+1−Jk=yk,(xk+1−xk)=skthenGk+1yk=sk(2)andyk=Gk+1−1sk(3)G_{k+1}(J_{k+1}-J_{k})=({\bf x_{k+1} - x_k}) \\ set \quad J_{k+1}-J_{k}=y_k,({\bf x_{k+1} - x_k})=s_k \\ \quad \\ then \quad G_{k+1}y_k=s_k \quad (2) \\ and \quad y_k = G_{k+1}^{-1}s_k \quad (3) \\ Gk+1(Jk+1−Jk)=(xk+1−xk)setJk+1−Jk=yk,(xk+1−xk)=skthenGk+1yk=sk(2)andyk=Gk+1−1sk(3)
注意：用xk+1,Jk+1{\bf x_{k+1}},J_{k+1}xk+1,Jk+1迭代Gk+1G_{k+1}Gk+1似乎只能近似第kkk次迭代的海森矩阵H(xk)H(\bf x_k)H(xk)。不过实际上，如果把式(1-1)改写为以下形式：
J(xk−1)−J(xk)=H(xk)(xk−1−xk)(1−2)→J(xk)−J(xk−1)=H(xk)(xk−xk−1)→J(xk+1)−J(xk)=H(xk+1)(xk+1−xk)J({\bf x_{k-1}})-J({\bf x_k}) = H({\bf x_{k}}) ({\bf x_{k-1} - x_k}) \quad (1-2) \\ \to J({\bf x_{k}})-J({\bf x_{k-1}}) = H({\bf x_{k}}) ({\bf x_{k} - x_{k-1}}) \\ \to J({\bf x_{k+1}})-J({\bf x_{k}}) = H({\bf x_{k+1}}) ({\bf x_{k+1} - x_{k}}) \\ J(xk−1)−J(xk)=H(xk)(xk−1−xk)(1−2)→J(xk)−J(xk−1)=H(xk)(xk−xk−1)→J(xk+1)−J(xk)=H(xk+1)(xk+1−xk)
这样，Gk+1G_{k+1}Gk+1就是H(xk+1)H(\bf x_{k+1})H(xk+1)的近似了。

几何解释

用一元函数来理解更简单，如下图，黑色线是优化函数的一阶导，A,BA,BA,B点是xk,xk+1x_k,x_{k+1}xk,xk+1，拟牛顿法的思想就是通过红色割线ABABAB的斜率来近似BBB点的二阶导（也就是B点的切线，绿线）；另一方面，红色切线对于AAA点的切线也能近似，这就是式(1-1),(1-2)的几何意义。

DFP算法

DFP算法（DFP是人名）通过迭代构造GkG_{k}Gk来近似Hk−1H_{k}^{-1}Hk−1。DFP公式推导如下：
Gk+1=Gk+Pk+QkGk+1yk=sk→Gkyk+Pkyk+Qkyk=sksetPkyk=sk,Qkyk=−GkykthenPk=skskTskTyk,Qk=−GkykykTGkykTGkykGk+1=Gk+skskTskTyk−GkykykTGkykTGkykG_{k+1}=G_k+P_k+Q_k \\ G_{k+1}y_k=s_k \to G_{k}y_k+P_ky_k+Q_ky_k=s_k \\ \quad \\ set \quad P_ky_k=s_k,\quad Q_ky_k=-G_ky_k \\ \quad \\ then \quad P_k= \frac{s_ks_k^T}{s_k^Ty_k},\quad Q_k=-\frac {G_ky_ky_k^TG_k}{y_k^TG_ky_k} \\ \quad \\ G_{k+1} = G_k+ \frac{s_ks_k^T}{s_k^Ty_k} - \frac {G_ky_ky_k^TG_k}{y_k^TG_ky_k} \\ Gk+1=Gk+Pk+QkGk+1yk=sk→Gkyk+Pkyk+Qkyk=sksetPkyk=sk,Qkyk=−GkykthenPk=skTykskskT,Qk=−ykTGkykGkykykTGkGk+1=Gk+skTykskskT−ykTGkykGkykykTGk

DFP算法中，当初始GGG正定时，迭代中的每个GGG都是正定的。

BFGS算法

BFGS算法（BFGS是四位作者的首字母）通过构造BkB_kBk来近似HkH_kHk，推导方式与DFP算法相似：
Bk+1=Bk+Pk+QkBk+1sk=yk→Bksk+Pksk+Qksk=yksetPksk=yk,Qksk=−BkskthenPk=ykykTykTsk,Qk=−BkskskTBkskTBkskBk+1=Bk+ykykTykTsk−BkskskTBkskTBkskB_{k+1}=B_k+P_k+Q_k \\ B_{k+1}s_k=y_k \to B_{k}s_k+P_ks_k+Q_ks_k=y_k \\ \quad \\ set \quad P_ks_k=y_k,\quad Q_ks_k=-B_ks_k \\ \quad \\ then \quad P_k= \frac{ y_k y_k^T}{ y_k^Ts_k},\quad Q_k=-\frac {B_ks_ks_k^TB_k}{s_k^TB_ks_k} \\ \quad \\ B_{k+1} = B_k+ \frac{ y_k y_k^T}{ y_k^Ts_k} - \frac {B_ks_ks_k^TB_k}{s_k^TB_ks_k} \\ Bk+1=Bk+Pk+QkBk+1sk=yk→Bksk+Pksk+Qksk=yksetPksk=yk,Qksk=−BkskthenPk=ykTskykykT,Qk=−skTBkskBkskskTBkBk+1=Bk+ykTskykykT−skTBkskBkskskTBk
BFGS算法中，当初始BBB正定时，迭代中的每个BBB都是正定的。

注意：由于原BFGS算法没有直接近似海森矩阵的逆，因此可以进一步应用Sherman-Morrison-Woodbury公式直接从Bk−1B_k^{-1}Bk−1得到Bk+1−1B_{k+1}^{-1}Bk+1−1。由于这个公式我还不大熟，因此这里给出结果，以后有时间再记录推导细节：
Bk+1−1=(I−skykTykTsk)Bk−1(I−ykskTykTsk)+skskTykTskB_{k+1}^{-1} = (I - \frac{s_ky_k^T}{y_k^Ts_k})B_k^{-1}(I-\frac{y_ks_k^T}{y^T_ks_k}) + \frac{s_ks_k^T}{y_k^Ts_k} Bk+1−1=(I−ykTskskykT)Bk−1(I−ykTskykskT)+ykTskskskT

Broyden类算法

上面的DFP和BFGS都能得到Hk−1H_k^{-1}Hk−1的迭代近似Gk−DFP,Gk−BFGSG_{k-DFP},G_{k-BFGS}Gk−DFP,Gk−BFGS，那么DFP和BFGS迭代的线性组合也满足拟牛顿条件，并且迭代保持正定：
Gk=λGk−DFP+(1−λ)Gk−BFGS,0≤λ≤1G_{k}=\lambda G_{k-DFP}+(1-\lambda)G_{k-BFGS},\quad 0\le\lambda\le 1 Gk=λGk−DFP+(1−λ)Gk−BFGS,0≤λ≤1
这被称为Broyden类算法。

代码示例

下面是我写的通过拟牛顿法DFP求二元函数极小值的python代码：

import numpy as np
import scipy.optimize
import time
import mathdef 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 non_linear_func_5(x, y):fxy = math.pow(math.exp(x) + math.exp(0.5 * y) + x, 2)return fxydef quasi_newton_optim(x, y, F, G, lr=0.1):dx, dy = partial_derivate_xy(x, y, F)grad = np.array([[dx], [dy]])vec_delta = - lr * np.matmul(G, grad)vec_opt = np.array([[x], [y]]) + vec_deltax_opt = vec_opt[0][0]y_opt = vec_opt[1][0]dxx, dyy = partial_derivate_xy(x_opt, y_opt, F)grad_delta = np.array([[dxx - dx], [dyy - dy]])G_update = G + np.matmul(vec_delta, vec_delta.T)/np.dot(vec_delta.T, grad_delta) \- np.matmul(np.matmul(np.matmul(G, grad_delta), grad_delta.T), G) / np.matmul(np.matmul(grad_delta.T, G), grad_delta)return x_opt, y_opt, grad, G_updatedef quasi_newton_optim_correct(x, y, F, G):dx, dy = partial_derivate_xy(x, y, F)grad = np.array([[dx], [dy]])vec_delta = - np.matmul(G, grad)vec_opt = np.array([[x], [y]]) + vec_deltax_opt = vec_opt[0][0]y_opt = vec_opt[1][0]dxx, dyy = partial_derivate_xy(x_opt, y_opt, F)grad_delta = np.array([[dxx - dx], [dyy - dy]])Gy = np.matmul(G, grad_delta)yG = np.matmul(grad_delta.T, G)yGy = np.matmul(np.matmul(grad_delta.T, G), grad_delta)G_update = G + np.matmul(vec_delta, vec_delta.T)/np.dot(vec_delta.T, grad_delta) \- np.matmul(Gy, yG) / yGyreturn x_opt, y_opt, grad, G_updatedef optimizer(x0, y0, F, th=0.01):x = x0y = y0G = np.eye(2)counter = 0while True:x_opt, y_opt, grad, G = quasi_newton_optim(x, y, F, G)if np.linalg.norm(grad) < th:breakx = x_opty = y_optcounter = counter + 1print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y))return x, ydef verify_min(x, y, F):fx = F(x, y)deltax = np.linspace(-0.1, 0.1, 100)deltay = np.linspace(-0.1, 0.1, 100)x_range = x + deltaxy_range = y + deltaycounter = 0for i in range(100):for j in range(100):f_range = F(x_range[i], y_range[j])f_delta = fx - f_rangeif f_delta < 0:counter += 1print('counter: {}'.format(counter))def scipyF(X):x, y = Xfxy = x ** 4 + 2 * y ** 4 + 3 * x ** 2 * y ** 2 + 4 * x * y ** 2 + x * y + x + 2 * y + 0.5return fxyif __name__ == '__main__':x0 = 2.y0 = 2.start = time.time()for i in range(1000):result_x, result_y = optimizer(x0, y0, non_linear_func_5)if i == 0:breakend = time.time()print(result_x, result_y, 'cost time: {}'.format(end - start))print(partial_derivate_xy(result_x, result_y, non_linear_func_5))verify_min(result_x, result_y, non_linear_func_5)# scipyRes = scipy.optimize.fmin_cg(scipyF, np.array([0, 0]))# print(scipyRes)

后记

牛顿法与拟牛顿法只剩一个LBFGS算法了。为了加深印象和实际应用，下一篇LBFGS我会把这一块的代码做成一个类来使用，并且添加可视化结果。

数值计算之拟牛顿法相关推荐

数值计算之线搜索法，Armijo，Wolfe，Goldstein条件，回溯法
数值计算之线搜索法,Wolfe conditions 前言梯度法的步长线搜索法非精确线搜索法 Armijo条件 Wolfe条件 Goldstein条件回溯法后记前言本篇是基于梯度的优化 ...
数值计算之 LBFGS
数值计算之 LBFGS 前言拟牛顿法 LBFGS算法代码示例后记前言本篇是非线性优化方法中的最后一个部分,LBFGS算法. 拟牛顿法上篇记录了拟牛顿法的思路:通过迭代矩阵GkG_{k}Gk ...
数值计算之牛顿法与函数极值
数值计算之牛顿法与函数极值前言最速下降法牛顿法牛顿法分析代码示例后记前言本篇继续优化理论的算法学习,牛顿法. 最速下降法首先回顾上次提到的梯度下降法(其实就是最速下降法):通过求取 ...
深度学习之数学基础（数值计算）
信息论是应用数学的一个分支,主要研究的是对一个信号能够提供信息的多少进行量化.如果说概率使我们能够做出不确定性的陈述以及在不确定性存在的情况下进行推理,那信息论就是使我们能够量化概率分布中不确定性的总 ...
java数值计算算法编程,Java数值计算算法编程
第1章 Java与数值计算. 1．1 数值计算中存在的问题 1．2 用Java实现数值计算算法的要点 1．3 实数类设计与实现第2章复数运算 2．1 复数类设计 2．2 复数乘法 2．3 复数除法 ...
计算机数值计算原理,C#数值计算算法编程
第1章 C#与数值计算. 1．1 数值计算中存在的问题 1．2 用C#实现数值计算算法的要点第2章复数运算 2．1 复数类设计 2．2 复数乘法 2．3 复数除法 2．4 复数的模 2．5 复数的 ...
【opencv】(2) 图像处理：边界填充、图像融合、图像阈值、数值计算
主要内容有:边界填充 cv2.copyMakeBorder(),数值计算 cv2.add(),改变尺寸 cv2.resize(),图像融合 cv2.addWeighted(),图像阈值 cv2.thr ...
什么是牛顿法(Newton methods)？什么是拟牛顿法(Quasi Newton methods)？牛顿法和梯度下降法的区别是什么？
什么是牛顿法(Newton methods)?什么是拟牛顿法(Quasi Newton methods)?牛顿法和梯度下降法的区别是什么? 牛顿法的最初提出是用来求解方程的根的.对于最优化问题,其极值 ...
梯度下降之模拟退火、梯度下降之学习计划、牛顿法、拟牛顿法、共轭梯度法
梯度下降之模拟退火.梯度下降之学习计划.牛顿法.拟牛顿法.共轭梯度法目录

数值计算之拟牛顿法

数值计算之拟牛顿法

前言

拟牛顿条件

几何解释

DFP算法

BFGS算法

Broyden类算法

代码示例

后记

数值计算之拟牛顿法相关推荐

最新文章

热门文章

数值计算之 拟牛顿法

数值计算之 拟牛顿法

前言

拟牛顿条件

几何解释

DFP算法

BFGS算法

Broyden类算法

代码示例

后记

数值计算之 拟牛顿法相关推荐

最新文章

热门文章

数值计算之拟牛顿法

数值计算之拟牛顿法

数值计算之拟牛顿法相关推荐