第四章 共轭梯度法和拟牛顿法

4 共轭梯度法和拟牛顿法

本章讨论地球物理反演中常用的两种算法，共轭梯度法和拟牛顿法。与牛顿法不同的是，这两种方法不需要显式计算Hessian的逆，而是沿着能够减小数据误差的方向迭代地移动下降方向。每次迭代的计算复杂度为$O(N^2)$次矩阵乘法。这些方法的另一个优势在于，二者均无需存储或显式计算Hessian。它们的缺点是无法保证快速收敛。但是，他们仍然比无预处理的最速下降法更快。

4.1 共轭梯度法

Hestenes和Stiefel（1952）提出的共轭梯度法（Conjugate-gradient，CG）是一种沿着相互共轭的方向进行线性最小化的有效方法。相比于牛顿方法，其每次迭代只需$O(N^2)$次代数运算，对于well-conditioned的正定矩阵$\rm \textbf H$，通常在远少于N次的迭代中即可实现收敛，其还具有编程简单，每次迭代只占用$O(N)$存储的优势。

共轭定义

当矢量$\rm \textbf {u, v}$与矩阵$\rm \textbf H$满足$(\rm \textbf u, \rm \textbf {Hv})=0$时称$\rm \textbf {u, v}$关于矩阵$\rm \textbf H$共轭。相互共轭的向量是线性独立的，并且跨越了一个$N$维的欧几里得空间。因此，任何$N\times 1$的解向量都可以扩展为N个独立共轭向量的线性组合。为了理解共轭方向的物理意义，我们首先定义拟牛顿公式。

拟牛顿公式

二阶导数的有限差分近似值与相邻点的一阶导数之差成正比：

$$\frac{d^2f(x)}{dx^2} \Delta x \approx \frac{df(x+\Delta x}{dx} -\frac{df(x)}{dx}, \tag{4.1}$$

其中，$\Delta x$是两点之间的空间增量。上述方程是拟牛顿公式的特例，其将Hessian矩阵与相邻点之间的梯度联系在了一起。拟牛顿（QN）公式定义为：

$${\rm \textbf H}({\rm \textbf x}^{k+1}-{\rm \textbf x}^{(k)}) ={\rm \textbf g}^{(k+1)}-{\rm \textbf g}^{(k)}, \tag{4.2}$$

上述公式可以通过式2.9推导得到：

$$\overset {\rm\textbf g^{(k+1)}} {\nabla f ({\rm\textbf x}_{0}+\Delta{\rm\textbf x})} =\overset {\rm\textbf g^{(k)}}{\rm\textbf g(\rm\textbf x_0)} + {\rm\textbf H\Delta {\rm\textbf x}} \tag{4.3}$$

定义${\rm \textbf x}_o \rightarrow {\rm \textbf x}^{(k)}$以及$\Delta {\rm \textbf x} \rightarrow {\rm \textbf x}^{(k+1)} - {\rm \textbf x}^{(k)}$公式4.3即变为4.2。
如果将${\rm \textbf x}^{(k+1)} - {\rm \textbf x}^{(k)}$归一化为单位矢量，那么沿此方向的函数曲率为$({\rm \textbf x}^{(k+1)} - {\rm \textbf x}^{(k)})^T {\rm \textbf H} ({\rm \textbf x}^{(k+1)} - {\rm \textbf x}^{(k)})$。QN公式将用于推导共轭梯度法。

共轭矢量的几何解释

对于SPD的Hessian，可将共轭矢量解释为指向图4.1b中${\rm \textbf d}^{(k-1)}$和${\rm \textbf g}^{(k)}$这两个矢量所跨越的平面上的损失靶心。
为了证明这一点，我们首先应注意到之前的搜索方向${\rm \textbf d}^{(k-1)}$和它的相关梯度${\rm \textbf g}^{(k)}$是相互垂直的，因此：

$$({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k)})=0 \tag{4.4}$$

