令 X=(x1,x2,⋯,xN)T∈RNX=(x_1,x_2,\cdots,x_N)^T \in {\bf R}^NX=(x1,x2,⋯,xN)T∈RN，目标函数 f:RN→Rf:{\bf R}^N \rightarrow {\bf R}f:RN→R， fff 为凸函数，且二阶连续可微，我们希望求解如下的无约束极小化问题：

min⁡Xf(X)\min_X f(X)Xminf(X)

1 牛顿法

1.1 N=1N=1N=1 时的迭代公式

为了简单起见，这里先考虑 N=1N=1N=1 的情形，此时目标函数 f(X)f(X)f(X) 变为 f(x)f(x)f(x)。

牛顿法的基本思想是： 在现有极小点估计值得附近对 f(x)f(x)f(x) 做二阶泰勒展开，进而找到极小点的下一个估计值。假设 xkx_kxk 是当前的极小点估计值，则：

φ(x)=f(xk)+f′(xk)(x−xk)+12f′′(xk)(x−xk)2\varphi (x) = f(x_k)+f'(x_k)(x-x_k)+\frac{1}{2}f''(x_k)(x-x_k)^2φ(x)=f(xk)+f′(xk)(x−xk)+21f′′(xk)(x−xk)2

表示 f(x)f(x)f(x) 在 xkx_kxk 附近的二阶泰勒展开式（其中略去了关于 x−xkx-x_kx−xk 的高阶项）。因为我们的目标是求最值，由极值的必要条件可知，φ(x)\varphi (x)φ(x) 应该满足：

φ′(x)=f′(xk)+f′′(xk)(x−xk)=0\varphi '(x) = f'(x_k)+f''(x_k)(x-x_k)=0φ′(x)=f′(xk)+f′′(xk)(x−xk)=0

从而有：

x=xk−f′(xk)f′′(xk)x=x_k - \frac{f'(x_k)}{f''(x_k)}x=xk−f′′(xk)f′(xk)

于是，若给定初始值 x0x_0x0，则可以按照下面的迭代公式

xk+1=xk−f′(xk)f′′(xk),k=0,1,⋯x_{k+1}=x_k - \frac{f'(x_k)}{f''(x_k)}, \ k=0,1,\cdotsxk+1=xk−f′′(xk)f′(xk), k=0,1,⋯

产生序列 {xk}\{x_k\}{xk} 来逼近 f(x)f(x)f(x) 的极小值点。

在一定的条件下，{xk}\{x_k\}{xk} 可以收敛到 f(x)f(x)f(x) 的极小值点。

1.2 N>1N>1N>1 时的迭代公式

当 N>1N>1N>1 时，二阶泰勒展示式写作：

φ(X)=f(Xk)+∇f(Xk)⋅(X−Xk)+12⋅(X−Xk)T⋅∇2f(Xk)⋅(X−Xk)\varphi (X) = f(X_k)+\nabla f(X_k)\cdot (X-X_k)+\frac{1}{2} \cdot (X-X_k)^T \cdot \nabla ^2 f(X_k) \cdot (X-X_k)φ(X)=f(Xk)+∇f(Xk)⋅(X−Xk)+21⋅(X−Xk)T⋅∇2f(Xk)⋅(X−Xk)

其中，∇f\nabla f∇f 为 fff 的梯度向量，∇2f\nabla ^2 f∇2f 为 fff 的海森矩阵（Hessian Matrix），其定义分别为：

∇f=[∂f∂x1∂f∂x2⋮∂f∂xN]\nabla f=\left[ \begin{matrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_N} \end{matrix} \right]∇f=⎣⎢⎢⎢⎢⎡∂x1∂f∂x2∂f⋮∂xN∂f⎦⎥⎥⎥⎥⎤

∇2f=[∂2f∂x12∂2f∂x1∂x2⋯∂2f∂x1∂xN∂2f∂x2∂x1∂2f∂x22⋯∂2f∂x2∂xN⋮⋮⋱⋮∂2f∂xN∂x1∂2f∂xN∂x2⋯∂2f∂xN2]N×N\nabla ^2 f=\left[ \begin{matrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_N} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_N} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_N \partial x_1} & \frac{\partial^2 f}{\partial x_N \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_N^2} \end{matrix} \right]_{N \times N}∇2f=⎣⎢⎢⎢⎢⎢⎡∂x12∂2f∂x2∂x1∂2f⋮∂xN∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xN∂x2∂2f⋯⋯⋱⋯∂x1∂xN∂2f∂x2∂xN∂2f⋮∂xN2∂2f⎦⎥⎥⎥⎥⎥⎤N×N

