Numerical Optimization Ch9. Derivative-Free Optimization

第九章: 无导数优化

许多实际的问题给出来时 并不自带(或不易求得)导数信息. 原则上来说, 这类问题可以用 第八章的 有限差分法近似梯度(甚至是Hessian矩阵), 再把这些近似的导数用进更早一些章节的算法和框架中. 不过尽管有限差分法在某些应用上很有成效, 但它却 不能被当成一种无导数优化的普适方法. 原因有二:

1. 可能需要计算过多的函数值;
2. 数据中噪声(本章中我们定义噪声为计算函数时不准确度的度量)的存在, 使得这一方法不再可靠.

因此, 现已研究出多种不需要近似梯度的算法. 它们仅需要采样点集上的函数值以通过其他的方式得到新的迭代点.
无导数优化(Derivative-free optimization, DFO)算法区别于它们如何使用采样的函数值决定新的迭代点. 其中一类通过构建目标函数的线性或二次模型, 并在一个信赖域中极小化之求得. 我们将特别关注这些基于模型的算法, 这是因为它们与我们之前讨论的无约束极小化算法联系紧密. 其他的广为使用的DFO算法包括Nelder-Mead的单纯形反射法, 模式搜索法, 共轭方向法和模拟退火法. 本章中我们仅简要地介绍这些方法(除了模拟退火法, 这是因为它是一种非确定性的方法, 也即随机方法).
无导数优化算法并不如之前的那些基于梯度的算法发展得好. 当前的算法仅对小规模问题有效. 尽管大多数DFO算法已被推广至能够处理带简单约束的问题, 例如边界约束, 但如何处理带一般类型约束的问题尚在研究中. 这里我们也仅限讨论无约束优化问题$$\underset{x\in {\mathbb{R}}^n}{\min }f(x).$$
无导数优化算法通常(往往效果喜忧参半)用于极小化不可导的函数或者是尝试定位一个函数的全局极小点. 由于我们不牵涉对非光滑优化或全局优化的讨论, 我们的讲述将仅限于光滑问题, 即目标函数$f$具有连续导数的情形. 我们会在第1,6节讨论噪声的影响.

1. 有限差分与噪声

一种显然的DFO算法是用有限差分法近似估计梯度, 之后再使用基于梯度的算法. 这方法有时会很有效, 并且我们也应当首先考虑这样的方法, 但有限差分估计在目标函数包含噪声的时候可能会变得不够准确. 我们在本节中尝试量化噪声的影响.

函数计算过程的噪声有许多的来源.

• 函数$f(x)$依赖于随机模拟, 则在计算函数时就会带上随机误差;
• $f$来源于求解某个微分方程或其他的复杂数值程序, 则在计算时所涉及的容忍限也会给$f$的值带来误差;
• …

于是, 目标函数$f$具有形式$$f(x)=h(x)+\varphi (x),$$其中$h$为一光滑函数, $\varphi$表示噪声. 注意我们把$\varphi$写成$x$的函数但实际上是不必要的. 例如, 若$f$的计算依赖与随机模拟, 则$\varphi$的值在每次计算中都会不一样, 即使是在同样的$x$上. 不过上面的表示在阐释噪声带来的负面影响以及发展无导数优化的算法时是非常有用的.

给定$ϵ$, 则$\nabla f$在$x$处的中心差分近似为$${\nabla }_{ϵ}f(x)={\left[\frac{f(x+ϵe_i)-f(x-ϵe_i)}{2ϵ}\right]}_{i=1,2,\dots ,n},$$其中$e_i$为第$i$个单位向量. 我们想要联系${\nabla }_{ϵ}f(x)$与潜在的光滑函数$h(x)$. 定义噪声水平$\eta$为$\varphi$在以$x$为中心$2ϵ$为边长的方形闭区域中的最大值, 即$$\eta (x;ϵ)=\underset{\Vert z-x{\Vert }_{\infty }\le ϵ}{\sup }|\varphi (z)|.$$将中心差分公式代入可得以下结果.

引理1 设${\nabla }^{2}h$在方形邻域 { z ∣ ∥ z − x ∥ ∞ ≤ ϵ } \{z|\Vert z-x\Vert_{\infty}\le\epsilon\} {z∣∥z−x∥∞​≤ϵ}内Lipschitz连续, 其中Lipschitz常数为$L_h$. 则我们有$$\Vert {\nabla }_{ϵ}f(x)-\nabla h(x){\Vert }_{\infty }\le L_h{ϵ}^2+\frac{\eta (x;ϵ)}{ϵ}.$$

因此中心差分的误差来源有两个:

• 有限差分本身的误差$O({ϵ}^2)$; 以及
• 噪声带来的误差$\eta (x;ϵ)/ϵ$.

若噪声极大地影响了差分区域, 我们对${\nabla }_{ϵ}f(x)$就不能抱有太大的幻想了. 这时候其实若$-{\nabla }_{ϵ}f(x)$是下降方向就已经够幸运了.
除了有限差分法中使用的当前点附近的较小区域的函数值, 我们其实可以将选取的点分的更开并使用它们构建目标函数的模型. 这一方法(我们将在下一节和第6节讨论)可能对于噪声更加强健.

2. 基于模型的方法

前面几章中一些无约束优化的高效算法往往是先对目标函数$f$构建一个二次模型, 再尝试极小化之以获取搜索方向. 这一模型基于当前迭代点的函数值和导数值信息. 当导数无法获取时, 我们可以定义模型$m_k$为目标函数$f$在适当选定的采样点上函数值所定义的插值二次函数. 这样的模型通常是非凸的, 因此本章中基于模型的算法将使用信赖域的框架计算搜索方向.

假设在当前迭代点$x_k$处我们有采样点集 Y = { y 1 , y 2 , … , y q } , y i ∈ R n , i = 1 , 2 , … , q Y=\{y^1,y^2,\ldots,y^q\},y^i\in\mathbb{R}^n,i=1,2,\ldots,q Y={y1,y2,…,yq},yi∈Rn,i=1,2,…,q. 假设$x_k$为集合中的一个元素且在$Y$中没有任何元素比$x_k$有更小的函数值. 我们希望建立二次模型$$m_k(x_k+p)=c+g^{T}p+\frac{1}{2}{p}^{T}Gp.$$我们不能直接定义$g=\nabla f(x_k),G={\nabla }^{2}f(x_k)$, 而是通过以下插值条件间接决定(包括$c$):$$m_k(y^l)=f(y^l),l=1,2,\dots ,q.$$由于模型中有$1+n+\frac{1}{2}n(n+1)=\frac{1}{2}(n+1)(n+2)$个未知量, 因此满足插值条件的$m_k$唯一仅当$$q=\frac{1}{2}(n+1)(n+2).$$此时插值条件可写作未知量的线性方程组. 若选取的插值点$y^1,y^2,\dots ,y^q$使得线性系统的系数矩阵非奇异, 则理论上$m_k$是唯一确定的.

