【最优化】基于导数的最优化方法

前言

在 【最优化】非线性优化及KKT条件 一文中，我们介绍了非线性优化问题的必要条件和充分条件，通过这些最优性条件，我们可以在数学上和理论上求解到显性的或者解析的最优解。但现实中待优化的目标函数可能较为复杂，并且考虑到希望利用计算机来求解的工程场景，基于迭代的方法来求出最优解不失为最佳的选择。

本文将着眼于无约束的非线性最优化问题：

$$\min\;f(\boldsymbol x),\quad\boldsymbol x\in\Bbb R^n$$

其中，$f(\boldsymbol x)$ 是定义在$\Bbb R^n$ 上的具有一阶连续偏导数的实函数。

介绍基于目标函数的导数，通过一系列一维搜索/线搜索（见文末附赠章节），不断迭代最终达到最优值的方法。这类方法的核心思想总是：希望从某一点出发，选择一个使得当前目标函数值下降最快的方向对起始点进行更新，以尽快到达极小点。

最速下降法

最速下降法（steepest descent）最早由法国数学家 Augustin-Louis Cauchy¹ 在 1847 年提出，Jacques Hadamard 也于 1907 年独立地提出了一种类似的方法。在 1944 年Haskell Curry 首次研究了该方法在非线性优化问题中的收敛性，此后该方法得到了越来越多的研究和应用。

最速下降方向

最速下降法的核心思想是找到在当前点$\boldsymbol x$ 处使得$f(\boldsymbol x)$ 下降最快的方向$\boldsymbol d$ . 因为函数沿方向$\boldsymbol d$ 的变化率可以用方向导数$\text{D}f(\boldsymbol x;\boldsymbol d)$ 表示，并且：

$$\text{D}f(\boldsymbol x;\boldsymbol d)=\nabla f(\boldsymbol x)^\top\boldsymbol d$$

所以，要找到下降最快的方向等价于最小化方向导数。由于只是找方向（对其模长不做要求），因此不是一般性地，我们可以约束$\Vert\boldsymbol d\Vert\leq1$。即找最速方向归结于求解另一个非线性优化问题：

$$\begin{aligned} \min \quad&\nabla f(\boldsymbol x)^\top\boldsymbol d\\ \text{s.t.}\quad&\Vert\boldsymbol d\Vert\leq1 \end{aligned}$$

【求解】利用 Schwartz 不等式：

$$|\nabla f(\boldsymbol x)^\top\boldsymbol d|\leq\Vert\nabla f(\boldsymbol x)\Vert\ \Vert\boldsymbol d\Vert\leq\Vert\nabla f(\boldsymbol x)\Vert$$

去掉绝对值，有

$$\nabla f(\boldsymbol x)^\top\boldsymbol d\geq-\Vert\nabla f(\boldsymbol x)\Vert$$

从而当

$$\boldsymbol d=-\frac{\nabla f(\boldsymbol x)}{\Vert\nabla f(\boldsymbol x)\Vert}$$

成立时，方向导数最小，从而函数值下降最快。

因为模长对方向没有影响，所以可以下结论：负梯度方向就是最速下降方向。

🔔上述推导中假设了约束中的范数为欧氏范数，而在其他度量空间下时最速下降方向的表达式就会发生变化，不一定再是负梯度方向。例如在$\bf A$ 范数下$\Vert\boldsymbol d\Vert_{\bf A}=(\boldsymbol d^\top\mathbf A\boldsymbol d)^{1/2}$，此时求解上述问题得到的最速方向就不是负梯度方向了。

当取 欧氏范数/$L2$ 范数 时，即最速下降方向是负梯度方向时，最速下降法就成为了梯度下降法。

最速下降算法步骤

最速下降法从初始点$\boldsymbol x^{(1)}$ 出发，每一次都通过当前最速下降方向$\boldsymbol d^{(k)}$ 和步长因子$\lambda_k$ 更新点：$\boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}+\lambda_k\boldsymbol d^{(k)}$ 。当经过多次迭代达到允许误差范围内后，停止计算。

其中的步长因子$\lambda_k$ 并不是一个固定的值，而是通过最优化单变量函数得到的，这个优化过程通常是利用一维搜索。令$\varphi(\lambda)=f(\boldsymbol x^{(k)}+\lambda\boldsymbol d^{(k)})$，则

$$\lambda_k=\arg\min_{\lambda}\varphi(\lambda)$$