矢量${\rm \textbf d}^{(k-1)}$和${\rm \textbf g}^{(k)}$跨越了如图4.1b所示的平面，损失等值线在这个平面上形成了一个局部靶心，这个局部靶心是全局最小值${\rm \textbf x}^{\ast}$的垂直投影。如果定义第$k+1$次的梯度${\rm \textbf g}^{(k+1)}$在局部靶心上，这时梯度将不存在沿着${\rm \textbf d}^{(k-1)}$和${\rm \textbf g}^{(k)}$所跨越平面的分量，所以${\rm \textbf g}^{(k+1)}$垂直于平面中的任何矢量，此时有：

$$({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k+1)})=0 \tag{4.5}$$

将${\rm \textbf d}^{(k-1)}$的内积带入拟牛顿公式4.2可得

$$({\rm \textbf d}^{(k-1)}, {\rm \textbf {Hd}}^{(k)}) =\overset {e.q. 4.5} {({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k+1)})} -\overset {e.q. 4.4} {({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k)})} =0, \tag{4.6}$$

其中${\rm \textbf x}^{(k+1)}-{\rm \textbf x}^{(k)} \rightarrow {\rm \textbf d}^{(k)}$。这个方程说明，一个终点在局部靶心处的搜索方向${\rm \textbf d}^{(k)}$，它与之前的搜索方向${\rm \textbf d}^{(k-1)}$是共轭的，其中靶心在${\rm \textbf d}^{(k-1)}$和${\rm \textbf g}^{(k)}$所跨越的平面内，如图4.1b所示。

4.1.1 共轭梯度算法

我们现在可以使用正定矩阵$\rm \textbf H$来构造一系列相互共轭的搜索向量以最小化二次损失函数。以下即为共轭梯度法，其使用最速下降方向作为初始搜索方向来寻找最优解。

1. 定义$k=1$时的初始模型${\rm \textbf x}^{(k-1)} \rightarrow {\rm \textbf x}^{(0)}$（如图4.1b），同时定义初始下降矢量${\rm \textbf d}^{(k-1)} \rightarrow {\rm \textbf d}^{(0)} = -{\rm \textbf g}^{(0)}$，迭代公式为：

   $${\mathbf x}^{(k)} ={\mathbf x}^{(k-1)}+ \alpha {\mathbf d}^{(k-1)}. \tag{4.7}$$

   其中，$\alpha$通过线最小化方法来使$f({\rm \textbf x}^{(k-1)} + \alpha {\rm \textbf d}^{(k-1)})$最小求得。根据公式3.10，$\alpha$可进一步表示为：

   $$\alpha=- \frac{({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k-1)})} {({\rm \textbf d}^{(k-1)}, {\rm \textbf {Hd}}^{(k-1)})} \tag{4.8}$$

2. 以${\rm \textbf x}^{(k)}$为初始模型，下次移动方向${\rm \textbf d}^{(k)}$需要与${\rm \textbf d}^{(k-1)}$共轭，并且其需要在${\rm \textbf g}^{(k)}=\nabla f({\rm \textbf x}^{(k)})$和${\rm \textbf d}^{(k-1)}$所跨越的平面内，即：

$${\rm \textbf d}^{(k)} ={\rm \textbf g}^{(k-1)} + \beta {\rm \textbf d}^{(k-1)}. \tag{4.9}$$

其中，$\beta$用于保证${\rm \textbf d}^{(k)}$与${\rm \textbf d}^{(k-1)}$共轭。将共轭定义带入公式4.9可得$\beta$为：

$$\beta= \frac{({\rm \textbf d}^{(k-1)}, {\rm \textbf {Hg}}^{(k)})} {({\rm \textbf d}^{(k-1)}, {\rm \textbf {Hd}}^{(k-1)})} \tag{4.10}$$

根据公式4.7可知，${\rm \textbf {Hd}}^{(k-1)}=({\rm \textbf g}^{(k)} - {\rm \textbf g}^{(k-1)})/\alpha$，结合公式4.2，公式4.10可以写为：