一旦有了$m_k$, 我们就可通过求解信赖域子问题$$\underset{p}{\min }{m}_k(x_k+p),\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert p{\Vert }_2\le \Delta$$计算搜索方向. 求解方法可见第四章. 若$x_k+p$在目标函数上取得了充分的下降, 则新的迭代点就取为$x_{k+1}=x_k+p$, 并更新信赖域半径, 进行新一步迭代. 否则拒绝这一搜索方向, 此时需要改进插值点集$Y$或者缩小信赖域.

为减小算法的计算耗费, 每步迭代中我们将更新$m_k$而不是从头重新计算. 实际上, 我们会选取二次多项式空间中的一组便利的基, 比如Lagrange多项式和Newton多项式. 这些基的性质即可以用来度量采样点集$Y$是否合适, 也可以在需要的时候对其作出一些改变. 将所有这些问题考虑进去并有效解决的完整算法要比第六章中讨论的拟牛顿法复杂得多. 因此, 对于基于模型的DFO方法, 我们仅提供大致的框架.

正如一般的信赖域算法, 步长接收和信赖域更新的策略均基于函数值实际的减少与模型预测的减少的比值, 即$$\rho =\frac{f(x_k)-f(x_k^{+})}{m_k(x_k)-m_k(x_k^{+})},$$其中$x_k^{+}$表示试探点. 本节中整数$q$定义为$\frac{1}{2}(n+1)(n+2)$.

算法1 (Model-Based Derivative-Free Method)
Choose an interpolation set Y = { y 1 , y 2 , … , y q } Y=\{y^1,y^2,\ldots,y^q\} Y={y1,y2,…,yq} such that the linear system is nonsingular, and select $x_0$ as a point in this set such that $f(x_0)\le f(y^i)$ for all $y^i\in Y$. Choose an initial trust region radius ${\Delta }_0$, a constant $eta\in (0,1)$, and set $k\leftarrow 0$.

repeat until a convergence test is satisfied:
\quad\quad Form the quadratic model $m_k(x_k+p)$ that satisfies the interpolation conditions;
\quad\quad Compute a step $p$ by approximately solving subproblem;
\quad\quad Define the trial point as $x_k^{+}=x_k+p$;
\quad\quad Compute the ratio $\rho$;
\quad\quad if $\rho \ge \eta$
\quad\quad\quad\quad Replace an element of $Y$ by $x_k^{+}$;
\quad\quad\quad\quad Choose ${\Delta }_{k+1}\ge {\Delta }_k$;
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k^{+}$;
\quad\quad\quad\quad Set $k\leftarrow k+1$ and go to the next iterations;
\quad\quad else if the set $Y$ need not be improved
\quad\quad\quad\quad Choose ${\Delta }_{k+1}<{\Delta }_k$;
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k$;
\quad\quad\quad\quad Set $k\leftarrow k+1$ and go to the next iteration;
\quad\quad end (if)
\quad\quad Invoke a geometry-improving procedure to update $Y$:
\quad\quad\quad\quad at least one of the points in $Y$ is replaced by some other point, with the goal of improving the interpolation condition;
\quad\quad Set ${\Delta }_{k+1}\leftarrow {\Delta }_k$;
\quad\quad Choose $\hat{x}$ as an element in $Y$ with lowest function value;
\quad\quad Set $x_k^{+}\leftarrow \hat{x}$ and recompute $\rho$;
\quad\quad if $\rho \ge \eta$
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k^{+}$;
\quad\quad else
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k$;
\quad\quad end (if)
\quad\quad Set $k\leftarrow k+1$;
end (repeat)

算法说明:

1. $\rho \ge \eta$的情形表示我们在价值函数上获得了充分的下降. 此时我们将接受$x_k^{+}$为新的迭代点并加入$Y$, 再从$Y$中移除一个元素. 而当下降不足$(\rho <\eta )$, 此时就有两种诱因: 插值点集$Y$选取不当或信赖域过大. 前者可能是因为迭代点都差不多在${\mathbb{R}}^n$的一个不包含解的低维曲面上. 此时算法可能会收敛到这个子集中的极小点. 这种情形可通过监测线性系统(具体说, 是插值条件的满足情况, 后面会说明)的条件数发现. 若条件数过高, 我们就要改进$Y$, 通常是以一个新元素替代$Y$中现有的一个元素使得插值系统尽可能远离奇异. 若$Y$的选取适当(即条件数适宜), 我们仅缩小信赖域半径$\Delta$, 这就与第四章相同了.
2. 一种较好的初始$Y$是${\mathbb{R}}^n$中单纯形的顶点和边上的中点.
3. 二次模型的使用限制了可以求解的问题的规模. 光是为启动算法就需要达$O(n^2)$次函数值的计算, 这对较小的$n$(比如$n=50$)也是很繁重的. 另外, 迭代过程的计算量也很大. 即使每步对$m_k$更新而不是重新计算, 用于构建$m_k$和计算步长的计算量也将达$O(n^4)$.
   为缓解这些困难, 我们可以线性模型代替二次模型, 即设$G=O$. 这样的模型只有$n+1$个参数, 从而只需$Y$有$n+1$个插值点, 此时每步迭代的耗费为$O(n^3)$. 而算法1对于线性模型也有较好的兼容度(仅需较小修改), 但这样一来就牺牲了收敛速度, 这是因为线性模型难以捕捉问题的曲率信息. 因此, 有些基于模型的算法以$n+1$个初始点出发, 用线性模型计算步长, 但当$q=\frac{1}{2}(n+1)(n+2)$个函数值变得可计算时便转而使用二次模型.

2.1 插值与多项式基

下面我们详细介绍如何利用插值构建目标函数的模型.