∇f\nabla f∇f 与 ∇2f\nabla ^2 f∇2f 中的函数均为关于 XXX 的函数，分别记为 ggg 和 HHH（ ggg 表示 gradient，HHH 表示 Hessian）。特别地，若 fff 的混合偏导数可交换次序（即对 ∀i,j\forall i,j∀i,j，有 ∂2f∂xi∂xj=∂2f∂xj∂xi\frac{\partial^2 f}{\partial x_i \partial x_j} = \frac{\partial^2 f}{\partial x_j \partial x_i}∂xi∂xj∂2f=∂xj∂xi∂2f），则海森矩阵 HHH 为 对称矩阵，而 ∇f(Xk)\nabla f(X_k)∇f(Xk) 与 ∇2f(Xk)\nabla ^2 f(X_k)∇2f(Xk) 表示将 XXX 取为 XkX_kXk 后得到的实值向量和矩阵，以下分别记为 gkg_kgk 和 HkH_kHk。

同样，因为是求极小点，极值的必要条件要求它是 φ(X)\varphi (X)φ(X) 的驻点，即：

∇φ(X)=gk+Hk⋅(X−Xk)=0\nabla \varphi (X) = g_k + H_k \cdot (X-X_k)=0∇φ(X)=gk+Hk⋅(X−Xk)=0

进一步，若矩阵 HkH_kHk 非奇异，则可解得：

X=Xk−Hk−1⋅gkX=X_k-H_k^{-1} \cdot g_kX=Xk−Hk−1⋅gk

于是，若给定初始值 X0X_0X0，则可以构造出类似的迭代公式：

Xk+1=Xk−Hk−1⋅gk,k=0,1,⋯(1-1)X_{k+1}=X_k-H_k^{-1} \cdot g_k, \ k=0,1,\cdots \tag{1-1}Xk+1=Xk−Hk−1⋅gk, k=0,1,⋯(1-1)

这就是 原始的牛顿迭代法，其迭代公式中的搜索方向 dk=−Hk−1⋅gkd_k=-H_k^{-1} \cdot g_kdk=−Hk−1⋅gk 就称为是 牛顿方向。

1.3 完整的算法描述

算法的伪代码描述为：

给定初始值 X0X_0X0 和精度阈值 ε\varepsilonε，并令 k:=0k:=0k:=0
计算 gkg_kgk 和 HkH_kHk
若 ∥gk∥<ε\|g_k\|<\varepsilon∥gk∥<ε，则停止迭代；否则计算搜索方向 dk=−Hk−1⋅gkd_k=-H_k^{-1} \cdot g_kdk=−Hk−1⋅gk
计算新的迭代点 Xk+1:=Xk+dkX_{k+1} := X_k+d_kXk+1:=Xk+dk
令 k:=k+1k:= k+1k:=k+1，跳到第2步

优点

当目标函数是二次函数时，由于二次泰勒展开函数是与原目标函数完全相同的二次式，海森矩阵退化成一个常数矩阵，从任一初始点出发，利用公式（1-1）只需一步迭代就可以达到 f(x)f(x)f(x) 的极小点 x∗x^*x∗，因此牛顿法是一种具有 二次收敛性 的算法。

对于非二次函数，若函数的二次性态较强，或迭代点已进入极小点的邻域，则其收敛速度也是很快的，这是牛顿法的主要优点。

缺点

由于原始牛顿法的迭代公式中没有步长因子，而是定长迭代，对于非二次型目标函数，有时会使函数值上升，即出现 f(Xk+1)>f(Xk)f(X_{k+1}) > f(X_k)f(Xk+1)>f(Xk) 的情况，这表明 原始牛顿法不能保证函数值稳定的下降，在严重的情况下甚至可能造成迭代点列 {Xk}\{X_k\}{Xk} 的发散而导致计算失败。

2 阻尼牛顿法

为了消除原始牛顿法的缺点，人们提出了阻尼牛顿法。

对于牛顿法，确定了迭代方向之后，迭代步长默认为1，但是这个迭代方向并不一定是朝着函数值下降的方向。阻尼牛顿法每次的迭代方向仍采用 dkd_kdk，但每次迭代需沿此方向做一维搜索（line search），寻找最优的步长因子 λk\lambda_kλk，即：

λk=argmin⁡λ∈Rf(Xk+λdk)(2-1)\lambda_k=arg \min_{\lambda \in R}f(X_k+\lambda d_k) \tag{2-1}λk=argλ∈Rminf(Xk+λdk)(2-1)

算法的伪代码描述为：

