优化理论之牛顿法

Nov 4, 2019 · 优化理论 ·

分享到:

牛顿法

基础牛顿法
牛顿法求解函数根
牛顿法求根的几何解释
函数最优值与求根
- 一元函数情况
- 多元函数的情况
实现算法
牛顿法修正
Goldstein-Price修正法(G-P法)
负曲率方向法
Goldfeld修正
Gill-Murray改进Cholesky分解法的稳定牛顿法
信赖域法
拟牛顿法
DFP算法
SR1算法
BFGS算法
L-BFGS算法
共轭梯度法
总结
参考文献

牛顿法和梯度下降法都是求解非线性无约束最优化问题的常用方法，但是效果速度比梯度法更好（但是计算量增加了，需要二阶信息）。牛顿法是迭代算法，每一步需要求解目标函数的 Hessian 矩阵的逆矩阵，计算过程复杂。虽然牛顿法收敛速度很快，但是有两个很难处理的问题。一是要保证Hessian 矩阵正定。二是求Hessian 矩阵是很耗计算量的。因此，有许多修正牛顿法来保证矩阵正定或是在矩阵不正定时提供候选方向。同时为了简化计算，诞生了拟牛顿法通过正定矩阵近似 Hessian 矩阵的逆矩阵或 Hessian 矩阵，简化计算过程。DFP方法、BFGS方法和L-BFGS方法都属于拟牛顿法。

牛顿法、修正的牛顿法和拟牛顿法解决的问题可以用如下无约束最小化问题概括： \[\mathop{min}\limits_{x}\quad f(x)\] 其中，\(x=(x_1，x_2，\ldots，x_N)^T\in\mathbb{R}^N\)。由于本文不讨论收敛性，没有特殊说明时，假设\(f\)为凸函数，且两阶连续可微，此外，记极小问题的解是\(x^\ast\)。

基础牛顿法

牛顿法求解函数根

