求局域最小值的方法 学习笔记

在之前的leetcode周赛里有这么一题，要求求平面上某点，到各个给定点的欧几里得距离之和最小，即求给定函数$f(x,y)={\sum }_{i=1}^n\sqrt{(x-x_i{)}^2+(y-y_i{)}^2}$在平面上的最小值。频道里的大佬说了能直接调包来做，查了一下，python有scipy的optimize库可以用不同的算法来求解函数局域最小值。

好的，非常感谢你的什么也没讲。那么就从这些缩写开始出发，开始摸一下求局域最小值。

最速下降法

最速下降法顾名思义，就是在$(x_k,y_k)$处，沿着函数$f(x,y)$下降最快的方向，即负梯度方向走上一段。迭代公式为$${\vec{x}}_{k+1}={\vec{x}}_k-{\rho }_k{\vec{g}}_k$$.
步长${\rho }_k$需满足$f({\vec{x}}_{k+1})<f(\vec{x})$
我们不妨先试着使用Leetcode上AC了的算法，使用固定步长的方法来求解上述Rosenbrock函数的最小值，即初始步长固定为${\rho }_0=0.03$,当$f({\vec{x}}_{k+1})>f({\vec{x}}_k)$时，将$\rho$减小一半直至步长满足函数下降要求。这里有一点需要注意，因为Rosenbrock函数的在山谷之外其梯度的模特别大，所以我们应该使用的是单位长度的梯度${\vec{g}}_k$来计算下一步的位置。(参考：$x=1.05,y=1.05$时，对应梯度的模为$|{\vec{g}}_k|=24.5127$)
结果：

不出意料的是在迭代10万次之后，函数甚至离目标值还有很大一段差距，在查看步长迭代的情况后，猜测是因为步长变化太剧烈，使得函数下降太慢：

那么接下来对步长提出要求，不仅需要使得函数严格递减，还要使函数递减到该方向上的最小值。也就是找到$\rho$使得$min(f({\vec{x}}_k-\rho {\vec{g}}_k))$。
令$g(\rho )=f({\vec{x}}_k-\rho {\vec{g}}_k)$，则$g(\rho )$的最小值在$$\frac{d}{d\rho }g(x)=\nabla f({\vec{x}}_k-\rho {\vec{g}}_k{)}^T{\vec{g}}_k=0$$处取到。这只需我们解一个一元方程即可得到最适宜步长$\rho$。
在求解$\rho$时，有个较为实用的方法，即将$f(x)$展开为二次型，$$f^{\ast }({\vec{x}}_k+\vec{\delta })=\frac{1}{2}{\vec{\delta }}^{T}H\vec{\delta }+{\vec{g}}^T\vec{\delta }+f({\vec{x}}_k)$$
若海森矩阵$H$为正定矩阵，则
$$\begin{array}{rl}\nabla & f^{\ast }({\vec{x}}_k-\rho {\vec{g}}_k{)}^T{\vec{g}}_k\\ =& {\vec{g}}_{k+k}^T\cdot {\vec{g}}_k\\ =& (H(\vec{x}-\rho {\vec{g}}_k)+{\vec{g}}_k{)}^T{\vec{g}}_k\\ =& (H\vec{x}+\vec{b}-\rho H{\vec{g}}_k{)}^T{\vec{g}}_k\\ =& ({\vec{g}}_k-\rho H{\vec{g}}_k{)}^T{\vec{g}}_k=0\end{array}$$
于是有$\rho =\frac{{\vec{g}}_k^T{\vec{g}}_k}{{\vec{g}}_k^{T}H{\vec{g}}_k}$，此时不再需要计算单位长度的梯度，而是使用${\vec{g}}_k$直接参与计算，在基础上重新计算Rosenbrock函数的最小值，则有：

可以很明显的看出在只使用原来的1/4迭代次数，就达到了最小值。在某些点处可以甚至有一个很明显的跳跃。对比了步长$\rho (k)$，发现固定步长法输的透彻。

牛顿法

在上述讨论最速下降法的时候，我们很容易发现面对这种存在盆地的函数，最速下降法收敛速度及其缓慢，那么从上面的函数展开式出发我们引入牛顿法。
将$f(\vec{x})$在${\vec{x}}_k$处展开为二次型，$$f^{\ast }({\vec{x}}_k+\vec{\delta })=\frac{1}{2}{\vec{\delta }}^{T}H\vec{\delta }+{\vec{g}}^T\vec{\delta }+f({\vec{x}}_k)$$
在函数局域最小值点有
$$\begin{array}{rl}\nabla & f^{\ast }({\vec{x}}_{k+1})=0\\ & =H_k\vec{\delta }+{\vec{g}}_k\\ & =0\end{array}$$
得迭代公式${\vec{x}}_{k+1}={\vec{x}}_k-H_k^{-1}{\vec{g}}_k$，尝试一下：

下降速度之快出乎我的意料，观察一下函数曲线$f({\vec{x}}_k)$：

第三次迭代的时候出现了一个上升，那么得到的第一个结论就是：牛顿法并不是一个下降算法，即他并不保证每一步的搜索都能保证函数值下降。
具体到这个问题中，原因是虽然搜索的方向${\vec{d}}_k$确实是函数下降方向，但是由于牛顿法的步长默认为1，出现了虽然是展开函数$f^{\ast }(\vec{x})$的最小值点${\vec{x}}^{\ast }$，但是不能使原函数$f(\vec{x})$下降。
另一个可能的会使函数值不降反升的可能是，当前点远离最小值点${\vec{x}}^{\ast }$使得海森矩阵$H$并非正定矩阵，则搜索方向$d_k=-H_k^{-1}{g}_k$不一定是函数的下降方向，即$g_k^T{d}_k<0$
还有一个关于牛顿法的较为严重的问题是，函数可能不能二阶可微，此时无法通过牛顿法求解。
为了解决上述几个问题，分别引入了带阻尼的牛顿法（解决步长问题），拟牛顿法（解决求海森矩阵以及海森矩阵不正定的问题）