数据、特征和数值优化算法是机器学习的核心，而牛顿法及其改良（拟牛顿法）是机器最常用的一类数字优化算法，今天就从牛顿法开始，介绍几个拟牛顿法算法。本博文只介绍算法的思想，具体的数学推导过程不做介绍。
1. 牛顿法
牛顿法的核心思想是”利用函数在当前点的一阶导数，以及二阶导数，寻找搜寻方向“（回想一下更简单的梯度下降法，她只用了当前点一阶导数信息决定搜索方向）。
牛顿法的迭代公式是（稍微有修改，最原始的牛顿法\(\gamma=1\)：
\[{{\bf{x}}_{n + 1}} = {{\bf{x}}_n} - \gamma {[Hf({{\bf{x}}_n})]^{ - 1}}\nabla f({{\bf{x}}_n})\]
其中，\({[Hf({{\bf{x}}_n})]^{ - 1}}\nabla f({{\bf{x}}_n})\)是线搜索方向。这个方向的含义是什么呢。有两种物理解释：
1. 一阶导数\(\nabla f({\bf{x}}){\mkern 1mu} \)当前搜寻点 与 \(\nabla f({\bf{x}}){\mkern 1mu} = 0\)连线的方向。
2. 当前点泰勒展开（舍弃二阶以上项）函数中 当前搜寻点 与 泰勒展开函数极小值连线方向。
2. 拟牛顿算法
上述的牛顿法需要计算Hessian矩阵的逆矩阵，运算复杂度太高。在动辄百亿、千亿量级特征的大数据时代，模型训练耗时太久。因此，很多牛顿算法的变形出现了，这类变形统称拟牛顿算法。拟牛顿算法的核心思想用一个近似矩阵\(B\)替代逆Hessian矩阵\({H^{ - 1}}\)。不同算法的矩阵\(B\)的计算有差异，但大多算法都是采用迭代更新的思想在tranning的没一轮更新矩阵\(B\)。
下面介绍两种最常用的拟牛顿算法：
- BFGS
BFGS的算法流程如下：
1. 初始化：初始点\({x_0}\)以及近似逆Hessian矩阵\(B_0^{ - 1}\)。通常，\({B_0} = I\),既为单位矩阵。
2. 计算线搜索方向：\({{\bf{p}}_k} = - B_k^{ - 1}\nabla f({{\bf{x}}_k})\)
3. 用”Backtracking line search“算法沿搜索方向找到下一个迭代点：\({{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\alpha _k}{{\bf{p}}_k}\)
4. 根据Armijo–Goldstein 准则，判断是否停止。
5. 计算\({{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\alpha _k}{{\bf{p}}_k}\); 以及 \({{\bf{y}}_k} = \nabla f({{\bf{x}}_{k + 1}}) - \nabla f({{\bf{x}}_k})\)
6. 迭代近似逆Hessian矩阵：
\[B_{k + 1}^{ - 1} = \left( {I - \frac{{{s_k}y_k^T}}{{y_k^T{s_k}}}} \right)B_k^{ - 1}\left( {I - \frac{{{y_k}s_k^T}}{{y_k^T{s_k}}}} \right) + \frac{{{s_k}s_k^T}}{{y_k^T{\mkern 1mu} {s_k}}}\]
上式5中的推到方法比较复杂，有兴趣的可以搜一下相关文献。
- L-BFGS
BFGS算法需要存储近似逆Hessian矩阵\(B_0^{ - 1}\)。对于很多应用来说（比如百度的CTR预估），千亿量级的feature数需要\({10^{16}}\)p存储。显然，目前的计算机技术还很难满足这样的空间需求。因此，内存受限的BFGS算法(Limited-memory BFGS)就诞生了。
L-BFGS算法不存储近似逆Hessian矩阵\(B_0^{ - 1}\)， 而是直接通过迭代算法获取本轮的搜索方向\({{\bf{p}}_k} = - B_k^{ - 1}\nabla f({{\bf{x}}_k})\)。
虽然L-BFGS算法不需要保存\(B_0^{ - 1}\)矩阵，但要保存每次迭代的中间信息，不过都是一维数组，而迭代次数有不会很多（100次以内），所以存储要求大大降低。
L-BFGS算法的流程如下：
定义:
\[{s_k} = {x_{k + 1}} - {x_k}{\mkern 1mu} \]
\[{y_k} = {g_{k + 1}} - {g_k}{\mkern 1mu} \]
以上两个一维数组需要保存每次迭代的结果。
再定义：
\[{\rho _k} = \frac{1}{{y_k^{\rm{T}}{s_k}}}\]
\[{g_k} \equiv \nabla f({x_k})\]
L-BFGS算法如下：
1. 初始化：初始点\({x_0}\)以及近似逆Hessian矩阵\(B_0^{ - 1}\)。通常，\({B_0} = I\),既为单位矩阵。
2. 用”Backtracking line search“算法沿搜索方向找到下一个迭代点：\({{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\alpha _k}{{\bf{p}}_k}\)。
3. 根据Armijo–Goldstein 准则，判断是否停止。
4. 更新搜索方向 \({{\bf{p}}_k}\); 用下面的”two loop recursion"算法。
two loop recursion算法：
\[\begin{array}{l}
q = {g_k}{\mkern 1mu} \\
For(i = k - 1,k - 2, \ldots ,k - m):\\
{\alpha _i} = {\rho _i}s_i^{\rm{T}}q{\mkern 1mu} \\
q = q - {\alpha _i}{y_i}{\mkern 1mu} \\
{H_k} = y_{k - 1}^{\rm{T}}{s_{k - 1}}/y_{k - 1}^{\rm{T}}{y_{k - 1}}\\
z = {H_k}q\\
done\\
For{\rm{(}}i = k - m,k - m + 1, \ldots ,k - 1):\\
{\beta _i} = {\rho _i}y_i^{\rm{T}}z{\mkern 1mu} \\
z = z + {s_i}({\alpha _i} - {\beta _i}){\mkern 1mu} \\
done\\
{{\bf{p}}_k} = B_k^{ - 1}{g_k} = z
\end{array}\]
参考文献：
1. http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
2. http://en.wikipedia.org/wiki/Limited-memory_BFGS
3. http://aria42.com/blog/2014/12/understanding-lbfgs/