$$\begin{align} \notag \beta &= ({\rm \textbf g}^{(k)} - {\rm \textbf g}^{(k-1)}, {\rm \textbf g}^{(k)}) / (\alpha {\rm \textbf d}^{(k-1)}, {\rm \textbf {Hd}}^{(k-1)}) \\ \notag &=({\rm \textbf g}^{(k)} - {\rm \textbf g}^{(k-1)}, {\rm \textbf g}^{(k)}) / ({\rm \textbf d}^{(k-1)}, {\rm \textbf g}^{(k)} - {\rm \textbf g}^{(k-1)}) \\ \notag &=({\rm \textbf g}^{(k)} - {\rm \textbf g}^{(k-1)}, {\rm \textbf g}^{(k)}) / ({\rm \textbf g}^{(k-1)}, {\rm \textbf g}^{(k-1)}), \end{align} \tag{4.11}$$

在4.11最后一行公式中，我们需要将${\rm \textbf d}^{(k-1)} = {\rm \textbf g}^{(k-1)} + \beta^{‘} {\rm \textbf d}^{(k-2)}$带入倒数第二行公式中，如图4.2b所示，第k次梯度垂直于k-1次的搜索方向，所以化简后即可得到最后一行公式。

公式4.7-4.11可以总结为如下：

$$\alpha ^{(k)}=- \frac{({\rm \textbf d}^{(k)}, {\rm \textbf g}^{(k)})} {({\rm \textbf d}^{(k)}, {\rm \textbf {Hd}}^{(k)})} \tag{4.12}$$ $${\rm \textbf x}^{(k+1)} ={\rm \textbf x}^{(k)} + \alpha ^{(k)} {\rm \textbf d}^{(k)}. \tag{4.13}$$ $${\rm \textbf g}^{(k+1)}=\nabla f({\rm \textbf x}^{(k+1)}) \tag{4.14}$$ $${\rm \textbf d}^{(k+1)} =-{\rm \textbf g}^{(k+1)} + \beta {\rm \textbf d}^{(k)}. \tag{4.15}$$ $$\beta = ({\rm \textbf g}^{(k+1)} - {\rm \textbf g}^{(k)}, {\rm \textbf g}^{(k+1)}) / ({\rm \textbf g}^{(k)}, {\rm \textbf g}^{(k)}), \tag{4.16}$$

CG法的几何解释如图4.1b所示，共轭方向${\rm \textbf d}^{(k)}$在${\rm \textbf d}^{(k-1)}$和${\rm \textbf g}^{(k)}$所跨越的平面内，且这个共轭矢量指向靶心。这里的靶心并非是全局最小，而是局部平面内的局部极小。

公式4.13-4.16描述的是Polak-Ribiere CG法（Gill et al., 1981）。非常明显的是，由于${\rm \textbf g}^{(k)}$和${\mathbf g}^{(k+1)}$是正交的(Nocedal and Wright，1999)，公式4.16中的$\beta$也可以写为：

$$\beta = ({\rm \textbf g}^{(k+1)}, {\rm \textbf g}^{(k+1)}) / ({\rm \textbf g}^{(k)}, {\rm \textbf g}^{(k)}), \tag{4.17}$$

将4.17代替4.16即可得到Fletcher-Reeves CG法（Gill et al., 1981）。
另一种求取$\alpha$的方法是Hestenes & Stiefel(1952)法：

$$\alpha ^{(k)}=- \frac{({\rm \textbf g}^{(k+1)}, {\rm \textbf g}^{(k+1)})} {({\rm \textbf g}^{(k+1)}, {\rm \textbf {Hg}}^{(k)})} \tag{4.18}$$

对于严格的二次函数以及无限精度的算法，Polak-Ribiere和Fletcher-Reeves CG法是等价的。但是，对于非二次函数，${\rm \textbf d}^{(k)}$也可能与${\rm \textbf g}^{(k)}$正交（由于计算机的舍入误差也会导致丧失正交性）。这种情况下，${\rm \textbf g}^{(k+1)} \approx {\rm \textbf g}^{(k)}$，所以公式4.16中所示的Polak-Ribiere方法会使得$\beta =0$且${\rm \textbf d}^{(k)}=-{\rm \textbf g}^{(k)}$，此时Polak-Ribiere CG法将会重置为最速下降法，这也使得Polak-Ribiere方法更稳定。
所以，当算法迭代速度较慢时，${\rm \textbf g}^{(k+1)} \approx {\rm \textbf g}^{(k)}$以致于$\beta =0$并且公式4.15变为${\rm \textbf d}^{(k+1)}=-{\rm \textbf g}^{(k+1)}$，算法将自动重置为最速下降法。Polak-Ribiere CG要比Fletcher-Reeves更好一些（Nocedal & Wright, 1999）。

