线性方程组的迭代方法

目录

直接方法与迭代方法

常规迭代算法

选择迭代求解器

预条件子

预条件子示例

均衡和重新排序

使用线性运算函数取代矩阵

---

        数值线性代数最重要也是最常见的应用之一是可求解以 A*x = b 形式表示的线性方程组。当 A 为大型稀疏矩阵时，您可以使用迭代方法求解线性方程组，使用这一方法，您可在计算的运行时间与解的精度之间进行权衡。本主题介绍 MATLAB® 中可用于求解方程 A*x = b 的迭代方法。

直接方法与迭代方法

求解线性方程 A*x = b 的方法有两种：

        直接方法是高斯消去法的变体。这些方法通过矩阵运算（例如 LU、QR 或 Cholesky 分解）直接使用各个矩阵元素。您可以使用直接方法来求解具有高精度水平的线性方程，但在大型稀疏矩阵上运行这些方法时，速度可能会很慢。用直接方法求解线性方程组的速度很大程度上取决于系数矩阵的密度和填充模式。

例如，以下代码用于求解一个较小的线性方程组。

A = magic(5);
b = sum(A,2);
x = A\b;
norm(A*x-b)
ans =

   1.4211e-14

        MATLAB 通过矩阵除法运算符 / 和 \ 以及 decomposition、lsqminnorm 和 linsolve 等函数实现直接的方法。

        迭代方法在经过有限的步数之后生成线性方程组的逼近解。这些方法对大型方程组非常有用，在求解这些方程组时，通过牺牲精度来缩短运行时间是合理的。迭代方法仅通过矩阵-向量积或抽象的线性运算函数间接使用系数矩阵。迭代方法可用于任何矩阵，但它们通常应用于直接求解速度慢的大型稀疏矩阵。使用间接方法求解线性方程组的速度不像直接方法那样严重依赖于系数矩阵的填充模式。但是，使用迭代方法通常需要针对每个特定问题调优参数。

        例如，以下代码用于求解一个具有对称正定系数矩阵的大型稀疏线性方程组。

A = delsq(numgrid('L',400));
b = ones(size(A,1),1);
x = pcg(A,b,[],1000);
norm(b-A*x)
pcg converged at iteration 796 to a solution with relative residual 9.9e-07.

ans =

   3.4285e-04

• MATLAB 根据系数矩阵 A 的属性，实现了多种具有不同优缺点的迭代方法。

如果提供实施这些方法所需的足够空间，直接方法与间接方法相比，通常前者更快，适用性更广。通常，您应该先尝试使用 x = A\b。如果直接求解的速度太慢，则可以尝试使用迭代方法。

常规迭代算法

求解线性方程组的大多数迭代算法都遵循类似的过程：

1. 首先获取解向量 x0 的初始估计值。（除非您指定了更好的估计值，否则通常为零元素向量。）

2. 计算残差范数 res = norm(b-A*x0)。

3. 将残差与指定的容差进行比较。如果 res <= tol，则结束计算并返回 x0 的计算答案。

4. 应用 A*x0 并根据残差值和其他计算出的量来更新向量 x0 的幅值和方向。这是完成大部分计算的步骤。

5. 重复步骤 2 到 4，直到 x0 的值足以满足容差要求。

        迭代方法在步骤 4 中更新 x0 的幅值和方向的方式各有不同，而且某些迭代方法在步骤 2 和 3 中的收敛条件也略有差别，以下概括了所有迭代求解器遵循的基本过程。

选择迭代求解器

        以下是 MATLAB 中的迭代求解器的流程图，大致概括了每种求解器适用的情况。一般而言，您可以对几乎所有的二次方非对称问题使用 gmres。在一些情况下，双共轭梯度算法（bicg、bicgstab、cgs 等）的效率高于 gmres，但它们的收敛行为难以预测，因此 gmres 往往是更好的初始选择。

如图所示：

预条件子

        迭代方法的收敛速度取决于系数矩阵的频谱（特征值）。因此，可以通过将线性方程组变换为具有更好的频谱（聚类特征值或接近 1 的条件数）来改善大多数迭代方法的收敛性和稳定性。这种变换的执行方法是将第二个被称为预条件子的矩阵应用于方程组。此过程将线性方程组

变换为等效方程组

        理想的预条件子将系数矩阵 A 变换为单位矩阵，因为任何迭代方法都将在使用这种预条件子的一次迭代中实现收敛。实际上，寻找一个好的预条件子需要权衡。采用以下三种方式之一执行变换：左侧预条件、右侧预条件或拆分预条件。

        第一种情况称为左侧预条件，因为预条件子矩阵 M 出现在 A 的左侧：

以下迭代求解器使用左侧预条件：

• ​bicg​

• ​gmres

• qmrqmr

        在右侧预条件中，M 出现在 A 的右侧：

以下迭代求解器使用右侧预条件：

• ​lsqr​

• bicgstabbicgstab

• bicgstablbicgstabl

• cgscgscgs
   

• tfqmr

        最后，对于对称系数矩阵A，拆分预条件可确保变换后的方程组仍是对称的。预条件子 M=HHT 被拆分，因子出现在 A 的不同侧：

​

        经过拆分预条件后的方程组的求解器算法以上述方程为基础，但实际上无需计算 H。求解器算法直接与 M 相乘并求解。

以下迭代求解器使用拆分预条件：

• pcg

• minres