整个算法步骤如下：

此处的步长因子在梯度下降法领域被称为学习率，常用$\eta$ 表示。梯度下降法不同的变体对其的设计也有所不同，未必是利用一维搜索实现。详见本站文章：【最优化】梯度下降算法及其变体总结

锯齿现象

容易证明，用最速下降法极小化目标函数时，相邻的两个搜索方向是正交的。说明如下：

在搜索最优步长时，理想情况下为极小化$\varphi(\lambda)$ ，令$\varphi'(\lambda)=0$ 可求得最优$\lambda_k$，此时

$$\begin{aligned} \varphi'(\lambda)&=\nabla f(\boldsymbol x^{(k)}+\lambda\boldsymbol d^{(k)})^\top\cdot\boldsymbol d^{(k)}=0\\\\ \Rightarrow &\nabla f(\boldsymbol x^{(k+1)})^\top\cdot\boldsymbol d^{(k)}=0\quad \color{teal}//\because \boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}+\lambda_k\boldsymbol d^{(k)}\\ \Rightarrow &{\boldsymbol d^{(k+1)}}^\top\cdot\boldsymbol d^{(k)}=0\quad\color{teal}//\because \boldsymbol d^{(k+1)}=-\nabla f(\boldsymbol x^{(k+1)}) \end{aligned}$$

这表明迭代产生的序列$\{\boldsymbol x^{(k)}\}$ 所走的路径呈 “之” 字形（如下图所示），当$\boldsymbol x^{(k)}$ 接近极小点时，步长会很小，锯齿现象（zigzagging phenomenon）尤为明显，这影响了算法是收敛效率。

综上所述，最速下降法反映了目标函数的一种局部性质。局部上看，最速下降方向对当前的搜索有利；但全局上看，向着极小点移动的过程中需要经历不少弯路，收敛速率越来越慢。因此，最速下降法更适用于计算过程的前期迭代部分或间插步骤。

Lemaréchal, C. (2012). “Cauchy and the Gradient Method” (PDF). Doc Math Extra: 251–254.

共轭梯度法

共轭梯度法（conjugate gradient method）最早由 Magnus Hestenes 和 Eduard Stiefel 在1952 年提出。一开始它是用于特定线性方程组——系数矩阵为对称正定阵的线性方程组$\bf Ax=b$ 的数值求解的一种迭代算法。其算法步骤如下：

之后，共轭梯度法被发现还可用于求解无约束优化问题，现已成为了一种重要的最优化方法。

共轭方向

设$A$ 为$n\times n$ 的对称正定矩阵，若$\exists \;\boldsymbol{d}^{(1)},\boldsymbol{d}^{(2)}\in\Bbb R^n$ 使得${\boldsymbol{d}^{(1)}}^\top A\boldsymbol{d}^{(2)}=0$ ，则称它们关于$A$ 共轭。

若$\exists \;\boldsymbol{d}^{(1)},\boldsymbol{d}^{(2)},..,\boldsymbol{d}^{(k)}\in\Bbb R^n$ 两两关于$A$ 共轭，则这组方向关于$A$ 共轭，或称它们为$A$ 的$k$ 个共轭方向。若它们均为非零向量，则它们构成的向量组线性无关。

共轭梯度算法原理

共轭梯度法在最优化问题上的思想是将共轭性与最速下降法结合，利用已知点出的梯度构造一组共轭方向，然后沿着这组方向进行搜索。

共轭梯度法由两个部分组成：（1）一维搜索最优$\lambda_k$ 使$f$ 下降最大；（2）搜索最优$\beta_k$ 使$\boldsymbol d$ 与之前的方向共轭。方向更新公式为：

$$\boldsymbol d^{(k+1)}=-\boldsymbol g_{k+1}+\beta_k\boldsymbol d^{(k)}$$

其中，$\boldsymbol g_k\xlongequal{\rm def}\nabla f(\boldsymbol x^{(k)})$. // 为书写方便，此后的梯度我们也用此变量表示

现在我们先考虑二次凸函数问题：

$$\min\;f(\boldsymbol x)\xlongequal{\rm def}\frac12\boldsymbol x^\top A\boldsymbol x+\boldsymbol b^\top\boldsymbol x+c$$

其中$A$ 为对称正定矩阵。

第一部分

和最速下降法一样，我们先确定步长因子。令$\varphi'(\lambda)=0$ ，即