4.1.2 共轭梯度算法求解线性方程系统

共轭梯度法同样可以求解$N \times N$线性方程组。将正定线性方程组表示为：

$$\mathbf {Lx} = \mathbf t, \tag{4.19}$$

其中，$\mathbf t$和$\mathbf x$均为$N \times 1$矢量，$\mathbf L$为$N \times N$正定矩阵。求解以上方程组与最小化以下损失函数问题等价：

$$f(\mathbf x) = \frac{1}{2} \mathbf x^T \mathbf {Lx} - \mathbf x^T \mathbf t, \tag{4.20}$$

上述损失函数的梯度可表示为：

$$\nabla f(\mathbf x) = \mathbf {Lx} - \mathbf t. \tag{4.2}$$

能够最小化方程4.20的矢量$\mathbf x^{\ast}$同样也满足$\nabla f(\mathbf x) = \mathbf {Lx} - \mathbf t=0$，同时其也是方程组4.19的解。因此，可以沿着相互共轭的方向来迭代地最小化损失函数进而求解$\mathbf {Lx}=\mathbf t$。

4.1.3 Matlab共轭梯度代码

% CG Iterations to solve SPD eqn Lx=t
% where gradient gk = L*xk-t;
% and iterate soln xk1 = xk + alpha*dk1
%
clear all;
L = [400 1; 1 4]; % Simple 2x2
t = [5 2]’; % problem
[m,n] = size(L);
nit = 3;
%
xk = zeros(n,1); % Starting point at origin
gk = (L*xk-t);dk=-gk;
%
for q = 1:nit;
    alpha = -dk’*gk/(dk’*L*dk);
    xk1 = xk + alpha*dk; % Compute new iterate
    gk1 = gk+alpha*L*dk; % Find new gradient
    beta = gk1’*gk1/(gk’*gk); % Fletcher-Reeves
    beta = gk1’*(gk1-gk)/(gk’*gk); % Polak-Ribiere
    dk1 = -gk1+beta*dk; % Find new conj. direction
    dk = dk1; % Refresh iterates
    gk = gk1;
    xk = xk1
end
%

如果方程组是超定的，问题转化为求解$\mathbf x$使其满足一般方程$\mathbf L^T \mathbf {Lx}=\mathbf L ^T \mathbf t$。与基于最速下降法求解超定方程组相类似，需要迭代地最小化损失函数（误差的平方和）来求解$\mathbf x$，损失函数可表示为：

$$\begin{align} \notag f(\mathbf x) &= \frac{1}{2}(\mathbf {Lx-t})^T(\mathbf {Lx-t}) \\ \notag &=\frac{1}{2} \mathbf x^T \mathbf L^T \mathbf {Lx} - \mathbf t^T \mathbf {Lx} + \frac{1}{2} \mathbf t^T \mathbf t. \end{align} \tag{4.22}$$

$f(\mathbf x)$是关于$\mathbf x$的二次函数，等号右侧公式可视为损失函数在原点关于$\mathbf x$的展开，梯度可表示为：

$$\nabla f(\mathbf x) = \mathbf L^T (\mathbf {Lx} - \mathbf t). \tag{4.23}$$

可通过沿着相互共轭的方向迭代最小化来求解。

计算梯度中的矩阵乘法时，应首先计算残差$\mathbf r$（即$\mathbf {r=Lx-t}$）,然后再将其乘以$\mathbf L^T$。相比于矩阵-矩阵乘法$\mathbf L^T \mathbf L$，由于只计算了矩阵-矢量的乘法，所以可以减少计算量并且减少舍入误差。

以下为CG法求解超定线性方程组的Matlab代码：