• symmlq

        在所有情况下，选用预条件子 M 以加快迭代方法的收敛。当迭代解的残差停滞不前或两次迭代之间进展甚微时，通常意味着您需要生成一个预条件子矩阵以加入问题中。

        MATLAB 中的迭代求解器允许您指定单个预条件子矩阵 M，或两个预条件子矩阵因子，使得 M = M1M2。这使得以分解形式指定预条件子变得容易，例如 M = LU。请注意，在同时保留了 M = HHT 的拆分预条件情况下，M1 和 M2 输入与 H 因子之间没有关联。

        在某些情况下，预条件子在给定问题的数学模型中自然发生。在没有自然预条件子的情况下，可以使用下表中的不完全分解之一来生成预条件子矩阵。不完全分解本质上是不完全直接求解，计算起来很快。

预条件子示例

        以二维方形域中的拉普拉斯方程的五点有限差分逼近方程为例。以下命令使用预条件共轭梯度法 (PCG) 和预条件子 M = L*L'，其中 L 是 A 的零填充不完全 Cholesky 因子。对于此方程组，pcg 无法在不指定预条件子矩阵的情况下求得解。

A = delsq(numgrid('S',250));
b = ones(size(A,1),1);
tol = 1e-3;
maxit = 100;
L = ichol(A);
x = pcg(A,b,tol,maxit,L,L');



pcg converged at iteration 92 to a solution with relative residual 0.00076.

        pcg 需要 92 次迭代才能达到指定的容差。但是，使用其他预条件子可以产生更好的结果。例如，使用 ichol 构造修正的不完全 Cholesky，pcg 只需经过 39 次迭代便可达到指定的容差。

L = ichol(A,struct('type','nofill','michol','on'));
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 39 to a solution with relative residual 0.00098.

均衡和重新排序

        对于难以计算的问题，您可能需要比 ilu 或 ichol 直接生成的预条件子更好的预条件子。例如，可能希望生成质量更好的预条件子，或最小化要执行的计算量。在这些情况下，可以使用均衡来使系数矩阵更加对角占优（这可以生成质量更好的预条件子），并使用重新排序来最小化矩阵因子中的非零元素数量（这可以降低内存需求并可能提高后续计算的效率）。

如果同时使用均衡和重新排序来生成预条件子，则该过程为：

1. ​对系数矩阵使用 equilibrate。

2. ​使用稀疏矩阵重新排序函数（例如 dissect 或 symrcm）对均衡后的矩阵重新排序。

3. ​使用 ilu 或 ichol 生成最终预条件子。

以下示例使用均衡和重新排序生成了用于稀疏系数矩阵的预条件子。

        创建系数矩阵 A 和全 1 向量 b 作为线性方程的右侧。计算 A 的条件数估计值。

load west0479;
A = west0479;
b = ones(size(A,1),1);
condest(A)
ans =

   1.4244e+12

使用 equilibrate 来改善系数矩阵的条件数。

[P,R,C] = equilibrate(A);
Anew = R*P*A*C;
bnew = R*P*b;
condest(Anew)
ans =

   5.1042e+04

        使用 dissect 对均衡后的矩阵重新排序。

q = dissect(Anew);
Anew = Anew(q,q);
bnew = bnew(q);

使用不完全 LU 分解生成预条件子。

[L,U] = ilu(Anew);

        使用 gmres，通过预条件子矩阵、容差 1e-10、最大外部迭代次数 50 和内部迭代次数 30，求解线性方程组。

tol = 1e-10;
maxit = 50;
restart = 30;
[xnew, flag, relres] = gmres(Anew,bnew,restart,tol,maxit,L,U);
x(q) = xnew;
x = C*x(:);

        现在，将 gmres（其中包括预条件子）返回的 relres 相对残差与没有使用预条件子的相对残差 resnew 和没有经过均衡的相对残差 res 进行比较。结果表明，尽管线性方程组都是等效的，但不同的方法会将不同的权重应用于每个元素，而这可能显著影响残差值。

relres
resnew = norm(Anew*xnew - bnew) / norm(bnew)
res = norm(A*x - b) / norm(b)
relres =
   8.7537e-11
resnew =
   3.6805e-08
res =
   5.1415e-04

使用线性运算函数取代矩阵

        MATLAB 中的迭代求解器不要求为 A 提供数值矩阵。由于求解器执行的计算使用矩阵-向量乘法 A*x 或 A'*x 的结果，因此您可以改而提供一个函数来计算这些线性运算的结果。计算这些量的函数通常被称为线性运算函数。

        除了使用线性运算函数取代系数矩阵 A，还可以使用线性运算函数取代预条件子矩阵 M。在这种情况下，函数需要计算 M\x 或 M'\x，如求解器参考页所示。

        通过使用线性运算函数，您可以利用 A 或 M 中的模式来计算线性运算的值，这比求解器显式使用矩阵执行满矩阵-向量乘法的效率更高。这也意味着，不需要内存来存储系数矩阵或预条件子矩阵，因为线性运算函数通常计算矩阵-向量乘法的结果而根本不用构建矩阵。

例如，有以下系数矩阵

A = [2 -1  0  0  0  0;
    -1  2 -1  0  0  0;
     0 -1  2 -1  0  0;
     0  0 -1  2 -1  0;
     0  0  0 -1  2 -1;
     0  0  0  0 -1  2];

        当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。因此，对于给定的向量 x，线性运算函数只需要将三个向量相加便可计算 A*x 的值：

function y = linearOperatorA(x)
y = -1*[0; x(1:end-1)] ...
  + 2*x ...
  + -1*[x(2:end); 0];
end

        大多数迭代求解器要求线性运算函数基于 A 返回 A*x 的值。同样，对于预条件子矩阵 M，该函数通常必须计算 M\x。对于求解器 lsqr、qmr 和 bicg，如果请求了返回 A'*x 或 M'\x 的值，线性运算函数还必须返回这些值。

参考

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

[2] Saad, Yousef, Iterative Methods for Sparse Linear Equations. PWS Publishing Company, 1996.