1. 线性模型.
   我们先从线性模型开始:$$m_k(x_k+p)=f(x_k)+g^{T}p.$$为确定向量$g\in {\mathbb{R}}^n$, 我们施加插值条件$m_k(y^l)=f(y^l),l=1,2,\dots ,n$, 也即$$(s^l{)}^{T}g=f(y^l)-f(x_k),l=1,2,\dots ,n,$$其中$$s^l=y^l-x_k,l=1,2,\dots ,n.$$于是便有$g$满足的线性方程组, 其系数矩阵的行为$(s^l{)}^T$. 因此模型被唯一确定当且仅当插值点 { y 1 , y 2 , … , y n } \{y^1,y^2,\ldots,y^n\} {y1,y2,…,yn}使得向量组 { s l : l = 1 , 2 , … , n } \{s^l:l=1,2,\ldots,n\} {sl:l=1,2,…,n}是线性无关的. 若这一条件成立, 则我们称由$x_k,y^1,y^2,\dots ,y^n$构成的单纯形是非退化的.
2. 二次模型.
   接下来考虑构建二次模型. 重写模型如下:$$\begin{array}{rl}{m}_k(x_k+p)& =f(x_k)+g^{T}p+\sum _{i<j}{G}_{ij}{p}_i{p}_j+\frac{1}{2}\sum _i{G}_{ii}{p}_i^2\\ & \stackrel{def}{=}f(x_k)+{\hat{g}}^T\hat{p},\end{array}$$其中 g ^ ≡ ( g T , { G i j } i &lt; j , { 1 2 G i i } ) T , p ^ ≡ ( p T , { p i p j } i &lt; j , { 1 2 p i 2 } ) T . \hat{g}\equiv\left(g^T,\{G_{ij}\}_{i&lt;j},\left\{\frac{1}{\sqrt{2}}G_{ii}\right\}\right)^T,\quad \hat{p}\equiv\left(p^T,\{p_ip_j\}_{i&lt;j},\left\{\frac{1}{\sqrt{2}}p_i^2\right\}\right)^T. g^​≡(gT,{Gij​}i<j​,{2 ​1​Gii​})T,p^​≡(pT,{pi​pj​}i<j​,{2 ​1​pi2​})T.此时二次模型与线性模型具有同样的形式, 因此亦可用同样的方法确定未知的系数$\hat{g}$.

多元二次函数可以用多种方式表示. 单项式基的表示让我们能够通过设置$G$中一些元素为0来表现已知的Hessian结构. 而其他的基则在规避线性系统奇异性时显得更加的便利.
我们以 { ϕ i ( ⋅ ) } i = 1 q \{\phi_i(\cdot)\}_{i=1}^q {ϕi​(⋅)}i=1q​表示$n$维二次函数空间中的一组基. 因此二次模型可表示成$$m_k(x)=\sum _{i=1}^q{\alpha }_i{\varphi }_i(x).$$当行列式$$\delta (Y)\stackrel{def}{=}\det \left(\begin{array}{ccc}{\varphi }_1(y^1)& \cdots & {\varphi }_1(y^q)\\ ⋮& & ⋮\\ {\varphi }_q(y^1)& \cdots & {\varphi }_q(y^q)\end{array}\right)$$不为0时, 插值点集 Y = { y 1 , y 2 , … , y q } Y=\{y^1,y^2,\ldots,y^q\} Y={y1,y2,…,yq}唯一决定系数${\alpha }_i$.
随着基于模型的算法行进, 行列式$\delta (Y)$可能会接近0, 从而引发数值上的不稳定. 因此一些算法包含了保持插值点处于合适位置的机制. 下面我们介绍其中一种.

2.2 更新插值点集

与其等着行列式$\delta (Y)$小于某个阈值, 我们可以在一个试探点不能提供$f$的充分下降时启动一种几何改进的程序, 其目的在于替换插值点集中的某个点使得行列式获得量级上的增加. 为说明, 我们先讨论$\delta (Y)$的性质, 期间使用Lagrange函数的术语.
对每个$y\in Y$, 我们定义Lagrange函数$L(\cdot ,y)$为次数至多为2的多项式, 使其满足$L(y,y)=1,L(\hat{y},y)=0,\hat{y}̸=y,\hat{y}\in Y$. 假设$Y$的更新是通过以$y_{+}$替换$y_{-}$完成的, 从而得到新集$Y^{+}$. 我们可以证明(在适当的规范化与一定的条件之下)有$$|\delta (Y^{+})|\le |L(y_{+},y_{-})||\delta (Y)|.$$算法1可充分利用这一不等式来更新插值集合:

1. 试探点$x^{+}$提供了目标函数的充分下降的情形$(\rho \ge \eta )$. 我们将$x^{+}$放入$Y$中并从中移出$y_{-}$. 受上面不等式的启发, 我们选取移出的点$y_{-}$为$$y_{-}=\arg \underset{y\in Y}{\max }|L(x^{+},y)|.$$
2. $x^{+}$不能提供目标函数的充分下降$(\rho <\eta )$, 我们就有两种情形需要考虑.
   • $Y$需要改进. 我们定义$Y$在当前迭代点$x_k$处是合适的, 若对于所有$y^i\in Y:\Vert x_k-y^i\Vert \le \Delta$, 我们以任一信赖域中的点$y$代替其中的某一插值点$y^i$, $|\delta (Y)|$不会翻倍. 若$Y$不合适, 则启动几何改进机制. 我们选取$y_{-}\in Y$并以另一点$y^{+}$代替之. 其中$y^{+}$的选取基于对行列式的改进. 对每个点$y^i\in Y$, 我们定义其潜在的替换点$y_r^i$为$$y_r^i=\arg \underset{\Vert y-x_k\Vert \le \Delta }{\max }|L(y,y^i)|.$$此时选取$y_{-}$为使得$|L(y_r^i,y^i)|$最大的$y^i\in Y$, $y^{+}$就为对应的$y_r^i$.
   • $Y$不需改进. 即$Y$是合适的, 此时需缩小信赖域的半径并重新开始迭代.

高效地实施这些规则并不直接, 而且注意我们并没有考虑所有可能的问题. 改进插值点集的策略仍然有待研究.

2.3 基于最小改变量更新的算法

我们下面考虑一种可视为是拟牛顿法扩展的方法. 此方法使用二次模型而每次迭代仅需$O(n^3)$的运算量, 这要比之前$O(n^4)$量级的计算量要小得多. 为此, 这一方法只利用$O(n)$个点上的插值条件, 而剩下的自由度则以最小改变量为约束消去. 最小改变量性质是拟牛顿法的关键. 当然这里还有其他的要求, 比如模型拟合最近两个迭代点处的梯度$\nabla f$. 下面我们介绍结合最小该变量以及函数值插值的算法.