% CG Iterations to solve overdetermined eqns Lx=b
% where gradient gk = L’*(L*xk-b);
% and iterate soln xk1 = xk + alpha*dk1
%
clear all;
L = [400 1; 1 4; 0 1]; % Simple 3x2
b = [5 2 .5]’; % problem
[m,n] = size(L);
nit = 3;
%
xk = zeros(n,1); % Starting point at origin
gk = L’*(L*xk-b);dk=-gk;
%
for q = 1:nit;
    res = L*dk;
    res1 = L’*res;
    alpha = -dk’*gk/(res’*res);
    xk1 = xk + alpha*dk; % Compute new iterate
    gk1 = gk+alpha*res1; % Find new gradient
    beta = gk1’*gk1/(gk’*gk); % Fletcher-Reeves
    beta = gk1’*(gk1-gk)/(gk’*gk); % Polak-Ribiere
    dk1 = -gk1+beta*dk; % Find new conj. direction
    dk = dk1; % Refresh iterates
    gk = gk1;
    xk = xk1
end

如果$N\times N$矩阵$\mathbf L$的特征值是$\lambda_i$，那么$\mathbf L^T \mathbf L$特征值为$\lambda_i^2$。这也意味着$\mathbf L^T \mathbf L$的条件数将会是$\mathbf L$的平方，这将会降低收敛效率，所以需要预处理。当方程组是ill-conditioned时，也可以使用LSQR算法来求解（Paige and Saunders，1982）。

Figure 3.4显示了使用Polak Ribiere共轭梯度法求解Rosenbrock函数的例子。这个例子中，CG方法显著优于最速下降法，但是不如牛顿法。对于非线性泛函，应该在第N次迭代时或之前将CG方向重置为梯度方向，其中$N$是向量$\mathbf x$的维数；每隔几次迭代就需要将CG方向重置为最速下降方向，因为搜索方向总是会失去相互共轭性。

4.1.4 收敛速度

对于线性问题，如果$N \times N$Hessian是正定的，那么CG算法至多需要$N$次就会收敛。如果$\mathbf H$的特征值在一个值附近，那么就可以加速收敛。在Nocedal & Wright（1999）的定理5.4中，如果$\mathbf H$有$r$个特征值，那么CG算法至多需要$r$次迭代就会停止。例如，如果一个$5\times 5$矩阵$\mathbf H$的特征值为$(5, 5, 2, 2, 2)$，那么只需要两次迭代即可。如果特征值的分布是$(5, 5.1, 4.8, 2, 2)$，则可能需要3次迭代才可收敛到理想的结果，其中两次迭代以特征值5为中心，一次迭代以2为中心。对此的直观解释是，与聚集特征值相关的损失函数等值面形成了一个超球体，其中一个搜索方向指向该球体的中心。

Nocedal and Wright (1999)定理5.5中，CG第$k+1$次迭代的结果$\mathbf x^{(k+1)}$与真实解$\mathbf x^{\ast}$的差可表示为：

$$||\mathbf x^{(k+1)} - \mathbf x^{\ast}||\approx \epsilon ||\mathbf x^{(0)}- \mathbf x^{\ast}|| \tag{4.24}$$

其中$\epsilon = \sigma_{N-k}- \sigma_1$是$\mathbf H$的第$N-k$个特征值和第一个特征值（最小）（特征值按升序排列）。对于一个较小的$\epsilon$，这个公式说明基于CG法，数次迭代就可以提供一个较好的估计。图4.2示意性地说明了四个不同特征值的收敛速度，其余特征值围绕一个点聚集。

共轭梯度法的收敛速度随着Hessian矩阵条件数的减少而增加。事实上，如果Hessian是SPD并且只有r个不连续的特征值，CG迭代将终止于最多r次迭代的解（Nocedal & Wright，1999，定理5.4）。如果Hessian的特征值是相等的，则损失函数的等值线为环形，因此CG方法的第一次迭代沿着指向靶心的最速下降路径。对于高度椭圆的等值线，收敛可能很慢，因此应使用预处理或正则化。

4.2 拟牛顿法

牛顿法的主要问题在于$f(\mathbf x)$的二阶导数或者Hessian的逆非常难求，为了避免这些问题，可以使用拟牛顿法，其迭代公式为：

