优化算法--牛顿迭代法

作者: MunCN | 来源:发表于2019-03-25 22:50 被阅读0次

优化算法--牛顿迭代法
L-BFGS算法
无约束凸优化算法
记--平方根的算法
每日一问之初识牛顿迭代法（Newton's method）
牛顿迭代法-最优化方法
1.3求根之牛顿迭代法
Logistic回归与最大熵模型-优化算法
牛顿法开根
吹水牛顿迭代法

博客搬家至 Mun： https://kiddie92.github.io

简书同步更新

牛顿法给出了任意方程求根的数值解法，而最优化问题一般会转换为求函数之间在"赋范线性空间"的距离最小点，所以，利用牛顿法去求解任意目标函数的极值点是个不错的思路。

方程求根

对于一元二次方程，求根其实很简单，只要套用求根公式就行了，但找到一个方程的求根公式（解析解）其实是很困难的，可以证明5次方程以上便没有解析解了，参考维基百科五次方程。其他的复杂方程如偏微分方程求解更是超级困难。好在随着计算机技术的发展，解析解变的不再那么重要（至少是在工程上），取而代之的方法便是数值解法，牛顿法便是众多数值解法中的一个。
数值法求解又叫做数值分析，主要利用逼近的思想来使数值解通过迭代计算不断接近解析解，而得出来得解就叫做数值解，在工程上，数值解只要是在精度要求范围内满足方程便是有用的。

牛顿迭代法

sqrt2.png

先考虑一个小问题：求解方程 $x^2-2=0$ 的根，也即求解 $\sqrt 2$ 。牛顿迭代法的思想从几何的角度很好理解，如上图所示（画图的脚本在这里）方程的根就是函数 $y=x^2-2$ 与 $x$ 轴的交点处横坐标的值。从图中 $x_n$ 点出发，计算函数在 $x_n$ 点处的切线，再计算切线和 $x$ 轴的交点得到 $x_{n+1}$ ，再计算函数在 $x_{n+1}$ 点处的切线... 一直这样迭代下去，可以发现 $x_{n}$ 会越来越接近方程的根。

上述思路的数学表达：
由 $x_{n}$ 计算 $y_{n}$

$f(x_n)=y_n$

得到切线方程：

$y-y_n= \left. f(x)' \right | _{x=x_n}(x-x_n)$

切线和 $x$ 轴的交点，也即，当 $y=0$ 时，

$0-y_n=\left. f(x)' \right | _{x=x_n}(x-x_n)$