在算法的第$k$步迭代, 我们将建立新的二次模型$m_{k+1}$. 模型中的系数$f_{k+1},g_{k+1},G_{k+1}$则通过求解下面的问题获取:$$\begin{array}{rl}{\min }_{f,g,G}& \Vert G-G_k{\Vert }_F^2\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& G\mathrm{s}\mathrm{y}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{c}\\ & m(y^l)=f(y^l)l=1,2,\dots ,\hat{q},\end{array}$$其中$G_k$为之前模型$m_k$中的Hessian, $\hat{q}$为依赖于$n$的整数. 我们可以证明整数$\hat{q}$必须大于$n+1$, 以保证$G_{k+1}$不等于$G_k$. 实际中合适的取法为$\hat{q}=2n+1$, 此时插值点的数量大约为线性模型中的两倍.

上述优化问题是一个带等式约束的二次规划, 其KKT条件可用等式加以叙述. 一旦$m_{k+1}$确定下来, 我们就求解信赖域子问题获取新步长. 这一方法同样需要保证插值点集$Y$的几何形状是合适的. 因此我们将再加上两个极小的约束:

1. $Y$应当使得$m(y^l)=f(y^l)$的成立不依赖于右端;
2. $y^i$不应全部位于一张超平面上.

若这两个条件成立, 则上述优化问题就有唯一解.

基于以上的算法与算法1类似: 它也需要产生新迭代点以及改进$Y$的几何形状. 它的优势之一为, 它仅需$O(n)$个插值点就能产生较好的搜索方向. 实际上, 此法往往只需远远少于$\frac{1}{2}(n+1)(n+2)$次函数值计算便可逼近解. 但由于它仅仅是最近提出, 从而没有充分的数值实验能够彰显它潜在的威力.

3. 坐标和模式搜索方法

坐标搜索和模式搜索方法并不从函数值中显式地构建$f$的模型, 而是从当前迭代点出发沿着特定的方向寻找具有更小函数值的点. 若找到了这样的点, 它们就步进并重复之前的过程, 当然此前可能会对之前的搜索方向做一些改动. 当没有满意的点存在时, 算法就可能会调整沿着当前搜索方向前进的步长, 或直接产生新的搜索方向.
我们先介绍一种实际操作中常用的简单方法. 而后再考虑推广的、更加高效的以及具有更强大的理论性质的方法.

3.1 坐标搜索方法

坐标搜索方法(也称为坐标下降法或交替变量法)循环遍历$n$个坐标方向$e_1,e_2,\dots ,e_n$, 在每个方向上均以线搜索获取新的迭代点. 具体说来, 第一次迭代中, 我们仅改变$x$的第一个分量$x_1$并寻找极小(或至少减小)目标函数的这一分量的新值. 下一步迭代中, 则对第二个分量$x_2$进行同样的操作. 在$n$步迭代后, 我们又回到第一个变量重复遍历. 尽管简单直观, 但此法在实际操作中却可能十分低效, 就比如下图中表示的, 我们将坐标搜索方法应用于带两个变量的二次函数上.

图中显示, 经几步迭代后, $x_1,x_2$向解的行进都明显变慢.

一般说来, 坐标搜索法甚至可以无限地迭代而从不逼近目标函数的稳定点, 这在精确线搜索下依然如此(相比较而言, 在一定条件下最速下降法能产生迭代序列 { x k } \{x_k\} {xk​}使得$\Vert \nabla f_k\Vert \to 0$). 事实上, 循环遍历任何一组线性无关的搜索方向均不能保证全局收敛性. 理论上, 这是由于最速下降方向$-\nabla f_k$变得与搜索方向垂直, 此时Zoutendijk条件在$\nabla f_k$并不趋于0时就满足了.

就算坐标搜索收敛于一解, 一般来说它也要比最速下降法慢得多, 并且这一差距随着变量数的增加而越变越大. 不过由于坐标搜索不需要计算梯度$\nabla f_k$, 且它在变量耦合程度不高时的效果还不错, 因此它仍是优化目标函数的有效方法.

至今已有坐标搜索算法的变体, 其中的一些(被证明)具有全局收敛性. 一种简单的变体是"后退-前进"方法, 其中我们依次沿着以下方向进行搜索:$$e_1,e_2,\dots ,e_{n-1},e_n,e_{n-1},\dots ,e_2,e_1,e_2,\dots ,(\mathrm{r}\mathrm{e}\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{t}).$$而还有一种方法则是受上面图中的启发, 先进行一轮坐标下降, 再沿着连接上一轮中第一个和(指向)最后一个点的方向搜索. 许多算法(如Hooke和Jeeves)均基于这样的想法.
下面的模式搜索算法则推广了坐标搜索: 它在每步迭代中允许考虑和使用更大的搜索方向集.

3.2 模式搜索方法

我们考虑这样的模式搜索算法: 在每个迭代点上选取特定的搜索方向并在每个方向上计算特定步长下的$f$. 这些候选点就构成了一个当前迭代点附近的"框架"或"模式". 若搜寻过程中找到了显著较低的函数值, 则就将对应的点作为新的迭代点, 从而框架的中心也移到这个点上. 事实上不论移动与否, 框架都会以一定的方式改变(如可能改变搜索方向的集合, 亦或者会增大或减小步长), 而后重复之前的过程. 对于特定的模式搜索方法, 我们是有全局收敛性的——特别地, 是当算法有稳定的聚点时.

噪声或其他形式的不精确性(如函数值的计算)可能会影响模式搜索算法的效果, 从而影响收敛. 我们也可以用简单的例子说明, 非光滑性亦可能导致较差的算法表现, 尽管在非光滑问题上往往能得到令人满意的全局收敛.

为引入模式搜索算法, 我们做一些符号说明. 对于当前迭代点$x_k$, 我们定义${\mathcal{D}}_k$为可能的搜索方向集, ${\gamma }_k$为线搜索的参数. 从而框架由点$x_k+{\gamma }_k{p}_k,p_k\in {\mathcal{D}}_k$组成. 当框架中其中一点带来了函数值的显著下降, 我们便移步且亦可能增大${\gamma }_k$, 从而扩大下一次迭代的框架. 若框架中无点达到要求, 就减小${\gamma }_k$(缩小框架), 并设$x_{k+1}=x_k$, 之后重复. 注意不论在哪种情形, 我们都在下一步迭代开始之前改变方向集${\mathcal{D}}_k$, 其中可能需要受限于特定的约束.