给定初始值 X0X_0X0 和精度阈值 ε\varepsilonε，并令 k:=0k:=0k:=0
计算 gkg_kgk 和 HkH_kHk
若 ∥gk∥<ε\|g_k\|<\varepsilon∥gk∥<ε，则停止迭代；否则计算搜索方向 dk=−Hk−1⋅gkd_k=-H_k^{-1} \cdot g_kdk=−Hk−1⋅gk
利用公式（2-1）得到步长λk\lambda_kλk，并计算新的迭代点 Xk+1:=Xk+λkdkX_{k+1} := X_k+\lambda_k d_kXk+1:=Xk+λkdk
令 k:=k+1k:= k+1k:=k+1，跳到第2步

可以看到，阻尼牛顿法相比于牛顿法，在每次参数更新之前，利用一维搜索法计算更新步长，确保优化方向为下降方向。

3 拟牛顿法

原始牛顿法虽然收敛速度快，但是需要计算海森矩阵的逆矩阵 H−1H^{-1}H−1，而且有时目标函数的海森矩阵无法保持正定，从而使得原始牛顿法失效。为了克服这两个问题，人们提出了拟牛顿法。这个方法的 基本思想 是：不用二阶偏导数，而是构造出一个可以近似海森矩阵（或海森矩阵的逆）的正定对称阵。不同的构造方法就产生了不同的拟牛顿法。

下面我们先推导一下拟牛顿条件，它给“对海森矩阵（或海森矩阵的逆）做近似”提供了理论指导，指出了用来近似的矩阵应该满足的条件。

3.1 拟牛顿条件

对 ∇f(x)\nabla f(x)∇f(x) 在 xkx_kxk 做泰勒展开我们可以得到以下近似：

∇f(x)=gk+Hk(x−xk)\nabla f(x)=g_k+H_k(x-x_k)∇f(x)=gk+Hk(x−xk)

取 x=xk+1x=x_{k+1}x=xk+1，则有：

∇f(xk+1)=gk+1=gk+Hk(xk+1−xk)\nabla f(x_{k+1})=g_{k+1} = g_k+H_k(x_{k+1}-x_k)∇f(xk+1)=gk+1=gk+Hk(xk+1−xk)

即：

gk+1−gk=Hk(xk+1−xk)g_{k+1} - g_k=H_k(x_{k+1}-x_k)gk+1−gk=Hk(xk+1−xk)

记 yk=gk+1−gky_k=g_{k+1}-g_kyk=gk+1−gk，δk=xk+1−xk\delta_k=x_{k+1}-x_kδk=xk+1−xk，则有：

yk=Hkδk(3-1)y_k=H_k \delta_k \tag{3-1}yk=Hkδk(3-1)

或

Hk−1yk=δk(3-2)H_k^{-1}y_k=\delta_k \tag{3-2}Hk−1yk=δk(3-2)

以上即为拟牛顿条件。

常用的拟牛顿法有DFP、BFGS、L-BFGS，区别在于如何选取替代矩阵。

3.2 DFP算法

GFP算法用于近似拟牛顿条件（3-2），这里用 GkG_kGk 代表对 Hk−1H_k^{-1}Hk−1 的近似，下面直接给出计算公式：

Gk+1=Gk+δkδkTδkTyk−GkykykTGkykTGkykG_{k+1}=G_k+\frac{\delta_k \delta_k^T}{\delta_k^T y_k} - \frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k}Gk+1=Gk+δkTykδkδkT−ykTGkykGkykykTGk

可以证明，如果初始矩阵 G0G_0G0 是正定对称的，则迭代过程中的每个矩阵 GkG_kGk 都是正定对称的，一般取 G0=IG_0=IG0=I。

3.3 BFGS算法

BFGS 算法用于近似拟牛顿条件（3-1），与 DFP 相比，BFGS 的性能更佳。这里用 BkB_kBk 代表对 HkH_kHk 的近似，下面直接给出计算公式：

Bk+1=Bk+ykykTykTδk−BkδkδkTBkδkTBkδkB_{k+1}=B_k+\frac{y_k y_k^T}{y_k^T \delta_k} - \frac{B_k \delta_k \delta_k^T B_k}{\delta_k^T B_k \delta_k}Bk+1=Bk+ykTδkykykT−δkTBkδkBkδkδkTBk

可以证明，如果初始矩阵 B0B_0B0 是正定对称的，则迭代过程中的每个矩阵 BkB_kBk 都是正定对称的，一般取 B0=IB_0=IB0=I。

若记 Gk=Bk−1G_k=B_k^{-1}Gk=Bk−1， Gk+1=Bk+1−1G_{k+1}=B_{k+1}^{-1}Gk+1=Bk+1−1 ，那么应用Sherman-Morrison公式可以将上述迭代公式改写为：

