优化算法--牛顿迭代法

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

方程求根

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

牛顿迭代法

先考虑一个小问题：求解方程$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. 猜一个初始值，因为根大概是1点多吧，那就给个$x_0=2$好了；
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$$

算法优缺点分析

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

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

优化问题求解

优化问题从泛函的角度理解起来，就是计算函数之间的距离最小。对于距离的定义有很多，比较常用的是二范数，使二范数距离最小的求解过程就叫做最小二乘。对于$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

几个改进方法

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

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

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

1. 拟牛顿法：为了避免计算Hessian矩阵，不直接计算$H$，而是构造一个矩阵$K$来近似，$K$需要一直正定并且更新起来比较简单，此处可以查看相关文献，不赘述了；

2. 高斯牛顿法：将目标函数$\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$的思路就不行了，此时，高斯-牛顿法并不收敛。

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

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

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

Reference

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