算法2 (Pattern-Search)
Given convergence tolerance ${\gamma }_{\mathrm{t}\mathrm{o}\mathrm{l}}$, contraction parameter ${\theta }_{\max }$, sufficient decrease function $\rho :[0,\infty )\to \mathbb{R}$ with $\rho (t)$ an increasing function of $t$ and $\rho (t)/t\to 0$ as $t\downarrow 0$;
Choose initial point $x_0$, initial step length ${\gamma }_0>{\gamma }_{\mathrm{t}\mathrm{o}\mathrm{l}}$, initial direction set ${\mathcal{D}}_0$;
for $k=1,2,\dots$
\quad\quad if ${\gamma }_k\le {\gamma }_{\mathrm{t}\mathrm{o}\mathrm{l}}$
\quad\quad\quad\quad stop;
\quad\quad if $f(x_k+{\gamma }_k{p}_k)<f(x_k)-\rho ({\gamma }_k)$ for some $p_k\in {\mathcal{D}}_k$
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k+{\gamma }_k{p}_k$ for some such $p_k$;
\quad\quad\quad\quad Set ${\gamma }_{k+1}\leftarrow {\varphi }_k{\gamma }_k$for some ${\varphi }_k\ge 1$; (increase step length)
\quad\quad else
\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k$;
\quad\quad\quad\quad Set ${\gamma }_{k+1}\leftarrow {\theta }_k{\gamma }_k$, where$0<{\theta }_k\le {\theta }_{\max }<1$;
\quad\quad end (if)
end (for)

算法说明:

1. 方向集${\mathcal{D}}_k$的选择对于算法的表现以及理论性质的证明是至关重要的. 我们需要对其提两个条件:

   1. 只要$\nabla f(x_k)̸=0$(即$x_k$不是稳定点), 这一集合中至少有一个方向能够给$f$带来下降.
   2. 集合${\mathcal{D}}_k$中向量的长度差不多相同, 从而框架的直径可用步长参数${\gamma }_k$表示. 因此, 要求$${\beta }_{\min }\le \Vert p\Vert \le {\beta }_{\max },p\in {\mathcal{D}}_k,$$其中${\beta }_{\min },{\beta }_{\max }$为独立于$k$的正常数.

   下面具体说明. 我们同第三章中一样定义负梯度与搜索方向的夹角$$\cos \theta =\frac{-\nabla f_k^{T}p}{\Vert \nabla f_k\Vert \Vert p\Vert }.$$由Zoutendijk定理可知, 当在每个迭代点$x_k$处有$\cos \theta \ge \delta >0$, 一个线搜索算法就能全局收敛于稳定点, 其中线搜索需满足特定条件. 同样的, 我们选取${\mathcal{D}}_k$使得只要存在一个$p\in {\mathcal{D}}_k$使得$\cos \theta >\delta$, 这里不管$\nabla f_k$是如何:$$\kappa ({\mathcal{D}}_k)\stackrel{def}{=}\underset{v\in {\mathbb{R}}^n}{\min }\underset{p\in {\mathcal{D}}_k}{\max }\frac{v^{T}p}{\Vert v\Vert \Vert p\Vert }\ge \delta .$$若上述两个条件均满足, 我们就有对$\forall k$, $$-\nabla f_k^{T}p\ge \kappa ({\mathcal{D}}_k)\Vert \nabla f_k\Vert \Vert p\Vert \ge \delta {\beta }_{\min }\Vert \nabla f_k\Vert ,\exists p\in {\mathcal{D}}_k.$$我们只需选取符合条件$p$, 从而就有全局收敛.

   满足上述两个条件的${\mathcal{D}}_k$的例子包括坐标方向集 { e 1 , e 2 , … , e n , − e 1 , − e 2 , … , − e n } , \{e_1,e_2,\ldots,e_n,-e_1,-e_2,\ldots,-e_n\}, {e1​,e2​,…,en​,−e1​,−e2​,…,−en​},和$n+1$个向量$$p_i=\frac{1}{2n}e-e_i,i=1,2,\dots ,n;p_{n+1}=\frac{1}{2n}e$$构成的方向集. $n=3$的情形可见下图.
   

2. 坐标下降算法类似于算法2的特殊情形(每步迭代选取 i ∈ { 1 , 2 , … , n } i\in\{1,2,\ldots,n\} i∈{1,2,…,n}, D k = { e i , − e i } \mathcal{D}_k=\{e_i,-e_i\} Dk​={ei​,−ei​}). 注意这样的${\mathcal{D}}_k$有$\kappa ({\mathcal{D}}_k)=0,\forall k$. 因此$\cos \theta$在每步迭代中可以任意接近0.

3. 通常, 满足条件的方向仅是${\mathcal{D}}_k$的一个子集. 其他的方向可以根据$f$的信息、尺度或者根据先前的迭代启发式地选取. 也可选取它们为核心方向(保证$\delta >0$的那些)的线性组合.

4. 算法2并不需要我们选取具有最小目标函数值的$x_k+{\gamma }_k{p}_k,p_k\in {\mathcal{D}}_k$. 事实上, 我们可以由此节省些计算量.

5. 充分下降函数$\rho (t)$的选取也很重要. 若$\rho (\cdot )$选为0, 则任一带来下降的点都会被接收为新迭代点. 正如我们在第三章中所讨论的那般, 这样的弱条件一般并不能得到强全局收敛性. 更加常用且适宜的选择是$\rho (t)=M{t}^{3/2}$, 其中$M$为某个正常数.

4. 共轭方向法

从第五章中我们了解到对于严格凸的二次函数$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x,$$它的极小点是可以通过沿着$n$个共轭方向分别做$n$次一维极小化得到的. 这些方向在第五章中定义为梯度的线性组合. 而在这一节中, 我们将展示如何仅利用函数值构建共轭方向, 从而设计出仅需函数值的极小化二次函数的算法. 自然地, 我们也将考虑这一方法对于一般非线性目标函数$f$的延展.

我们将利用平行子空间性质. 以$n=2$为例, 考虑两条平行直线$l_1(\alpha )=x_1+\alpha p,l_2(\alpha )=x_2+\alpha p$, 其中$x_1,x_2,p$为${\mathbb{R}}^2$中的给定向量, $\alpha$为定义直线的标量参数. 设$x_1^{\ast },x_2^{\ast }$分别为$f(x)$沿$l_1,l_2$的极小点, 我们证明$x_1^{\ast }-x_2^{\ast }$与$p$相互共轭. 从而若再沿着连接$x_1^{\ast },x_2^{\ast }$的直线一维极小化, 便能达到$f$的极小点. 可见下图.