$$\mathbf x^{(k+1)} = \mathbf x^{(k)} - \mathbf H^{-1}_k \mathbf g^{(k)}, \tag{4.31}$$

其中，$\mathbf H^{-1}_k$是$\mathbf H^{-1}$的近似，其通过$f(\mathbf x)$的一阶导（而非二阶导）迭代更新求取。对于$N\times N$SPD的Hessian(Scales, 1985)：

$$\underset {k \rightarrow N} {\rm {lim}} \mathbf H^{-1}_k = \mathbf H^{-1}. \tag{4.32}$$

与牛顿方法相比，准牛顿（QN）方法的优点在于，每次更新只需要$O(N^2)$次代数运算，并且只计算一阶导数；与CG方法不同的是，其不需要周期性地重新开始迭代。然而，与CG方法的$O(N)$存储相比，仍然需要$O(N^2)$存储。降低存储的一种方法是使用limited-memory QN方法（Nocedal和Wright，1999），其只存储了几个向量来隐式构建$\mathbf H$的逆。

Rank-one 拟牛顿法

rank-one法中，设置4.31中初始近似$\mathbf H^{-1}$为$\mathbf H^{-1}_0 = \mathbf I$，此时有：

$$\mathbf x^{(1)} = \mathbf x^{(0)} -\alpha \mathbf H^{-1}_0 \mathbf g^{(0)} \tag{4.33}$$

其中，$\alpha$是通过先搜索方法最小化$f(\mathbf x)$得到的。
然后使更新后的$\mathbf H^{-1}_1$满足4.2准牛顿公式来更新$\mathbf H^{-1}_0$：

$$\mathbf x^{(1)} -\mathbf x^{(0)} = \mathbf H^{-1}_1 (\mathbf g^{(1)} -\mathbf g^{(0)}) \tag{4.34}$$

要注意的是，$\mathbf H^{-1}_0$不满足准牛顿公式，因为其只是Hessian逆的近似，为了保证4.34成立，通过rank-one更新构建$\mathbf H^{-1}_1$：

$$\mathbf H^{-1}_1 = \mathbf H^{-1}_0 + a \mathbf {uu}^T \tag{4.35}$$

其中$\mathbf u$是一个$N \times 1$矢量，$a$是常数。将4.35带入4.34即可得到$\mathbf u$和$a$：

$$\Delta \mathbf x^{(0)} = (\mathbf H^{-1}_0 + a \mathbf {uu}^T) \Delta \mathbf g^{(0)} \tag{4.36}$$

整理上式可得：

$$a \mathbf {uu}^T \Delta \mathbf g^{(0)} = \Delta \mathbf x^{(0)} - \mathbf H^{-1}_0 \Delta \mathbf g^{(0)} \tag{4.37}$$

其中，$\Delta \mathbf x^{(0)} = \mathbf x^{(1)} -\mathbf x^{(0)}$，$\Delta \mathbf g^{(0)} = \mathbf g^{(1)} -\mathbf g^{(0)}$。$\mathbf {uu}^T$形式的矩阵被称为仅具有一个线性独立行向量的秩为1的矩阵。$\mathbf {uu}^T$的每一列都是向量$\mathbf u$的倍数，每一行都是向量$\mathbf u^T$的倍数。$\mathbf I + a \mathbf {uu}^T$形式的矩阵被称为初等矩阵。
上述公式中有$N+1$个未知数，但只有$N$个方程。因此一个可能的解为（Gill et al., 1981）：

$$\mathbf u = \Delta \mathbf x ^{(0)} - \mathbf H^{-1}_0 \Delta \mathbf g ^{(0)} \tag{4.38}$$

以及

$$a = \frac{1}{\mathbf u ^T \Delta \mathbf g ^{(0)}} \tag{4.39}$$

将4.38和4.39代入公式4.35可得$\mathbf H^{-1}_1$的迭代解：

