2.3　常用经典数值解法

在数字图像处理中存在有大量的数值计算问题，例如求函数极值或求解线性方程组等等。但这些问题在计算机中的求解和我们平常在数学课上的处理方式是不同的。例如，很多时候，在利用计算机求解某个问题时，我们并不要求得出一个精确解，而是要求一个近似解。毕竟计算机所能表示的数值精度有限，很多精确解，计算机也无法求出或表示出来。计算机所擅长的（其实也是它仅能做的）就是大量地重复执行简单的任务。所以计算机中常用的经典数值解法本质上来说都是采用迭代求解的策略。

2.3.1　牛顿迭代法

牛顿迭代法（Newton method）又称为牛顿-拉夫逊方法（Newton-Raphson method），它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。

既然牛顿迭代法可以用来求解方程的根，那么不妨以方程x²=n为例，试着求解它的根。为此令f（x）=x²-n，也就相当于是求f（x）=0的解，如图2-2所示。首先随便找一个初始值x₀，如果x₀不是解，做一个经过（x₀，f（x₀））这个点的切线，与x轴的交点为x₁。同样道理，如果x₁不是解，做一个经过（x₁，f（x₁））这个点的切线，与x轴的交点为x₂。以此类推，不断迭代，所得到之x_i会无限趋近于f（x）=0的解。

判断x_i是否是f（x）=0的解有两种方法：一是直接计算f（x_i）的值判断是否为0，二是判断前后两个解x_i和x_i-1是否无限接近。

经过（x_i，f（x_i））这个点的切线方程为（注意这也是一元函数的一阶泰勒展式）

f（x）=f（x_i）+f′（x_i）（x-x_i）

其中，f′（x）为f（x）的导数，例如在当前所讨论的问题中为2x。令切线方程等于0，即可求出

于是乎就得到了一个迭代公式，而且它必然在f（x^∗）=0处收敛，其中x^∗就是方程的根，由此便可对方程进行迭代求根。

例如当前所讨论的问题，继续化简就会得到迭代公式为

基于上述迭代公式，我们其实给出了一个求平方根的算法。事实上，这也的确是很多计算机语言中内置的开平方函数的实现方法。

有时候方程的求根公式可能很复杂（甚至某些方程可能没有求根公式），导致求解困难。这时便可利用牛顿法进行迭代求解。上面的例子已经很好地演示了利用牛顿迭代法解方法的原理。但事实上，它的用途还不仅限于解方程。对上述方法稍加调整，就能得出一个强有力的优化工具。

假设当前任务是优化一个目标函数f，也就是求该函数的极大值或极小值问题，可以转化为求解函数f的导数f′=0的问题，这样求可以把优化问题看成方程f′=0求解问题。剩下的问题就和前面提到的牛顿迭代法解方程的过程很相似了。这次为了求解方程f′=0的根，把原函数f（x）的做泰勒展开，展开到二阶形式（注意之前是一阶）：

当且仅当Δx无线趋近于0时，（可以舍得后面的无穷小项）使得等式成立。此时上式等价于

注意因为Δx无限趋近于0，前面的常数1/2将不再起作用，可以将其一并忽略，即

求解得

也就可以得出迭代公式

事实上，最优化问题除了用牛顿法来解之外，还可以用梯度下降法来解。但是通常来说，牛顿法可以利用到曲线本身的信息，比梯度下降法更容易收敛，即迭代更少次数。

回想第1章中介绍的黑塞矩阵与多元函数极值问题，对于一个多维向量x，以及在点x₀的邻域内有连续二阶偏导数的多元函数f（x），可以写出该函数在点x₀处的（二阶）泰勒展开式

其中，是高阶无穷小表示的皮亚诺余项。而Hf（x₀）是一个黑塞矩阵。依据之前的思路，忽略掉无穷小项，写出迭代公式即为

由此可见，高维情况依然可以用牛顿迭代求解，但问题是黑塞矩阵引入的复杂性，使得牛顿迭代求解的难度大大增加。所以人们又提出了所谓的拟牛顿法（Quasi-Newton method），即不再直接计算黑塞矩阵，而是每一步的时候使用梯度向量更新黑塞矩阵的近似。拟牛顿法是针对牛顿法的弱点进行了改进，因而更具实际应用价值。限于篇幅，我们并不打算在此处展开，有兴趣的读者可以参阅相关资料以了解更多关于拟牛顿法的细节。

2.3.2　雅可比迭代

雅可比迭代法（Jacobi iterative method）是一种常用的求解线性方程组的方法。作为一个简单的示例，考虑如下来自文献[7]的方程组

上述方程可表示成如下形式

这样就提出了下列雅可比迭代过程

如果从P₀=（x₀，y₀，z₀）=（1，2，2）开始，则上式中的迭代将收敛到解（2，4，3）。

将x₀=1，y₀=2和z₀=2代入上式中每个方程的右边，即可得到如下新值

新的点P₁=（1.75，3.375，3.00）比P₀更接近（2，4，3）。使用上述雅可比迭代过程生成点的序列{P_k}将收敛到解（2，4，3）。