$$\begin{aligned} \varphi'(\lambda)&=\nabla f(\boldsymbol x^{(k)}+\lambda\boldsymbol d^{(k)})^\top\cdot\boldsymbol d^{(k)}\\ &=[A(\boldsymbol x^{(k)}+\lambda\boldsymbol d^{(k)})+\boldsymbol b]^\top\cdot\boldsymbol d^{(k)}\quad \color{teal}//\because \nabla f=A\boldsymbol x+\boldsymbol b\\ &=[\boldsymbol g_k+\lambda A\boldsymbol d^{(k)}]^\top\cdot\boldsymbol d^{(k)}\\&=0 \end{aligned}$$

可以解得

$$\begin{aligned} \lambda_k=-\frac{\boldsymbol g_k^\top\boldsymbol d^{(k)}}{\boldsymbol d^{(k)\top} A\boldsymbol d^{(k)}} \end{aligned}$$

这个$\lambda_k$ 就是当前的最优步长因子。

第二部分

将方向更新公式$\boldsymbol d^{(k+1)}=-\boldsymbol g_{k+1}+\beta_k\boldsymbol d^{(k)}$ 两边同时左乘${\boldsymbol d^{(k)}}^\top A$ 得：

$${\boldsymbol d^{(k)}}^\top A\boldsymbol d^{(k+1)}=-{\boldsymbol d^{(k)}}^\top A\boldsymbol g_{k+1}+\beta_k{\boldsymbol d^{(k)}}^\top A\boldsymbol d^{(k)}$$

为了让$\boldsymbol d^{(k+1)}$ 和$\boldsymbol d^{(k)}$ 共轭，我们令上式为零，从而解得：

$$\beta_k=-\frac{\boldsymbol d^{(k)\top} A\boldsymbol g_{k+1}}{\boldsymbol d^{(k)\top} A\boldsymbol d^{(k)}}$$

前提条件

可以证明，当初始搜索方向$\boldsymbol d^{(1)}$ 取$-\boldsymbol g_1$ 时（符合最速下降法），利用 FR法 构造的一系列方向关于$A$ 共轭，并且经过有限步迭代必达极小点。
但是，如果初始方向不遵循最速下降法的选取方式，FR法得到的后继方向就未必共轭！！

对一般函数

推广到任意函数$f(\boldsymbol x)$ 的极小化问题时，$\lambda_k$ 不再能显式地求解出来，应与最速下降法类似采用其他一维搜索策略，而使用到的正定矩阵$A$ 则不是固定不变，每一步都使用当前函数的 Hesse 矩阵$\nabla^2f(\boldsymbol x^{(k)})$ 代替。这是将函数的二阶 Taylor 展开式近似为一个二次凸函数推导得来的。

此外，对应一般函数往往共轭梯度法不能直接在有限步求得极小点，所以算法中往往会添加一个重置判断，达到重置条件时重新用当前最速下降方向作为迭代方向，而不再是用共轭梯度法的更新公式。

最后我们得到下列所示的计算步骤（来自：Open Problems in Nonlinear Conjugate Gradient Algorithms for Unconstrained Optimization | Semantic Scholar）：

图中的$\alpha_k$ 即是本文中所说的步长$\lambda_k$，所使用的一维搜索方法为 Wolfe line search；而关于$\beta$ 的选取，不同作者给出了不同的表达式：

对于二次凸函数作为目标函数的情况时，以上表达式均是等价的！但在应对一般函数时，性能有所差别。

牛顿法

牛顿法（Newton’s method）最初由英国科学家艾萨克·牛顿（Isaac Newton）在17世纪提出，首次发表于 1685 年 John Wallis 的《历史与实用代数论》（A Treatise of Algebra both Historical and Practical）中，用于求解多项式方程的根。

1690 年，约瑟夫·拉夫逊（Joseph Raphson）在《普遍等式分析》（Analysis aequationum universalis）一书中给出了更加简化的描述 。虽然他也只将该方法应用于多项式，但他从原始多项式中提取了每个连续修正，从而避免了牛顿原始方法繁琐的重写过程，从而为每个问题推导出一个可重复使用的迭代表达式。因此牛顿法有时也被称为牛顿-拉夫逊方法（Newton-Raphson method）。

最后，1740 年，托马斯·辛普森（Thomas Simpson）将牛顿法描述为一种使用微积分求解一般非线性方程的迭代方法，辛普森还将牛顿法推广到了二元方程组，并指出牛顿法可以通过将梯度设为零来解决优化问题。