这就启发我们设计如下针对二维二次函数$f$的算法:

1. 选取线性无关的方向, 如坐标方向$e_1,e_2$. 从任意初始点$x_0$出发, 沿着$e_2$极小化$f$得到点$x_1$.
2. 再从$x_1$出发, 分别依次沿着$e_1,e_2$极小化$f$得到$z$. 由平行子空间性质, $z-x_1$与$e_2$相互共轭, 这是因为$x_1,z$分别是平行于$e_2$的两条线上的极小点.
3. 最后再从$x_1$出发, 沿着$z-x_1$极小化即得到$f$的极小点.

下面我们叙述最一般情形的平行子空间极小化性质. 设$x_1,x_2$为${\mathbb{R}}^n$中两个不同的点, { p 1 , p 2 , … , p l } \{p_1,p_2,\ldots,p_l\} {p1​,p2​,…,pl​}为${\mathbb{R}}^n$中线性无关的一组向量. 定义两个平行体:$$\begin{gathered} S_1=\left\{x_1+\sum _{i=1}^l{\alpha }_i{p}_i|{\alpha }_i\in \mathbb{R},i=1,2,\dots ,l\right\}, \\ {S}_2=\left\{x_2+\sum _{i=1}^l{\alpha }_i{p}_i|{\alpha }_i\in \mathbb{R},i=1,2,\dots ,l\right\}. \end{gathered}$$记$f$在$S_1,S_2$上的极小点分别为$x_1^{\ast },x_2^{\ast }$, 则我们有$x_2^{\ast }-x_1^{\ast }$与$p_1,p_2,\dots ,p_l$相互共轭. 事实上, 由极小化性质我们有$$\begin{gathered} {\frac{\partial f(x_1^{\ast }+{\alpha }_i{p}_i)}{\partial {\alpha }_i}|}_{{\alpha }_i=0}=\nabla f(x_1^{\ast }{)}^T{p}_i=0,i=1,2,\dots ,l, \\ {\frac{\partial f(x_2^{\ast }+{\alpha }_i{p}_i)}{\partial {\alpha }_i}|}_{{\alpha }_i=0}=\nabla f(x_2^{\ast }{)}^T{p}_i=0,i=1,2,\dots ,l. \end{gathered}$$因此$$\begin{array}{rl}0& =(\nabla f(x_1^{\ast })-\nabla f(x_2^{\ast }){)}^T{p}_i\\ & =(A{x}_1^{\ast }-b-A{x}_2^{\ast }+b{)}^T{p}_i\\ & =(x_1^{\ast }-x_2^{\ast }{)}^{T}A{p}_i,i=1,2,\dots ,l.\end{array}$$我们以$n=3$为例描述如何利用平行子空间性质产生三个共轭方向.

1. 首先选取线性无关方向, 比如$e_1,e_2,e_3$. 从任一点$x_0$出发, 首先沿着$e_3$得到点$x_1$.
2. 之后从$x_1$出发, 依次沿着$e_1,e_2,e_3$得到$z$.
3. 接着从$x_1$出发, 沿着$p_1=z-x_1$得到$x_2$. 如前述, $p_1=z-x_1$与$e_3$共轭. 注意$x_2$事实上为$f$在集合 S 1 = { y + α 1 e 3 + α 2 p 1 ∣ α 1 ∈ R , α 2 ∈ R } S_1=\{y+\alpha_1e_3+\alpha_2p_1|\alpha_1\in\mathbb{R},\alpha_2\in\mathbb{R}\} S1​={y+α1​e3​+α2​p1​∣α1​∈R,α2​∈R}上的极小点, 其中$y$是沿着$e_1,e_2$极小化后得到的中间点.

展开新一轮的迭代. 我们舍去$e_1$, 并定义新的搜索方向集为 { e 2 , e 3 , p 1 } \{e_2,e_3,p_1\} {e2​,e3​,p1​}.
4. 从$x_2$出发, 依次沿着$e_2,e_3,p_1$极小化$f$得到$\hat{z}$. 注意$\hat{z}$可看作是$f$在集合 S 2 = { y ^ + α 1 e 3 + α 2 p 1 ∣ α 1 ∈ R , α 2 ∈ R } S_2=\{\hat{y}+\alpha_1e_3+\alpha_2p_1|\alpha_1\in\mathbb{R},\alpha_2\in\mathbb{R}\} S2​={y^​+α1​e3​+α2​p1​∣α1​∈R,α2​∈R}, 其中$\hat{y}$为中间点. 因此应用平行子空间极小化性质于$S_1,S_2$, 我们就有$p_2=\hat{z}-x_2$与$e_3,p_1$共轭.
5. 接着沿着$p_2$得到$x_3$, 而这就是$f$的极小点.

这样我们就得到了共轭方向$e_3,p_1,p_2$.

下面我们就可以叙述适用于一般情形的算法了, 它由内迭代与外迭代构成. 内迭代中, 我们将沿着线性无关的方向集做$n$次一维极小化. 内迭代完成后, 即产生新的共轭方向, 再以之替换先前储存的搜索方向.

算法2 (DFO Method of Conjugate Directions)
Choose an initial point $x_0$ and set $p_i=e_i$, for $i=1,2,\dots ,n$;
Compute $x_1$ as the minimizer of $f$ along the line $x_0+\alpha p_n$;
Set $k\leftarrow 1$

repeat until a convergence test is satisfied
\quad\quad Set $z_1\leftarrow x_k$;
\quad\quad for $j=1,2,\dots ,n$
\quad\quad\quad\quad Calculate ${\alpha }_j$ so that $f(z_j+{\alpha }_j{p}_j)$ is minimized;
\quad\quad\quad\quad Set $z_{j+1}\leftarrow z_j+{\alpha }_j{p}_j$;
\quad\quad end (for)
\quad\quad Set $p_j\leftarrow p_{j+1}$ for $j=1,2,\dots ,n-1$ and $p_n\leftarrow p_n\leftarrow z_{n+1}-z_1$;
\quad\quad Calculate ${\alpha }_n$ so that $f(z_{n+1}+{\alpha }_n{p}_n)$ is minimized;
\quad\quad Set $x_{k+1}\leftarrow z_{n+1}+{\alpha }_n{p}_n$;
\quad\quad Set $k\leftarrow k+1$;
end (repeat)

算法说明:

1. 线搜索可通过使用搜索方向上的三个函数值做二次插值实施. 由于严格凸二次函数限制在一条线上仍然是严格凸的二次函数, 因此在不计误差的情况下插值带来的是精确的模型, 从而能够计算精确的一维极小点. 注意在外迭代第$k$步结束后, $p_{n-k},p_{n-k+1},\dots ,p_n$相互共轭. 因此若共轭方向均非零, 算法在$n-1$步外迭代后终止于$f$的极小点. 不幸的是, 这种可能性无法排除. 因此我们必须采取一些防护措施以增强鲁棒性. 一般地, 算法3在$n-1$步外迭代后终止, 总共需要计算$O(n^2)$量阶的函数值.
2. 算法3可推广至极小非二次的目标函数的情形. 这里唯一的变动在于线搜索上: 我们只能利用插值近似求解. 由于可能非凸的存在, 因此我们需要特别关注一维搜索的过程.
3. 数值实验表明算法3的推广对于维数较小的问题表现良好, 不过有时 { p i } \{p_i\} {pi​}变得越来越线性相关. 一种改进的方法是, 先定义 { p i } \{p_i\} {pi​}共轭程度的一种度量. 定义缩放的方向$${\hat{p}}_i=\frac{p_i}{\sqrt{p_i^{T}A{p}_i}},i=1,2,\dots ,n.$$我们可以证明$$|\det ({\hat{p}}_1,{\hat{p}}_2,\dots ,{\hat{p}}_n)|$$达到极大当且仅当$p_i$相互共轭. 这一结论表明, 若以最近产生的共轭方向替代 { p 1 , p 2 , … , p n } \{p_1,p_2,\ldots,p_n\} {p1​,p2​,…,pn​}中的某一元会导致行列式下降, 就应当取消替代的操作.
   程序4针对二次目标函数的情形实施了上述策略. 代数运算使得我们无需Hessian矩阵$A$便能计算缩放的方向${\hat{p}}_i$, 这是因为$p_i^{T}A{p}_i$在沿着$p_i$做线搜索时就有了. 进一步地, 我们在判定行列式不增时仅需比较计算过的函数值. 下述程序将在算法3中内迭代完成后被立即调用.

程序4 (Updating of the Set of Directions)
Find the integer m ∈ { 1 , 2 , … , n } m\in\{1,2,\ldots,n\} m∈{1,2,…,n} such that ${\psi }_m=f(x_{m-1})-f(x_m)$ is maximized;
Let $f_1=f(z_1),f_2=f(z_{n+1}),f_3=f(2z_{n+1}-z_1)$;
if $f_3\ge f_1$ or $(f_1-2f_2+f_3)(f_1-f_2-{\psi }_m{)}^2\ge \frac{1}{2}{\psi }_m(f_1-f_3{)}^2$
\quad\quad Keep the set $p_1,p_2,\dots ,p_n$ unchanged and set $x_{k+1}\leftarrow z_{n+1}$;
else
\quad\quad Set $\hat{p}\leftarrow z_{n+1}-z_1$ and calculate $\hat{\alpha }$ so that $f(z_{n+1}+\hat{\alpha }\hat{p})$ is minimized;
\quad\quad Set $x_{k+1}\leftarrow z_{n+1}+\alpha \hat{p}$;
\quad\quad Remove $p_m$ from the set of directions and add $\hat{p}$ to this set;
end (if)

上述程序同样可应用于一般目标函数上非精确一维线搜索的情形. 最终得到的共轭梯度法对于求解小规模问题是很有用的.

5. Nelder-Mead方法

Nelder-Mead单纯形反射方法自其1965年被发明以来一直是广受欢迎的DFO算法. 其名称来源于在算法的任意阶段, 我们都保持对${\mathbb{R}}^n$中$n+1$个点的追索, 它们形成的凸腔构成单纯形. (这一方法与线性规划中的单纯形法没有任何关联.) 给定有顶点 { z 1 , z 2 , … , z n + 1 } \{z_1,z_2,\ldots,z_{n+1}\} {z1​,z2​,…,zn+1​}的单纯形$S$, 我们可定义关联矩阵$V(S)$为$$V(S)=[z_2-z_1,z_3-z_1,\dots ,z_{n+1}-z_1].$$若$V$非奇异, 则称单纯形非退化或非奇异(例如对于${\mathbb{R}}^3$中的单纯形, 若其四个顶点不共面, 则就非奇异).

Nelder-Mead算法中的单次迭代中, 我们移除具有最差函数值的顶点并代之以具更好值的点. 而新点则是通过沿着连接最差顶点与剩余顶点中心的线反射、扩张或收缩单纯形得到. 若以这种方式我们无法找到更好的点, 我们就仅保持具最好函数值的点不动, 而令其他的顶点向该值移动以收缩单纯形.

下面我们来详细描述算法的一次迭代过程. 先定义一些符号: 当前单纯形的$n+1$个顶点记作 { x 1 , x 2 , … , x n + 1 } \{x_1,x_2,\ldots,x_{n+1}\} {x1​,x2​,…,xn+1​}, 其中排列的顺序满足$$f(x_1)\le f(x_2)\le \cdots \le f(x_{n+1}).$$最佳的$n$个点的中心定义为$$\bar{x}=\sum _{i=1}^n{x}_i.$$连接$\bar{x}$与最差点$x_{n+1}$的线有参数化$$\bar{x}(t)=\bar{x}+t(x_{n+1}-\bar{x}).$$

程序5 (One Step of Nelder-Mead Simplex)
Compute the reflection point $\bar{x}(-1)$ and evaluate $f_{-1}=f(\bar{x}(-1))$;
if $f(x_1)\le f_{-1}<f(x_n)$
\quad\quad (reflect point is neither best nor worst in the new simplex)
\quad\quad replace $x_{n+1}$ by $\overline{-1}$ and go to next iteration;
else if $f_{-1}$<f(x_1)$
\quad\quad (reflect point is better than the current best; try to go farther along this direction)
\quad\quad Compute the expansion point $\bar{x}(-2)$ and evaluate $f_{-2}=f(\bar{x}(-2))$;
\quad\quad if $f_{-2}<f_{-1}$
\quad\quad\quad\quad replace $x_{n+1}$ by $x_{-2}$ and go to next iteration;
\quad\quad else
\quad\quad\quad\quad replace $x_{n+1}$ by $x_{-1}$ and go to next iteration;
else if $f_{-1}\ge f(x_n)$
\quad\quad (reflect point is still worse than $x_n$; contract)
\quad\quad if $f(x_n)\le f_{-1}<f(x_{n+1})$
\quad\quad\quad\quad (try to perform “outside” contraction)
\quad\quad\quad\quad evaluate $f_{-1/2}=\bar{x}(-1/2)$;
\quad\quad\quad\quad if $f_{-1/2}\le f_{-1}$
\quad\quad\quad\quad\quad\quad replace $x_{n+1}$ by $x_{-1/2}$ and go to next iteration;
\quad\quad else
\quad\quad\quad\quad (try to perform “inside” contraction)
\quad\quad\quad\quad evaluate $f_{1/2}=\bar{x}(1/2)$;
\quad\quad\quad\quad if $f_{1/2}<f_{n+1}$
\quad\quad\quad\quad\quad\quad replace $x_{n+1}$ by $x_{1/2}$ and go to next iteration;
\quad\quad (neither outside nor inside contraction was acceptable; shrink the simplex toward $x_1$)
\quad\quad replace $x_i\leftarrow (1/2)(x_1+x_i)$ for $i=2,3,\dots ,n+1$;