$\frac{-y_n}{\left. f'(x) \right | _{x=x_n}} = (x-x_n)$

当 $\left. f'(x) \right | _{x=x_n} \neq 0$ 时，

$x = x_n- \frac{y_n}{\left. f'(x) \right | _{x=x_n}}$

由 $y_n = f(x_n)$ ，得到：

$x = x_n- \frac{f(x_n)}{\left. f'(x) \right | _{x=x_n}}$

令 $x=x_n$ ，继续迭代，则得到迭代公式：

$x_{n+1} = x_n- \frac{f(x_n)}{\left. f'(x) \right | _{x=x_n}}$

推导过程还可以从函数泰勒展开的角度去理解，这在很多博客里有写，这里就不赘述了。

根据上面的迭代公式，可以计算方程 $x^2-2=0$ 的根了：

猜一个初始值，因为根大概是1点多吧，那就给个 $x_0=2$ 好了；
计算 $x_1$ ：
$x_{1} = x_0- \frac{f(x_0)}{\left. f'(x) \right | _{x=x_0}}= 2- \frac{f(2)}{\left. f'(x) \right | _{x=2}}=1.5$
$x_{2} = 1.5- \frac{f(1.5)}{\left. f'(x) \right | _{x=1.5}}=1.416667$
$x_{3} = 1.416667- \frac{f(1.416667)}{\left. f'(x) \right | _{x=1.416667}}=1.414216$

算法优缺点分析

牛顿法的优点当然就是提供了一种方程求根的数值解方法。而缺点也有几点：

首先算法是要求函数处处可导的，如果对于优化问题还需要导函数连续（因为要求处处存在二阶导数），否则算法就不能计算函数的根了，比如 $f(x)=x^{1/3}$ 就不能收敛，虽然函数的根为0，但是它在0处的导数是不存在的；
求出的解可能仅仅是众多解中的一个，这个比较依赖于初始值的选取，比如上面的问题，初始值为2，则收敛到了方程的正数解，要想得到负数解，则需要将初始值选在负数中，现实中的问题，很难去估计解的大小范围；
如果初始的估计值与根的距离太远收敛就会变的比较慢；
要求每次迭代是得到的切线导数不能为0，如推导过程所示；
如果方程本来就没有根，那牛顿法是不能收敛的；

优化问题求解

优化问题从泛函的角度理解起来，就是计算函数之间的距离最小。对于距离的定义有很多，比较常用的是二范数，使二范数距离最小的求解过程就叫做最小二乘。对于 $Gm=data_{predict}$ 这样的线性问题（非线程问题可以通过泰勒展开转换成线性问题），可以定义距离为 $\phi (m)=||Gm-data_{observation}||_2$ ，为了求距离最小值点，需要先求极值点，问题便转换为求解 $\phi '(m)=0$ 的根，这时候牛顿法便派上了用场。与之前问题不同的是，这里需要求 $\phi '(m)$ 的导数，也即求解 $\phi "(m)$ ，也即Hessian矩阵。假设，此处的参数 $m$ 是n维向量，则Hessian矩阵为：

$H = \begin{pmatrix} \frac{\partial ^2f}{\partial m_1^2} & \frac{\partial ^2f}{\partial m_1 \partial m_2} & \cdots & \frac{\partial ^2f}{\partial m_1 \partial m_n} \\ \frac{\partial ^2f}{\partial m_2 \partial m_1} & \frac{\partial ^2f}{\partial m_2^2} & \cdots & \frac{\partial ^2f}{\partial m_2 \partial m_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial ^2f}{\partial m_n \partial m_1} & \frac{\partial ^2f}{\partial m_n \partial m_2} & \cdots & \frac{\partial ^2f}{\partial m_n^2} \\ \end{pmatrix}$

所以，牛顿法求解最优化问题，需要先求目标函数的Jacobian矩阵和Hessian矩阵，计算量比较大的便是计算Hessian矩阵了，因为二阶导计算量成指数增长。

注意，这里若二阶导数是连续的，则 $H$ 是对称矩阵。

算法步骤

步骤1： 给定误差阈值 $0\leq\epsilon <<1$ ，初始模型 $m_0$ （也可以给定迭代次数）；
步骤2： 计算梯度 $g_k=\nabla f(m_k)$ ，若 $g_k\leq \epsilon$ ，停止计算，输出 $m^* \approx m_k$ ;
步骤3： 计算Hessian矩阵 $G_k=\nabla^2f(m_k)$ ，计算 $d_k=\frac {g_k}{G_k}$ ;
步骤4： 令 $x_{k+1}=x_{k}-d_{k}$ ，k=k+1，转到第2步。

示例

一个例子：求极小值 $f(m_1,m_2) = -m_1^3-m_2^3+3m_1^2+2m_2^2+m_1+m_2-1$

主要代码如下所示，完整代码请查看这里

while ((k < n) or np.sqrt(np.power(gk[0,0],2)+np.power(gk[1,0],2))) > e:
    #向前差分计算一阶导
    gk[0,0] = 1/m1stp*(func(m1+m1stp,m2)-func(m1,m2))
    gk[1,0] = 1/m2stp*(func(m1,m2+m2stp)-func(m1,m2))
    
    #向前差分计算海森矩阵,注意：以函数为二阶导连续为前提
    Gk[0,0] = 1/m1stp*(func(m1+m1stp,m2)-2*func(m1,m2)+func(m1-m1stp,m2))
    Gk[0,1] = 1/(m1stp*m2stp)*(func(m1+m1stp,m2+m2stp)-func(m1,m2+m2stp)-func(m1+m1stp,m2)+func(m1,m2))
    Gk[1,0] = Gk[0,1]
    Gk[1,1] = 1/m2stp*(func(m1,m2+m2stp)-2*func(m1,m2)+func(m1,m2-m2stp))

    dk = Gk.I*gk

    #修正模型
    m1 = m1-dk[0,0]
    m2 = m2-dk[1,0]

    k = k+1

几个改进方法

优化算化考虑重点包括算法的通用性、有效性、收敛性、效率，当然，这些都包括在时间复杂度和空间复杂度中。牛顿法存在几个问题需要考虑一下：

计算Hessian矩阵太耗资源和时间了；
牛顿法不稳定，只有 $H$ 正定时才收敛，也即要求目标函数的 Hessian 阵 $H_{i,j}$ 在每个迭代点 $m_i$ 处是正定的，否则难以保证牛顿法收敛的方向，实际上， $H$ 很可能是一个病态/奇异矩阵；
初始模型 $m_0$ 很重要，选的不好会迭代很多次，收敛比较慢；
初始模型的选取不在最小值附近，很容易让结果陷入局部极小值。

对此，大牛们提出了一些改进的方法：

拟牛顿法：为了避免计算Hessian矩阵，不直接计算 $H$ ，而是构造一个矩阵 $K$ 来近似， $K$ 需要一直正定并且更新起来比较简单，此处可以查看相关文献，不赘述了；
高斯牛顿法：将目标函数 $\phi (m)=||Gm-data_{observation}||_2$ 变换为 $\phi (m)=\frac {1}{2}||r(m)||_2$ ，其中 $r(m)$ 表示残差（residual），则根据chain rule，可以得到：
$\nabla^2 \phi (m)=\nabla r(m) \nabla^T r(m)+\sum_{i=1}^nr_i(m)\nabla^2 r_{i}(m)$
这里令 $Q(m)=\sum_{i=1}^nr_i(m)\nabla^2 r_{i}(m)$ ，若对于将要迭代的值 $m^*$ ，有 $r_{i}(m^*)=0$ 则 $Q(m)=0$ ；这样的话就不需要计算Hessian矩阵了。这个想法不错，当 $m^*$ 和极值点/最小值的距离比较近时，简直完美；但是，当初始值距离最小值较远时， $Q(m) \approx 0$ 的思路就不行了，此时，高斯-牛顿法并不收敛。

所以高斯-牛顿法也是极度依赖初始模型/初值的选取的

莱文贝格－马夸特方法(Levenberg–Marquardt algorithm)：该方法结合了高斯-牛顿法和最速下降法/梯度法，因为高斯-牛顿法比较依赖初始模型/初值，梯度法可以克服这个问题；而梯度法收敛速度要低于高斯-牛顿法，所以该方法能提供数非线性最小化（局部最小）的数值解。其实做法也很简单，就是在目标函数内加了一个参数 $\lambda$ ，所以该方法也叫做阻尼最小二乘法。类似的做法在Tikhonov正则化中也出现了。

所有这些方法都可能陷入局部极小值，而非找到全局极小值/最小值。要想克服这个问题，就需要启发式/非线性优化算法了。

Reference

Calculus- where Newton's method fail
Newton's Method
马昌凤, 《最优化方法及其 Matlab 程序设计》

优化算法--牛顿迭代法
博客搬家至 Mun： https://kiddie92.github.io 简书同步更新牛顿法给出了任意方程求根...
L-BFGS算法
BFGS算法是用来求解最优化问题的，在这个算法中，相对于普通的牛顿迭代法有很大的改进。链接：http://blog...
无约束凸优化算法
本章涉及知识点1、scipy库求解全局最优和局最优2、多元函数的极值求解算法3、牛顿迭代法算法4、牛顿迭代法求解多...
记--平方根的算法
Java实现牛顿迭代法: C/C++实现3d游戏引擎算法实现1/sqrt(x),改一下返回值成为sqrt()算法:...
每日一问之初识牛顿迭代法（Newton's method）
什么是牛顿迭代法？今天在刷 LeetCode 的 sqrt(x) 这道题的时候，看到别人的解法中有使用牛顿迭代法...
牛顿迭代法-最优化方法
原理五次及以上多项式方程没有根式解，这个是被伽罗瓦用群论做出的最著名的结论。但牛顿迭代法将探索高阶函数的问题简化...
1.3求根之牛顿迭代法
目录 [TOC] 前言今天我们讲的是具有收敛速度快，能求重根的解方程之法，牛顿迭代法。（一）牛顿迭代法的分析 ...
Logistic回归与最大熵模型-优化算法
Logistic回归与最大熵模型-理论推导中提到了4个优化算法：分别是：梯度下降算法拟牛顿法（牛顿法）通用迭...
牛顿法开根
牛顿迭代法（Newton's method）又称为牛顿-拉夫逊方法（Newton-Raphson method）。...
吹水牛顿迭代法
因为吹水的能力不佳，所以要先打个草稿，今天的吹水过程大概是：1、牛顿迭代法的演绎过程2、牛顿迭代法求n次方根3、牛...