牛顿法最核心的功能是迭代求函数根，其他函数优化功能都是从这里衍生出去的。我们知道并不是所有的方程都有求根公式，或者求根公式很复杂，导致求解困难。利用牛顿法，可以迭代求解。原理是利用泰勒公式。我们将任一函数在\(x_0\)处展开，且展开到一阶，即： \[f(x)=f(x_0)+(x–x_0)f'(x_0)+\underbrace{\frac{1}{2}(x-x_0)^2f''(\varepsilon)}_{余项}\] 其中，\(\varepsilon\)处于\(x,x_0\)之间，由于后面的余项是\((x-x_0)\)的高阶数，先忽略。因此要求解方程\(f(x)=0\)，即是\(f(x_0)+(x–x_0)f'(x_0)=0\)，解为： \[x=x_1=x_0–f(x_0)/f'(x_0)\] 因为这是利用泰勒公式的一阶展开，且忽略了余项，余项的存在导致\(f(x)与f(x_0)+(x–x_0)f'(x_0)\)在\(x_1\)并不是完全相等，而是近似相等。因此，这里求得的\(x_1\)并不能让\(f(x)=0\)，只能说\(f(x_1)\)的值比\(f(x_0)\)更接近\(f(x)=0\)，于是乎，迭代求解的想法就很自然了，可以进而推出 \[x_{n+1}=x_n–f(x_n)/f'(x_n)\] 通过迭代，这个式子必然在\(f(x^∗)=0\)的时候收敛. 整个过程如下图：

图1牛顿法求根

牛顿法求根的几何解释

牛顿法求根其实是利用了切线的性质。如果我们查看过点\(x_0\)处的切线方程： \[ y-f(x_0)=f'(x_0)(x-x_0)\\ \] 当\(y=0\)时，有： \[ -f(x_0)=f'(x_0)(x-x_0)\\ \Rightarrow x=x_0-\frac{f(x_0)}{f'(x_0)} \] 这个式子和牛顿法求根的迭代公式一致。说白了，牛顿法求根就是利用当前点所在直线与\(y=0\)的交点，这个交点\(x_1\)会比\(x_0\)更接近让函数值为0。

函数最优值与求根

一元函数情况

从简单的一元函数开始，我们先列出函数的二阶泰勒展开式： \[f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{1}{2}f''(x_0)(x-x_0)^2+R_2(x)\]

然后，我们再来讨论我们的问题，求函数最优值。我们知道一个可微函数的最值点，一定是它的驻点（或边界点，假设定义域范围暂不考虑）。那我们就要求驻点，即\(f'(x)=0\)的点。我们如果对\(f(x)\)和其二阶泰勒展开式求导（忽略余项）： \[f'(x)=f'(x_0)+f''(x_0)(x-x_0)\] 根据驻点性质，要求\(f'(x)=0\)，则根据上一节牛顿法求解函数根可知： \[x_1=x_0-\frac{f'(x_0)}{f''(x_0)}\] 这样我们就得到了下一点的位置\(x_1\)。\(f'(x_1)\)比\(f'(x_0)\)更接近0。接下来重复这个过程，直到到达导数为0的点，由此得到求极值的牛顿法的迭代公式： \[x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}\] 给定初始迭代点\(x_0\)，反复用上面的公式进行迭代，直到达到导数为0的点或者达到最大迭代次数。

多元函数的情况

下面推广到多元函数的情况，根据多元函数的泰勒展开公式（参考矩阵论教材），我们对目标函数在\(x_0\)点处做二阶泰勒展开，有（自变量皆为列向量，忽略余项）： \[f(x)=f(x_0)+\nabla f(x_0)^T(x-x_0)+\frac{1}{2}(x-x_0)^T\nabla^2f(x_0)(x-x_0)\] 其中，\(\nabla^2f(x_0)\)为多元实值函数\(f(x)\)在\(x_0\)处的海森矩阵（Hessian Matrix），我们以后用\(H\)表示。那么如果函数梯度为0，且海森矩阵\(H\)可逆，则有： \[\nabla f(x)=\vec{0} \Rightarrow \\ \nabla f(x_0)+\nabla^2f(x_0)(x-x_0)=\vec{0}\Rightarrow \\ x_1=x_0-(\nabla^2f(x_0))^{-1}\nabla f(x_0)\Rightarrow \\ x_1=x_0-H^{-1}[f(x_0)]\nabla f(x_0)\] 注意矩阵左右乘不可颠倒,\(H^{-1}[f(x_0)]\)，为海森矩阵的逆矩阵。特别地，若是极小值点则Hessian是正定矩阵。同理，我们可以得到如下递推公式： \[x_{n+1}=x_n-H^{-1}[f(x_n)]\nabla f(x_n) \\ 可简写为：x_{n+1}=x_n-H^{-1}_n{g_n}\] 最终会到达函数的驻点处。其中\(-H^{-1}g\)称为牛顿方向。迭代终止的条件是梯度的模接近于0，或者函数值下降小于指定阈值。

注意，我们目前默认牛顿法的步长是1，这样可以保证2阶收敛速度，但是并不是对所有函数都有效。现不加证明的指出，带步长的牛顿法是收敛的。

实现算法

1.给定初始值\(x_0\)和精度阈值\(\epsilon\)，设置\(k = 0\)

2.计算梯度\(g_k\)和矩阵\(H_k\)

3.如果\(\|g_k\|<\epsilon\)即在此点处梯度的值接近于0，则达到极值点处，停止迭代

4.计算搜索方向\(d_k=-H^{-1}_k{g_k}\)

5.计算新的迭代点\(x_{k+1}=x_k-\alpha H^{-1}_k{g_k}\)(一般可以取\(\alpha=1\))

6.令\(k = k + 1\)，返回步骤2

其中\(\epsilon\)是一个人工设定的接近于0的常数，和梯度下降法一样，需要这个参数的原因是保证\(x_{k+1}\)在\(x_k\)的邻域内，从而可以忽略泰勒展开的高次项。如果目标函数是二次函数，Hessian矩阵是一个常数矩阵，对于任意给定的初始点，牛顿法只需要一步迭代就可以收敛到极值点。

牛顿法收敛速度快，具有二阶收敛性。但是和梯度下降法一样，牛顿法寻找的也是导数为0的点，这不一定是极值点，因此会面临局部极小值和鞍点问题。另外，它需要Hessian矩阵，计算量大，并且Hessian矩阵并不是一直能够求出。因此出现了一些降低计算量得拟牛顿法。

除此之外，牛顿法在每次迭代时序列\(x_i\)可能不会收敛到一个最优解，它甚至不能保证函数值会按照这个序列递减。解决第一个问题可以通过调整牛顿方向的步长来实现，目前常用的方法有两种：直线搜索和可信区域法。加了步长搜索的牛顿法也叫阻尼牛顿法。

牛顿法修正

牛顿法的一大困难在于某一k步的Hessian 矩阵\(H_k\)可能不正定。这时二次型模型不一定有极小点，甚至没有平稳点。此时，牛顿方向\(-H_k^{-1}g_k\)不一定会让目标函数值下降。对应的各种改进方案应运而生，各位接着往下看。

Goldstein-Price修正法(G-P法)

针对\(H_k\)可能不正定的问题，如其不正定，就用“最速下降方向”来作为搜索方向。 \[d_{k+1}=\begin{cases}-H_kg_k，H_k正定\\-g_k，H_k非正定 \end{cases}\] 确定方向后，采用下列精确一维搜索（Armijo-Goldstein准则）（PS：其他准则也可以） \[1:f(\vec x_k+\alpha\vec d_k)≤f(\vec x_k)+\alpha c\nabla f(\vec x_k)\vec d_k\] \[2:f(\vec x_k+\alpha\vec d_k)≥f(\vec x_k)+\alpha (1-c)\nabla f(\vec x_k)\vec d_k\] 其中，\(c ∈ (0,0.5)\)

在一定条件下，G-P法全局收敛。但当海森矩阵\(H\)非正定情况较多时，收敛速度降为接近线性（梯度法）。

负曲率方向法

在鞍点处，我们有梯度为0，同时H矩阵不定，G-P等修正牛顿法无法在鞍点处继续，需要额外的方式：采用负曲率方向作为搜索方向，可使目标函数值下降。

负曲率方向，即\(d^T_k∇_2f(x_k)d_k<0\)。这里涉及到负曲率方向的定义，在张贤达的《矩阵分析与应用》中说到，当矩阵H是非线性函数\(f(x)\)的Hessian矩阵时，称满足\(p^HHp>0\) 的向量\(p\)为函数\(f\)的正曲率方向，满足\(p^HHp<0\)的向量\(p\)为函数\(f\)的负曲率方向。标量\(p^HHp\)称为函数\(f\)沿着方向\(p\)的曲率。沿着负曲率方向进行一维搜索必能使目标函数值下降。

因此，在鞍点时，我们只要找到任意负曲率方向，使函数值继续减小，就可以摆脱鞍点。

Goldfeld修正

与上面的Goldstein-Price修正的思路不同，Goldfeld在1966年也提出了一种方法，他的方法虽然还是在搜索方向\(d_k\) 上动手，但是当\(H_k\) 不正定时，他不是用最速下降方向\(-g_k\)来作为搜索方向，而是将\(d_k\)修正成下降方向——用下面的式子： \[{d_k} = - B_k^{ - 1}{g_k}\] 其中，\({B_k} = {H_k} + {E_k}\)是一个正定矩阵，\(E_k\)称为修正矩阵。在\(E_k\)满足一定条件的时候，（Goldfeld修正）牛顿法具有整体收敛性。

接下来我们使用\(B_k\)表示经过各种修正后的海森矩阵。

我们可以设\(E_k=v_kI\)，使得使得\(H_k+v_kI\)正定。比较理想的是，\(v_k\)为使\(H_k+v_kI\)正定的最小\(v\)。特征值的估计可以用盖尔圆盘定理或者Rayleigh商方法。

Gill-Murray改进Cholesky分解法的稳定牛顿法

前一个小节里说过加上一个\(v_kI\)使修正的Hessian矩阵正定，本节具体介绍一种高效稳定方法。Gill-Murray改进Cholesky分解法的稳定牛顿法可以说是修正牛顿法中的集大成。基本思想是当\(H_k\)不定矩阵时，采用修改Cholesky分解强迫矩阵正定；当\(g_k\)趋于0时，采用负曲率方向使函数值下降。如果矩阵正定且负曲率方向不存在，则是最优解。其大体流程如下

给定初始点\(x_0,\varepsilon>0,k=1\)

计算\(g_k,H_k\)

对\(H_k\)进行修改cholesky分解强迫矩阵正定，\(L_kD_kL_k^\ast=H_k+E_k\)

若\(||g_k||>\varepsilon\)（梯度不为0的计算机表达），则解方程\(L_kD_kL_k^\ast d_k=-g_k\Rightarrow d_k=-(H_k+E_k)^{-1}g_k\)，转步骤6

若\(||g_k||>\varepsilon\approx0\)，构造负曲率方向，如果构造不出负曲率方向，得到局部最优。

一维搜索步长，可用精确法或不精确准则得到步长\(\alpha_k\)，令\(x_{k+1}=x_k+\alpha_k d_k\)

若第5步满足终止条件，则OK；否则\(k=k+1\)，转2.

可以证明该方法总体收敛。现在还有其中最重要的一步：修正的Cholesky分解，请看笔记《海森矩阵校正》，更加详细的关于修正的Cholesky分解可见教程《Modified Cholesky Decomposition and Applications》

信赖域法

信赖域和line search同为最优化算法的基础算法，但是，从“Trust Region”这个名字你就可以看出，它是没有line search过程的，它是直接在一个region中“search”。在一维搜索中，从\(x_k\)点移动到下一个点的过程，可以描述为：\(x_k + α_k d_k\)，此处\(α_k d_k\)就是在\(d_k\)方向上的位移，可以记\(s_k\)，总的来说，线搜索把复杂问题改编成一系列简单的一维优化问题。而信赖域算法是根据一定的原则，直接确定位移\(s_k\)，同时，与一维搜索不同的是，它并没有先确定搜索方向\(d_k\)。如果根据“某种原则”确定的位移能使目标函数值充分下降，则扩大信赖域；若不能使目标函数值充分下降，则缩小信赖域。如此迭代下去，直到收敛。总的来说，信赖域方法是把复杂问题转化为一系列相对简单的局部优化问题。

信赖域算法的基本思想

在每次迭代中给出一个信赖域，这个信赖域一般是当前迭代点的一个小邻域。然后，在这个邻域内求解一个子问题，得到试探步长（trial step）\(s_k\)，接着用某一评价函数来决定是否接受该试探步以及确定下一次迭代的信赖域。如果试探步长被接受，则：\(x_{k+1}=x_k+s_k\)，否则，\(x_{k+1}=x_k\)。新的信赖域的大小取决于试探步的好坏，粗略地说，如果试探步较好，在下一步信赖域扩大或者保持不变，否则减小信赖域。

基本信赖域法

设当前点\(x_k\)的邻域定义为：\(\Omega_k\)。目标函数在极值点附近近似一个二次函数，因此对于无约束优化问题，利用二次逼近，构造如下信赖域子问题：设当前点\(x_k\)的邻域定义为：\(\Omega_k=\{x ∈ R^n|\Vert x-x_k\Vert≤\Delta_k\}\)，其中，\(\Delta_k\)称为信赖域半径。目标函数在极值点附近近似一个二次函数(有点类似插值法)，因此对于无约束优化问题，利用二次逼近，构造如下信赖域子问题： \[\min q^{(k)}(s)=f(x_k)+g_k^Ts+\frac{1}{2}(s^TB_ks)\\ s.t. \quad \|s\|_2≤\Delta_k\] 其中，\(s=x-s_k\)，\(g_k\)是目标函数\(f(x)\)在当前迭代点\(x_k\)处的梯度，\(B_k ∈ R^{n*n}\)对称，是\(f(x)\)在\(x_k\)处Hessian矩阵或Hessian矩阵的近似。

设\(s_k\)是信赖域子问题的解。目标函数\(f(x)\)在第k步的实际下降量（真实下降量）： \[Ared_k=f(x_k)-f(x_k+s_k)\] 二次模型函数\(q^{(k)}(s)\)的下降量（预测下降量）： \[Pred_k=q^{(k)}(0)-q^{(k)}(s_k)>0\] 定义比值： \[r=\frac{Ared_k}{Pred_k}\] 它衡量了二次模型与目标函数的逼近程度，\(r_k\)越接近于1，表明接近程度越好。因此，我们也用这个量来确定下次迭代的信赖域半径。

越接近与1，表明接近程度越好。这时可以增大\(\Delta_k\)以扩大信赖域；
\(r_k>0\)但是不接近于1，保持\(\Delta_k\)不变；
如果\(r_k\)接近于0，减小\(\Delta_k\)，缩小信赖域。 https://www.codelast.com/%E5%8E%9F%E5%88%9B%E4%BF%A1%E8%B5%96%E5%9F%9Ftrust-region%E7%AE%97%E6%B3%95%E6%98%AF%E6%80%8E%E4%B9%88%E4%B8%80%E5%9B%9E%E4%BA%8B/

Levenberg–Marquardt法

https://www.codelast.com/%e5%8e%9f%e5%88%9blm%e7%ae%97%e6%b3%95%e7%9a%84%e5%ae%9e%e7%8e%b0/

Dogleg法

https://zhuanlan.zhihu.com/p/99392484 https://zhuanlan.zhihu.com/p/101124802

拟牛顿法

拟牛顿法是求解非线性优化问题最有效的方法之一，于20世纪50年代由美国Argonne国家实验室的物理学家W.C.Davidon所提出来。Davidon设计的这种算法在当时看来是非线性优化领域最具创造性的发明之一。不久R. Fletcher和M. J. D. Powell证实了这种新的算法远比其他方法快速和可靠，使得非线性优化这门学科在一夜之间突飞猛进。

拟牛顿法的本质思想是改善牛顿法每次需要求解复杂的Hessian矩阵的逆矩阵的缺陷，它使用正定矩阵来近似Hessian矩阵的逆，从而简化了运算的复杂度。拟牛顿法和最速下降法一样只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足以产生超线性收敛性。这类方法大大优于最速下降法，尤其对于困难的问题。另外，因为拟牛顿法不需要二阶导数的信息，所以有时比牛顿法更为有效。如今，优化软件中包含了大量的拟牛顿算法用来解决无约束，约束，和大规模的优化问题。