程序5在$n=3$下的一个例子如下图所示.

其中最差点为$x_3$, 可能的替换点为$\bar{x}(-1),\bar{x}(-2),\bar{x}(-\frac{1}{2}),\bar{x}(\frac{1}{2})$. 若无一可用, 则收缩单纯形为图中虚线框定的小三角形, 其中保持$x_1$不动. 标量$t$在我们的讨论中只选定了一些特殊值$-1,-2,-\frac{1}{2},\frac{1}{2}$. 但这并不是说其他的选择就不好, 不过这需要另外加上特定的限制.

Nelder-Mead算法的实际表现通常是较为合理的, 尽管在非最优点处可能会有停滞. 此时重启是个较好的选择. 注意除非以收缩单纯形结束, 算法的平均函数值$$\frac{1}{n+1}\sum _{i=1}^{n+1}f(x_i)$$在每一步均会下降. 进一步, 若$f$凸, 收缩也不会增加平均函数值.
对于Nelder-Mead方法的收敛理论现仍较有限.

6. 隐式滤波法

最后我们介绍函数值计算基于模型$f(x)=h(x)+\varphi (x)$的算法——隐式滤波法——其中$f$是光滑的. 它的最简形式即为最速下降法的变体, 其中梯度$\nabla f_k$以有限差分估计代替, 差分参数$ϵ$也不必过小.

当函数随着迭代接近解时噪声水平不断下降, 隐式滤波法将发挥最佳的效果. 这种情形在我们控制噪声水平时是可能发生的, 例如当$f$是通过在人为给定的容忍限内求解微分方程得到的, 亦或者通过经人为指定数量的随机模拟后得到的(模拟的次数越多方差越小, 噪声水平越低). 隐式滤波算法在给定当前$x$处的噪声水平的情况下, 系统地(但是, 人们希望不会像误差衰减那么快)减小$ϵ$, 从而保证${\nabla }_{ϵ}f(x)$具有合理精度. 对每个$ϵ$的值, 它都将经过一次内循环——利用搜索方向$-{\nabla }_{ϵ}f(x)$的Armijo线搜索. 若内循环无法在至少$a_{\max }$步中得到较为满意的步长, 我们就回到外循环, 选取更小的$ϵ$重复以上过程. 算法的证实表述如下.

算法6 (Implicit Filtering)
Choose a sequence { ϵ k } ↓ 0 \{\epsilon_k\}\downarrow 0 {ϵk​}↓0, Armijo parameters $c$ and $\rho$ in $(0,1)$, maximum bactracking parameter $a_{\max }$;
Set $k\leftarrow 1$, Choose initial point $x=x_0$;
repeat
\quad\quad increment_k$\leftarrow$false;
\quad\quad repeat
\quad\quad\quad\quad Compute $f(x)$ and ${\nabla }_{{ϵ}_k}f(x)$;
\quad\quad\quad\quad if $\Vert {\nabla }_{{ϵ}_k}f(x)\Vert \le {ϵ}_k$
\quad\quad\quad\quad\quad\quad increment_k$\leftarrow$true;
\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad Find the smallest integet $m$ between 0 and $a_{\max }$ s. t.
\quad\quad\quad\quad\quad\quad $f(x-{\rho }^m{\nabla }_{{ϵ}_k}f(x))\le f(x)-c{\rho }^m\Vert {\nabla }_{{ϵ}_k}f(x){\Vert }_2^2$;
\quad\quad\quad\quad\quad\quad if no such $m$ exists
\quad\quad\quad\quad\quad\quad\quad\quad increment_k$\leftarrow$true;
\quad\quad\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad\quad\quad $x\leftarrow x-{\rho }^m{\nabla }_{{ϵ}_k}f(x)$;
\quad\quad until increment_k;
\quad\quad $x_k\leftarrow x;k\leftarrow k+1$;
until a termination test is satisfied.

注意到算法6的内循环本质上就是第三章中的回溯线搜索算法, 附加上了一个收敛的判定准则: 判定是否在差分参数${ϵ}_k$所对应的精度下求得了极小. 若梯度${\nabla }_{{ϵ}_k}f$很小, 亦或是线搜索无法找到令人满意的迭代点(说明当前的${\nabla }_{{ϵ}_k}f(x)$不足以产生充分的下降), 我们就减小差分参数至${ϵ}_{k+1}$继续内循环.
算法6的一个基本的收敛结果如下.

定理2设${\nabla }^{2}h$Lipschitz连续, 算法6产生无限的迭代序列 { x k } \{x_k\} {xk​}, 且$$\underset{k\to \infty }{\lim }{ϵ}_k^2+\frac{\eta (x_k;{ϵ}_k)}{{ϵ}_k}=0.$$另假设仅有有限个内循环不以$\Vert {\nabla }_{{ϵ}_k}f(x_k)\Vert \le {ϵ}_k$结束. 则序列 { x k } \{x_k\} {xk​}的所有聚点均是稳定点.
证明: 由于 { ϵ k } ↓ 0 \{\epsilon_k\}\downarrow0 {ϵk​}↓0, 所以基于假设条件我们有${\nabla }_{{ϵ}_k}f(x_k)\to 0$. 因此$\nabla h(x_k)\to 0$. 从而所有聚点满足$\nabla h(x)=0$. 证毕.

更加复杂的隐式滤波法则使用梯度估计${\nabla }_{{ϵ}_k}f$构建拟牛顿近似Hessian, 从而产生拟牛顿搜索方向而不是算法6中的近似负梯度方向.