这个过程称为雅可比迭代，可用来求解某些类型的线性方程组。从表2-1中可以看出，经过19步迭代，迭代过程收敛到一个精度为9位有效数字的近似值（2.00000000，4.00000000，3.00000000）。但有时雅可比迭代法是无效的。通过下面的例子可以看出，重新排列初始线性方程组后，应用雅可比迭代法可能会产生一个发散的点的序列。

设重新排列的线性方程组如下

这些方程可以表示为如下形式

用如下雅可比迭代过程求解

如果同样P₀=（x₀，y₀，z₀）=（1，2，2）开始，则上式中的迭代将对解（2，4，3）发散。将x₀=1，y₀=2和z₀=2代入上式中每个方程的右边，即可得到新值x₁，y₁和z₁

新的点P₁=（-1.5，3.375，5.00）比P₀更远地偏离（2，4，3）。如表2-2所示，使用上述迭代过程生成点的序列是发散的。

2.3.3　高斯迭代法

有时候可以通过一些技巧来加快常规的雅可比迭代的收敛速度。观察由雅可比迭代过程产生的3个序列{x_k}，{y_k}和{z_k}，它们分别收敛到2.00，4.00和3.00。由于x_k+1被认为是比x_k更好的x的近似值，所以在计算y_k+1时用x_k+1来替换x_k是合理的。同理，可用x_k+1和y_k+1计算z_k+1。这种方法被称为高斯-赛德尔（Gauss-Seidel）迭代法，简称高斯迭代法。

下面就以2.3.2节中给出的方程组为例演示高斯迭代法的具体执行步骤。首先，写出下列高斯迭代过程

如果从P₀=（x₀，y₀，z₀）=（1，2，2）开始，用上式中的迭代可收敛到解（2，4，3）。

将y₀=2和z₀=2代入上式中第一个方程可得

将x₁=1.75和z₀=2代入第二个方程可得

将x₁=1.75和y₁=3.75代入第三个方程可得

新的点P₁=（1.75，3.75，2.95）比P₀更接近解（2，4，3），而且比2.3.2节例子中的值更好。用高斯迭代过程生成序列{P_k}收敛到（2，4，3）。

正如前面讨论的，应用雅可比迭代法计算有时可能是发散的，所以有必要建立一些判定条件判断雅可比迭代是否收敛。在给出这个条件之前，先来看看严格对角占优矩阵的定义。设有N×N维矩阵A，如果

其中，i是行号，j是列号，则称该矩阵是严格对角占优矩阵。显然，严格对角占优的意思就是指对角线上元素的绝对值不小于所在行其他元素的绝对值和。

现在使雅可比迭代和高斯迭代过程一般化。设有如下线性方程组

设第k个点，则下一点。向量P_k的上标（k）可用来标识属于该点的坐标。迭代公式根据前面的值（，使用上述线性方程组中第j行求解x^（k+1）_j。

雅可比迭代

其中，j=1，2，…，N。

雅可比迭代使用所有旧坐标来生成所有新坐标，而高斯-赛德尔迭代尽可能使用新坐标得到更新的坐标。

高斯-赛德尔迭代

其中，j=1，2，…，N。

下面的定理给出了雅可比迭代收敛的充分条件：设矩阵A具有严格对角优势，则AX=B有唯一解X=P。利用前面给出的迭代式可产生一个向量序列{P_k}，而且对于任意初始向量{P₀}，向量序列都将收敛到P。

当矩阵A具有严格对角优势时，可证明高斯迭代法也会收敛。在大多数情况下，高斯迭代法比雅可比迭代法收敛得更快，因此通常会更倾向于使用高斯迭代法。但在某些情况下，雅可比迭代会收敛，而高斯迭代不会收敛。

2.3.4　托马斯算法

当线性联立方程组的系数矩阵为三对角阵时，可采用非常有效的托马斯（Thomas）算法来进行求解。为了帮助读者理解这种数值计算法方法，下面以一个6×6矩阵为例来介绍，例如有如下的一个线性方程组

将矩阵变为上三角矩阵，首先要把上面公式中的系数矩阵变为一个上三角矩阵。

对于第1行

b₁x₁+c₁x₂=d₁

将上式除以b₁，得

可写成

对于第2行

a₂x₁+b₂x₂+c₂x₃=d₂

将变换后的第1行乘以a₂，再与第2行相减，可消去x₁，得

（b₂-a₂γ₁）x₂+c₂x₃=d₂-a₂ρ₁

所以新的矩阵方程为

同理，可推得第3行至第6行分别为

最后得到新的上三角矩阵公式为

接下来就可以逆序求出结果，如

因此可以归纳总结出一般性公式，对于如下的线性方程组

注意使用托马斯算法求解时，系数矩阵需要是对角占优的。另外，从以上讨论可以看出，托马斯算法具有线性算法复杂性O（N），而高斯消元法的复杂性为O（N²）。还可以看到，它虽然含有对矩阵元素的运算，但并不需要以矩阵形式存储A，只需以矢量形式存储三对角元素即可。托马斯算法也常称为追赶法。因为前向代入过程可以称为“追”，反向代入可称为“赶”。