Gk+1=(I−δkykTδkTyk)Gk(I−δkykTδkTyk)T+δkδkTδkTykG_{k+1}=(I-\frac{\delta_k y_k^T}{\delta_k^T y_k}) G_k (I-\frac{\delta_k y_k^T}{\delta_k^T y_k})^T + \frac{\delta_k \delta_k^T}{\delta_k^T y_k}Gk+1=(I−δkTykδkykT)Gk(I−δkTykδkykT)T+δkTykδkδkT

上式就是 BFGS 算法关于 GkG_kGk 的迭代公式。

3.4 L-BFGS算法

在BFGS中，需要用到一个 NNN 阶矩阵 GkG_kGk ，当 NNN 很大时，存储这个矩阵将消耗大量的计算机资源。为了解决这个问题，减少 BFGS 迭代过程中所需的内存开销，就有了 L-BFGS。

L-BFGS（Limited-memory BFGS 或 Limited-storage BFGS）对 BFGS 进行了近似，其 基本思想 是：不再存储完整的矩阵 GkG_kGk ，而是存储计算过程中的向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk}{yk}，需要矩阵 GkG_kGk 时，利用向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk}{yk} 的计算来代替。而且，向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk}{yk} 也不是所有的都存储，而是保留最新的 mmm 个，每次计算 GkG_kGk 时，只利用最新的 mmm 个向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk}{yk}。这样一来，存储空间由 O(N2)O(N^2)O(N2) 降至 O(mN)O(mN)O(mN)。

牛顿法（Newton Methods）、阻尼牛顿法和拟牛顿法相关推荐

什么是牛顿法(Newton methods)？什么是拟牛顿法(Quasi Newton methods)？牛顿法和梯度下降法的区别是什么？
什么是牛顿法(Newton methods)?什么是拟牛顿法(Quasi Newton methods)?牛顿法和梯度下降法的区别是什么? 牛顿法的最初提出是用来求解方程的根的.对于最优化问题,其极值 ...
《统计学习方法》啃书辅助：附录B 牛顿法和拟牛顿法
<统计学习方法>啃书辅助:附录B 牛顿法和拟牛顿法 B.1 牛顿法 [基础知识]梯度 [基础知识]浅谈「正定矩阵」和「半正定矩阵」 - Xinyu Chen 的文章 - 知乎二阶泰勒展开 ...
机器学习优化算法中梯度下降,牛顿法和拟牛顿法的优缺点详细介绍
1.梯度下降法梯度下降法实现简单,当目标函数是凸函数时,梯度下降法的解是全局解.一般情况下,其解不保证是全局最优解,梯度下降法的速度也未必是最快的. 梯度下降法的优化思想:用当前位置负梯度方向作为搜 ...
梯度类算法原理：最速下降法、牛顿法和拟牛顿法
文章目录算法结构最速下降法牛顿法拟牛顿法算法结构梯度类算法有很多,本文主要学习最常见的3个算法:最速下降法.牛顿法和拟牛顿法.算法名称虽多,但是他们的算法结构都是一样的,可以描述为 (1) ...
二阶偏微分方程组龙格库塔法_牛顿法和拟牛顿法——（书中附录B）
牛顿法(Newton method)和拟牛顿法(quasi-Newton method)也是求解无约束最优化问题的常用方法,具有收敛速度快的优点. 牛顿法是迭代算法,每一步需要求解目标函数的海赛矩阵的 ...
梯度下降法、牛顿法和拟牛顿法——机器学习面试
梯度下降.牛顿.拟牛顿原理梯度下降牛顿法为Hesse矩阵参数更新的方程: 每一次迭代的更新方向都是当前点的牛顿方向,步长固定为1.每一次都需要计算一阶导数以及Hesse矩阵的逆矩阵,对于 ...
统计学习方法--牛顿法和拟牛顿法
与公众号同步更新,详细内容及相关ipynb文件在公众号中,公众号:AI入门小白文章目录牛顿法拟牛顿法的思路 DFP (Davidon-Fletcher- Powell) 算法(DFP algor ...
牛顿法(Newton‘s method)和拟牛顿法(quasi Newton method)
简述在看伊恩·古德费洛的深度学习,4.3节基于梯度的优化方法时提到仅使用梯度信息的优化算法称为一阶优化算法 ,如梯度下降. 使用Hessian矩阵的优化算法称为二阶最优化算法 ,如牛顿法. 牛 ...
八、梯度下降法和拟牛顿法
1.梯度 2.梯度上升和梯度下降 3.梯度下降算法详解 3.1 直观解释 3.2 梯度下降相关概念 3.3 梯度下降的矩阵描述 3.4 梯度下降的算法调优 4.梯度下降法大家族 5.梯度下降法和其他无 ...

牛顿法（Newton Methods）、阻尼牛顿法和拟牛顿法