$$\mathbf H^{-1}_{k+1} = \mathbf H^{-1}_{k} + \frac {(\Delta \mathbf x ^{(k)} - \mathbf H^{-1}_{k} \Delta \mathbf g ^{(k)}) (\Delta \mathbf x ^{(k)} - \mathbf H^{-1}_{k} \Delta \mathbf g ^{(k)}) ^T} {(\Delta \mathbf x ^{(k)} - \mathbf H^{-1}_{k} \Delta \mathbf g ^{(k)})^T \Delta \mathbf g ^{(k)}} \tag{4.40}$$

以及

$$\mathbf x ^{(k+1)} = \mathbf x ^{(k)} - \alpha \mathbf H ^{-1}_k \mathbf g ^{(k+1)} \tag{4.41}$$

rank-one拟牛顿法存在两个问题。第一，$\mathbf g ^{(k)}$通常并不能保证为正定矩阵，因此优化过程并不一定是下降的。地儿，4.40中的分母可能接近或者等于0，会导致不稳定问题（Fletcher，1980）。

Rank-two 拟牛顿法

Rank-two拟牛顿法是Rank-one拟牛顿法的改进，其Hessian的更新公式为：

$$\mathbf H^{-1}_1 = \mathbf H^{-1}_0 + a \mathbf {uu}^T+b\mathbf {vv}^T, \tag{4.42}$$

其中，$N \times 1$矢量$\mathbf u$，$\mathbf u$，以及常数$a$、$b$满足4.34所示的拟牛顿公式。将4.42带入4.34可得：

$$(a \mathbf {uu}^T+b\mathbf {vv}^T) \Delta \mathbf g^{(k)} = \Delta \mathbf x^{(k)} - \mathbf H ^{-1}_k \Delta \mathbf g^{(k)}, \tag{4.43}$$

一组非常明显的解为：

$$\mathbf u = \Delta \mathbf x ^{(k)}, \tag{4.44}$$ $$\mathbf v = - \mathbf H ^{-1}_k \Delta \mathbf g^{(k)}, \tag{4.45}$$ $$a = \frac{1}{\mathbf u^T \Delta \mathbf g ^{(k)}} \tag{4.46}$$ $$b = \frac{1}{\mathbf v^T \Delta \mathbf g ^{(k)}} \tag{4.47}$$

则Rank-two拟牛顿法的更新公式为：

$$\mathbf H^{-1}_{k+1} =\mathbf H^{-1}_{k} +\frac{\Delta\mathbf x^{(k)} {\Delta\mathbf x^{(k)}}^T} { {\Delta\mathbf x ^{(k)}}^T \Delta\mathbf g^{(k)} } -\frac {\mathbf H^{-1}_k \Delta\mathbf g^{(k)} (\Delta\mathbf g^{(k)} \mathbf H^{-1}_k)^T} { {\Delta\mathbf g^{(k)}}^T \mathbf H^{-1}_k \Delta\mathbf g^{(k)} } \tag{4.48}$$

4.48即为Davidon-Fletcher-Powell（DFP）公式。结合确定性线搜索，DFP公式要优于最速下降法，运算效率要比CG更高（Fletcher，1980）。如果$\mathbf H^{-1}_0=\mathbf I$，DFP公式即变为CG法。DFP更新公式非常有效，但很快就被BFGS公式所超越，而BFGS公式目前被认为是所有准牛顿更新公式中最有效的（Nocedal & Wright, p.197, 1999）。
当未知数数量较多时，QN法中是较难存储$\mathbf H^{-1}$的。这种情况下可以使用有限内存的QN法（Nocedal, 1980; Zhu et al., 1997）,其使用了不超过十次迭代的位置信息和梯度信息，通过矢量-矢量运算来构建$\mathbf H_k$。

二次模型假设曲率（即$\mathbf H$的元素）在任何地方都是常数。在实际地震问题中的泛函都是高度非线性的，因此曲率值与空间中的位置相关。这意味着线性CG或QN方法（两者都假设纯二次函数）无法保证N次迭代收敛。一种典型的解决方法是具有重置的非线性CG或QN，在CG中，新的搜索方向沿着局部梯度方向。例如，在损失函数停止下降之前一直使用CG法，然后将搜索方向重置为最速下降方向，以开始新一轮CG迭代。同样也可以对QN方法进行类似的重置。