牛顿法与拟牛顿法
本文首先从一个不带约束条件的最优化问题出发，导出用以解决此类问题的牛顿法，接着针对牛顿法涉及的计算复杂问题，阐述在此基础上改进得到的拟牛顿法的思路。
牛顿法
我们首先引入如下最优化问题：
$$min_{x\in{\mathbb{R}^n}}f(x)$$
若 $f(x)$ 在其定义域内具有二阶连续偏导数，在其定义域内任取一点 $x^{(k)}$，将 $f(x)$ 在点 $x^{(k)}$ 处进行二阶泰勒展开，得：
$$f(x)=f(x^{(k)})+g_k^T(x-x^{(k)})+\frac{1}{2}(x-x^{(k)})^TH_k(x-x^{(k)})$$
其中，${g_k}=\nabla_{x^{(k)}}f(x^{(k)})$，$H_k=(\frac{\partial^2 f(x^{(k)})}{\partial x_i \partial x_j})_{n \times n}$，该矩阵 $H_k$ 称为 Hessian 矩阵，当 $x \to x^{(k)}$ 时上式成立（即在 $x^{(k)}$ 的小邻域内成立）。由以上引入的最优化问题，我们需要得到 $x^{(k)}$ 小邻域内的极小值，而 $f(x)$ 在 $x^{(k)}$ 小邻域内取得极值的必要条件是 $\nabla_{x}f(x)=0$，在上式中令 $\nabla_{x}f(x)=0$，得：
$$\nabla_{x}f(x)={g_k}+H_k(x-x^{(k)})=0$$
当 Hessian 矩阵 $H_k$ 正定时（$H_k^{-1}$ 也正定），$f(x)$ 在 $x^{(k)}$ 小邻域内取得极小值，这也是牛顿法能够收敛的条件，由上式解得此时 $x=x^{(k)}-{H_k}^{-1}g_k$，将此时得到的 $x$ 置为新的 $x^{(k)}$，如此循环迭代，最后必能得到函数 $f(x)$ 的全局最小值。特别地，当 $x \in \mathbb{R}$ 时，以上过程如下图所示：
上图中的 $f(x_i),(i=1, 2, 3, 4, 5)$ 实际上指的是目标函数的一阶导数，而 $f’(x_i),(i=1,2,3,4,5)$ 则指的是目标函数的二阶导数。由上图可知，随着 $x$ 值的不断更新，使得一阶导数 $f(x_i)$ 逐渐接近其零点 $x_5$，而一阶导数为 $0$ 时，其对应的目标函数取得全局最优值。
我们将由牛顿法求目标函数全局最小值的过程概括成如下步骤：
输入：目标函数 $f(x)$，梯度 $g(x)$ 和 Hessian 矩阵 $H(x)$，精度 $\varepsilon$。
输出：目标函数 $f(x)$ 取得极小值时的最优值 $x^*$ 。
(1) 取初始点为 $x^{(0)}$，置 $k$ 为 $0$；
(2) 计算梯度 $g_k=g(x^{(k)})$，若 $\Arrowvert g_k \Arrowvert<\varepsilon$，则停止计算，令 $x^*=x^{(k)}$，否则继续下一步；
(3) 计算 Hessian 矩阵 $H_k=H(x^{(k)})$，按 $x^{(k+1)}=x^{(k)}-H_k^{-1}g_k$ 更新 $x$ 的值，令 $k=k+1$，返回步骤（2）。
泰勒公式反映了一种重要的微积分思想，那就是“以曲代曲”，即用简单曲线代替复杂曲线的求解过程。二阶泰勒展开是一种二次曲线，当目标函数也是二次函数时，或者目标函数不是二次函数，但是二次形态较强，亦或者迭代点已经进入极小点的邻域，那么牛顿法的收敛速度是非常快的。但是，对于非二次型的目标函数，有时会使函数值上升，即出现 $f(x^{(k+1)})>f(x^{(k)})$ 的情况，在严重的情况下甚至能够引起迭代点列 $\{x^{(k)}\}$ 发散而导致计算失败。针对这一问题，考虑在原始牛顿法的基础上引入迭代步长，这种方法称为 阻尼牛顿法，记最优的步长因子为 $\lambda_k(>0)$，其值通过以下过程获得：
$$\lambda_k=argmin_{\lambda \in \mathbb{R}}f(x^{(k)}-\lambda H_k^{-1}g_k)$$
将以上步骤 (3) 中的 $x^{(k+1)}=x^{(k)}-H_k^{-1}g_k$ 改为 $x^{(k+1)}=x^{(k)}-\lambda H_k^{-1}g_k$，原过程其他部分保持不变，即得阻尼牛顿法求目标函数全局最小值的过程。
Hessian 矩阵 $H_k$ 正定（$H_k^{-1}$ 也正定）是牛顿法收敛的重要条件，也就是说这样可以保证牛顿法的搜索方向是目标函数的下降方向，下面我们进一步说明这是为什么。若记 $p_k=-H_kg_k$，同时引入步长因子 $\lambda_k$，则 $x$ 的更新过程为：$x^{(k+1)}=x_k+\lambda_k p_k$，将此式代入目标函数的二阶泰勒展开式中，可得：
$$f(x^{(k+1)}) \thickapprox f(x^{(k)})-\lambda_k {g_k}^T {H_k}^{-1}g_k$$
由于 $H_k^{-1}$ 正定，由矩阵正定的定义可得 ${g_k}^T {H_k}^{-1}g_k > 0$，而步长因子 $\lambda_k > 0$，因此，$f(x^{(k+1)}) < f(x^{(k)})$，所以 Hessian 矩阵 $H_k$ 正定保证了牛顿法的搜索方向始终是目标函数的下降方向。
由以上过程可知，牛顿法在进行数值更新的时候，每次都要计算 Hessian 矩阵的逆 $H_k^{-1}$。这种计算是很复杂的。由此我们考虑是否可以用某个矩阵来代替这里的 Hessian 矩阵 $H_k$ 或者 $H_k^{-1}$，基于这种思想发展出了后来称之为拟牛顿法的优化方法。
拟牛顿法
我们分别考虑用矩阵 $G_k$ 来近似矩阵 $H_k^{-1}$，用矩阵 $B_k$ 来近似矩阵 $H_k$，分别得到 DFP 算法和 BFGS 算法。
DFP(Davidon-Fletcher-Powell) 算法
用矩阵 $G_k$ 来近似矩阵 $H_k^{-1}$，则由原始牛顿法的迭代过程可知，此时 $x^{(k+1)}=x^{(k)}-\lambda_kG_kg_k$，在不满足迭代终止条件的时候，为了继续得到 $x^{(k+2)}$ 使得迭代过程继续下去，那么接下来的一个问题是 $x=x^{(k+1)}$ 时 $G_{k+1}$ 是多少呢？不妨对目标函数 $f(x)$ 在 $x=x^{(k+1)}$ 处做二阶泰勒展开，得：
$$f(x)=f(x^{(k+1)})+g_{k+1}^T(x-x^{(k+1)})+\frac{1}{2}(x-x^{(k+1)})^TH_{k+1}(x-x^{(k+1)})$$
上式两边分别对 $x$ 求梯度，得：
$$\nabla_xf(x)=g_{k+1}+H_{k+1}(x-x^{(k+1)})$$
上式中，因为 $x^{(k)}$ 在 $x^{(k+1)}$ 的一个小邻域内，故可取 $x=x^{(k)}$，代入上式，可得：
$$g_{k+1}-g_k=H_{k+1}(x^{(k+1)}-x^{(k)})$$
记 $y_k=g_{k+1}-g_k$，$\delta_k=x^{(k+1)}-x^{(k)}$，上式变为：
$$y_k=H_{k+1}\delta_k$$
或：
$$\delta_k=H_{k+1}^{-1}y_k$$
以上两式称为 拟牛顿条件。
当使用 $G_{k+1}$ 来近似 $H_{k+1}^{-1}$ 时，它必须满足上文所述的拟牛顿条件，即 $\delta_k=G_{k+1} \cdot y_k$。对 $G_k$ 进行更新得到 $G_{k+1}$，令 $G_{k+1} = G_k + \Delta G$，考虑 $\Delta G$ 为两个附加项 $P_k$、$Q_k$ 之和，即有：
$$G_{k+1}=G_k+P_k+Q_k$$
上式等号两边同时乘以 $y_k$ 可得：
$$G_{k+1}y_k=G_k y_k + P_k y_k + Q_k y_k$$
上式等号左端为 $\delta_k$，若令右端 $P_k y_k = \delta_k$，则必须有 $Q_k y_k = -G_k y_k$，为满足上述条件，可取 $P_k=\frac{\delta_k{\delta_k}^T}{\delta_k^T y_k}$，$Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k}$，则有：
$$G_{k+1}=G_k+\frac{\delta_k{\delta_k}^T}{\delta_k^T y_k}-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k}$$
以上就是 DFP 算法用到的 $G_{k+1}$ 的迭代公式。细心的读者可能发现此时 $G_{k+1}$ 是否正定呢？事实上，只要初始的 $G_0$ 是正定的，后续迭代得到的每个 $G_{k+1}$ 都是正定的。
经过以上对 $G_k$ 迭代过程的分析，我们得到 DFP 算法的具体步骤如下：
输入：目标函数 $f(x)$，梯度 $g(x)$，精度 $\varepsilon$。
输出：目标函数 $f(x)$ 取得极小值时的最优值 $x^*$ 。
(1) 取初始点为 $x^{(0)}$，取 $G_0$ 为正定对称矩阵，置 $k$ 为 $0$；
(2) 计算梯度 $g_k=g(x^{(k)})$，若 $\Arrowvert g_k \Arrowvert<\varepsilon$，则停止计算，得 $x^*=x^{(k)}$，否则继续下一步；
(3) 令 $p_k=-G_kg_k$，执行一维搜索，求步长因子 $\lambda_k$ 使得 $\lambda_k=argmin_{\lambda}f(x^{(k)}+\lambda p_k)$；
(4) 计算 $x^{(k+1)}=x^{(k)}+\lambda_k p_k$ ；
(5) 计算梯度 $g_{k+1}=g(x^{(k+1)})$，若 $\Arrowvert g_{k+1} \Arrowvert<\varepsilon$，则停止计算，得 $x^*=x^{(k+1)}$，否则按式 $G_{k+1}=G_k+\frac{\delta_k{\delta_k}^T}{\delta_k^T y_k}-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k}$ 计算 $G_{k+1}$，令 $k=k+1$，返回步骤（3）。
BFGS(Broyden-Fletcher-Goldfarb-Shanno) 算法
用矩阵 $B_k$ 来近似矩阵 $H_k$ 的时候，同上文所述，需要 $B_{k+1}$ 满足以上提到的拟牛顿条件，即：$y_k=B_{k+1}\delta_k$。同理考虑：
$$B_{k+1}=B_k+P_k+Q_k$$
上式等号两边同乘 $y_k$，得：
$$B_{k+1}\delta_k=B_k\delta_k+P_k\delta_k+Q_k\delta_k$$
如上文所述一样，可令 $P_k\delta_k=y_k$，$Q_k\delta_k=-B_k\delta_k$，取 $P_k=\frac{y_ky_k^T}{y_k^T\delta_k}$，$Q_k=-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$，得：
$$B_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$$
以上即为 $B_{k+1}$ 更新的迭代公式。只要初始矩阵 $B_0$ 是正定的，这个迭代过程得到的每个 $B_{k+1}$ 都是正定的。
经过以上过程，得到 BFGS 算法如下：
输入：目标函数 $f(x)$，梯度 $g(x)$，精度 $\varepsilon$。
输出：目标函数 $f(x)$ 取得极小值时的最优值 $x^*$ 。
(1) 取初始点为 $x^{(0)}$，取 $B_0$ 为正定对称矩阵，置 $k$ 为 $0$；
(2) 计算梯度 $g_k=g(x^{(k)})$，若 $\Arrowvert g_k \Arrowvert<\varepsilon$，则停止计算，得 $x^*=x^{(k)}$，否则继续下一步；
(3) 令 $B_kp_k=-g_k$，求出 $p_k$，执行一维搜索，求步长因子 $\lambda_k$ 使得 $\lambda_k=argmin_{\lambda}f(x^{(k)}+\lambda p_k)$；
(4) 计算 $x^{(k+1)}=x^{(k)}+\lambda_k p_k$ ；
(5) 计算梯度 $g_{k+1}=g(x^{(k+1)})$，若 $\Arrowvert g_{k+1} \Arrowvert<\varepsilon$，则停止计算，得 $x^*=x^{(k+1)}$，否则按式 $B_{k+1}=B_k+\frac{y_k{y_k}^T}{y_k^T \delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$ 计算 $B_{k+1}$，令 $k=k+1$，返回步骤（3）。
Broyden 类算法
在介绍 Broyden 类算法之前，首先介绍一下 Sherman-Morrison 公式。
Sherman-Morrison 公式
假设 $A$ 是 n 阶可逆矩阵，$u、v$ 是 n 维向量，且 $A+uv^T$ 也是可逆矩阵，则：
$$(A+uv^T)^{-1}=A^{-1}-\frac{A^Tuv^TA^{-1}}{1+v^TA^{-1}u}$$
上式称为 Sherman-Morrison 公式。
对于上文所述的 DFP 算法 和 BFGS 算法 而言，容易验证：$G_k=B_k^{-1}$，$G_{k+1}=B_{k+1}^{-1}$，对 BFGS 算法 的迭代公式 $B_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$ 两次应用 Sherman-Morrison 公式可得：
$$G_{k+1}=(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})G_k(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})^T+\frac{\delta_k\delta_k^T}{\delta_k^Ty_k}$$
以上称为 BFGS 算法 关于 $G_{k+1}$ 的迭代公式。
若将 DFP 算法 中关于 $G_{k+1}$ 的迭代公式得到的 $G_{k+1}$ 记为 $G^{DFP}$，以上 BFGS 算法 得到的 $G_{k+1}$ 记为 $G^{BFGS}$，它们都满足拟牛顿条件，故而它们的线性组合：
$$G_{k+1}=\alpha G^{DFP}+(1-\alpha)G^{BFGS} (0\leq\alpha\leq1)$$
也满足拟牛顿条件，并且是正定的。这样便可以得到一类拟牛顿法，称之为 Broyden 类算法。