牛顿法求根

在数值分析中，牛顿法可以用于求出一元可微函数的零点（方程的根）的数值解。

首先，选取一个接近函数$f(x)$ 零点的初值$x_0$ ，计算对应的函数值$f(x_0)$ 和其切线斜率$f'(x_0)$ 。然后我们给出一条直线，其斜率为$f'(x_0)$ 并且过点$(x_0,f'(0))$ 。记该直线与$x$ 轴的交点坐标为$(x_1,0)$ ，那么其横坐标的值$x_1$ 通常就比$x_0$ 更接近函数的零点$x^*\quad(f(x^*)=0)$ 。

上述过程实际上就等价于求解下列方程：

$$0=(x-x_0)·f'(x_0)+f(x_0)$$

解得

$$x_1=x_0-\frac{f(x_0)}{f'(x_0)}$$

接下来，为了找到零点更精确的值，我们再用同样的方法从$x_1$ 出发寻找$x_2,x_3,...$ 直到满足一定的精度要求。整个过程如下图所示：

事实上，这个方法等价于在$x_0$ 处将原函数近似看作一次函数。这个近似由Taylor展开式忽略高阶项得出：

$$f(x)\approx f(x_0)+f'(x_0)(x-x_0)\xlongequal{\rm def}\varphi(x)$$

此时令$\varphi(x)=0$ 即可解得$x_1$ 。第$k$ 次迭代的更新公式就可以给出：

$$x_k=x_{k-1}-\frac{f(x_{k-1})}{f'(x_{k-1})}$$

优化中的牛顿法

在最优化领域，我们往往需要通过令可微函数的导数为零求其极值点。事实上求“导函数为零的点”就等价于求解“导函数的零点/方程的根”。

所以对于二阶可导的一元函数，我们可以利用牛顿法得到极值点的更新公式。

推广之，对于多元二次可微实函数$f(\boldsymbol x)$ ，通过其二阶 Taylor 展开式进行近似：

$$f(\boldsymbol x)\approx\phi(\boldsymbol x)=f(\boldsymbol x^{(k)})+\nabla f(\boldsymbol x^{(k)})^\top(\boldsymbol x-\boldsymbol x^{(k)})+\frac12(\boldsymbol x-\boldsymbol x^{(k)})^\top\nabla^2 f(\boldsymbol x^{(k)})(\boldsymbol x-\boldsymbol x^{(k)})$$

其中，$\nabla^2 f(\boldsymbol x^{(k)})$ 为 Hesse 矩阵。

令$\nabla\phi(\boldsymbol x)=\bf0$ 即可求得平稳点的近似解。可得迭代公式：

$$\boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}-\nabla^2 f(\boldsymbol x^{(k)})^{-1}\nabla f(\boldsymbol x^{(k)})$$

整个算法流程如下图所示：

二次终止性

对于$n$ 元正定二次凸函数$f(\boldsymbol x)=\frac12\boldsymbol x^\top A\boldsymbol x+\boldsymbol b^\top\boldsymbol x+c$ ，我们可以利用$\nabla f(\boldsymbol x)=\bf0$ 求出其极小点的解析解$\boldsymbol x^*=-A^{-1}\boldsymbol b$ 。

而当我们利用牛顿法来求解该问题时，任取一个初点$\boldsymbol x^{(1)}$ ，计算其后继点：

$$\begin{aligned} \boldsymbol x^{(2)}&=\boldsymbol x^{(1)}-\nabla^2 f(\boldsymbol x^{(1)})^{-1}\nabla f(\boldsymbol x^{(1)})\\ &=\boldsymbol x^{(1)}-A^{-1}(A\boldsymbol x^{(1)}+\boldsymbol b)\\ &=-A^{-1}\boldsymbol b \end{aligned}$$

可以发现，正好就是最优解。这说明牛顿迭代法具有二次终止性（也称二次收敛性）。

$\bf Define:$ 对于$n$ 元正定二次凸函数$f(\boldsymbol x)=\frac12\boldsymbol x^\top A\boldsymbol x+\boldsymbol b^\top\boldsymbol x+c$ 求极小点的问题，若某个迭代算法能从任意点出发，经过有限步迭代达到最优解，则称该算法具有二次终止性。

我们可以下结论：最速下降法不具备二次终止性，而共轭梯度法具有二次终止性。

此外，牛顿法在极小点附近具有很快的收敛速度，避免了最速下降法的锯齿现象。上述三种算法的收敛可视化如下图所示：

阻尼牛顿法

在使用牛顿迭代法求解最优化问题时，有一点值得注意：当初始点远离极小点时，牛顿法可能不收敛。

记牛顿方向$\boldsymbol d=-\nabla^2 f(\boldsymbol x^{(k)})^{-1}\nabla f(\boldsymbol x^{(k)})$ ，此方向并不一定是下降方向。也就是说利用牛顿法迭代，函数值甚至还有上升的可能。另外，即是它当前是下降方向，直接使用$\boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}+\boldsymbol d^{(k)}$ 来更新，得到的新点也未必是这个方向下最佳的选择。

因此，阻尼牛顿法（Damped Newton’s Method）提了出来，和最速下降法和共轭梯度法一样引入了可一维搜索的步长因子$\lambda_k$ 来对牛顿方向进行限制：即使牛顿方向使函数值上升，优化过程中也会因为一维搜索这一步迫使$\lambda_k\to0$ 来消减影响。

修正牛顿法

虽然阻尼牛顿法解决了原始牛顿法的一些问题，但是它们仍然拥有共同的缺点：

1. 可能出现 Hesse 矩阵 是奇异的情况，因其不可逆使更新公式失效；
2. Hesse 矩阵未必正定，此时牛顿方向未必下降，阻尼牛顿算法也会失效（$\lambda=0$ 时$\boldsymbol x$ 得不到更新）。

为解决上述问题，修正牛顿法（Modified Newton’s Method）提出构造一个对称正定矩阵$B_k$ 来取代 Hesse 矩阵$\nabla^2 f(\boldsymbol x^{(k)})$ .

构造$B_k$ 的其中一种方法是令

$$B_k=\nabla^2 f(\boldsymbol x^{(k)})+\epsilon_k\boldsymbol I$$

$\boldsymbol I$ 为$n$ 阶单位阵，$\epsilon_k$ 为一个适当的正数。 假设 Hesee 矩阵的特征值为$\alpha_k$，那么合适的$\epsilon_k$ 可以保证$B_k$ 的特征值$\alpha_k+\epsilon_k$ 均大于零，从而保证其正定性。

最后我们给出修正牛顿法的算法描述：

拟牛顿法

牛顿法最突出的特点就是收敛速度很快，但是它需要计算二阶偏导数，并且 Hesse 矩阵可能非正定。修正牛顿法虽然解决了 Hesse 矩阵的非正定问题，但是对目标函数的二阶可微性、计算二阶偏导数的计算代价 仍然有要求。为此，一类提出直接使用与二阶导数无关的矩阵来近似牛顿法中的$\nabla^2 f(\boldsymbol x^{(k)})^{-1}$ 的算法诞生了，这类算法我们就称为拟牛顿法。

拟牛顿条件

为构造$\nabla^2 f(\boldsymbol x^{(k)})^{-1}$ 的近似矩阵$G_k$ ，为了不使用到二阶导，我们需要考察$\nabla^2 f(\boldsymbol x^{(k)})^{-1}$ 与一阶导数的关系。

设在第$k$ 次迭代后得到了点$\boldsymbol x^{(k+1)}$ ，此时目标函数在此处的二阶 Taylor 近似如下：

$$f(\boldsymbol x)\approx f(\boldsymbol x^{(k+1)})+\nabla f(\boldsymbol x^{(k+1)})^\top(\boldsymbol x-\boldsymbol x^{(k+1)})+\frac12(\boldsymbol x-\boldsymbol x^{(k+1)})^\top\nabla^2 f(\boldsymbol x^{(k+1)})(\boldsymbol x-\boldsymbol x^{(k+1)})$$

进而梯度$\nabla f(\boldsymbol x)\approx \nabla f(\boldsymbol x^{(k+1)})+\nabla ^2f(\boldsymbol x^{(k+1)})(\boldsymbol x-\boldsymbol x^{(k+1)})$。

接下来取$\boldsymbol x=\boldsymbol x^{(k)}$ 带入上式，令$\boldsymbol s^{(k)}=\boldsymbol x^{(k+1)}-\boldsymbol x^{(k)}, \boldsymbol y^{(k)}=\nabla f(\boldsymbol x^{(k+1)})-\nabla f(\boldsymbol x^{(k)})$，得到

$$\boldsymbol y^{(k)}\approx\nabla^2 f(\boldsymbol x^{(k+1)})\boldsymbol s^{(k)}$$

当 Hesse 可逆时，$\boldsymbol s^{(k)}\approx\nabla^2 f(\boldsymbol x^{(k+1)})^{-1}\boldsymbol y^{(k)}$.

于是，如果我们需要找到一个矩阵$G_{k+1}$ 来近似 Hesse 矩阵的逆，相当于它需要满足下列等式：

$$\boldsymbol s^{(k)}=G_{k+1}\boldsymbol y^{(k)}$$

这个公式也被叫做拟牛顿条件(Secant equation).

构造近似矩阵

满足拟牛顿条件的 矩阵$G_{k}$ 应当是正定的以保证拟牛顿方向是下降方向。那么构造的一般策略是取$G_1$ 为任意$n$ 阶对称正定矩阵（通常取单位阵$I$），然后通过矫正矩阵$\Delta G_k$ 来修正$G_k$ 得到$G_{k+1}$：

$$G_{k+1}=G_k+\Delta G_k$$

接下来就是如何设计矫正矩阵了。不同的算法有着不同的设计方式。

秩一矫正

我们知道一个$n$ 阶方阵可以由一个$n$ 维列向量经过乘积得到，比如设：

$$\Delta G_k=\alpha_k\boldsymbol z^{(k)}{\boldsymbol z^{(k)}}^\top$$

其中，$\alpha_k$ 为一个常数，$\boldsymbol z^{(k)}$ 为$n$ 维向量。这样的矩阵正好是秩一矩阵。

将这样的矫正矩阵带入拟牛顿条件，可以得到

$$\Delta G_k=\frac{(\boldsymbol s^{(k)}-G_k\boldsymbol y^{(k)})(\boldsymbol s^{(k)}-G_k\boldsymbol y^{(k)})^\top}{\boldsymbol y^{(k)\top}(\boldsymbol s^{(k)}-G_k\boldsymbol y^{(k)})}$$

DFP算法

DFP算法是拟牛顿法中一种著名的算法，由 Davidon 首先提出，后经 Fletcher 和 Powell 改进，DFP的命名取自他们三人的名字 Davidon-Fletcher-Powell。 该算法也被称为 变尺度法。

DFP 设计矫正矩阵如下：

$$\Delta G_k=\alpha_k\boldsymbol u^{(k)}{\boldsymbol u^{(k)}}^\top+\beta_k\boldsymbol v^{(k)}{\boldsymbol v^{(k)}}^\top$$

并且，取$\boldsymbol u^{(k)}=\boldsymbol s^{(k)},\;\boldsymbol v^{(k)}=G_k\boldsymbol y^{(k)}$ 以及限制$\alpha_k=1/{\boldsymbol s^{(k)}}^\top\boldsymbol y^{(k)}$ 和$\beta_k=1/{\boldsymbol y^{(k)}}^\top G_k\boldsymbol y^{(k)}$.

最后得到

$$\Delta G_k=\frac{\boldsymbol s^{(k)}\boldsymbol s^{(k)\top}}{\boldsymbol s^{(k)\top}\boldsymbol y^{(k)}}-\frac{G_k\boldsymbol s^{(k)}(G_k\boldsymbol s^{(k)})^\top}{\boldsymbol y^{(k)\top} G_k\boldsymbol y^{(k)}}$$

Broyden族

前面我们提到的$G_k$ 都是用来近似$\nabla^2 f(\boldsymbol x^{(k)})^{-1}$ (有逆) 的。事实上我们还可以用矩阵$B_k$ 来近似$\nabla^2 f(\boldsymbol x^{(k)})$ (无逆)。

类似的方法我们有$B_{k+1}=B_k+\Delta B_k$ ，并且$\boldsymbol y^{(k)}=B_{k+1}\boldsymbol s^{(k)}$ 。

显然，要找到一个$\Delta B_k$ 的设计方法使上式成立，我们可以化用 DFP 的设计方法，只是将变量稍作修改而已。也确实存在这样的算法：Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm。

有研究表明，BFGS 和 DFP 的加权线性组合仍然满足拟牛顿公式，权重选取不同得到的修正公式就不同，全体的修正公式被称为Broyden族。

此外还有其他的设计方法：

信赖域方法

待更

附赠：一维搜索

待更