迭代技术 (GNU Octave (版本 11.1.0))

下一个: 使用稀疏矩阵的实际示例, 上一个: 稀疏矩阵的线性代数, 上级: 稀疏矩阵 [目录][索引]

22.3 应用于稀疏矩阵的迭代技术 ¶

上一节讨论的左除 \ 和右除 / 运算符使用直接求解器来求解形如 x = A \ b 或 x = b / A 的线性方程。Octave 还包含许多使用迭代技术求解稀疏线性方程的函数。

x = pcg (A, b, tol, maxit, m1, m2, x0, …) ¶

x = pcg (A, b, tol, maxit, M, [], x0, …) ¶

[x, flag, relres, iter, resvec, eigest] = pcg (A, b, …) ¶

通过预处理共轭梯度（Preconditioned Conjugate Gradient）迭代法求解线性方程组 A * x = b。

输入参数如下：

A 是线性系统的矩阵，且必须是方阵。A 可以作为矩阵、函数句柄或内联函数 Afcn 传入，满足 Afcn(x) = A * x。传递给 Afcn 的附加参数可以在 x0 之后传入。
A 必须是 Hermitian 正定矩阵（HPD）。如果 pcg 检测到 A 不是正定的，则会打印一条警告并设置 flag 输出。
b 是右侧向量。
tol 是残差误差所需的相对容差，即 b - A * x。当 norm (b - A * x) ≤ tol * norm (b) 时迭代停止。如果 tol 省略或为空，则使用默认容差 1e-6。
maxit 是允许的最大迭代次数；如果 maxit 省略或为空，则使用默认值 20。
m 是一个 HPD 预处理矩阵。对于任意分解 m = p1 * p2，使得 inv (p1) * A * inv (p2) 是 HPD 的，则共轭梯度法正式应用于线性系统 inv (p1) * A * inv (p2) * y = inv (p1) * b，其中 x = inv (p2) * y（分裂预处理）。实际上，在共轭梯度法的每次迭代中，都会使用 mldivide 求解一个以 m 为矩阵的线性系统。如果存在特定的分解 m = m1 * m2（例如，a 的不完全 Cholesky 分解），则可以传入两个矩阵 m1 和 m2，相应的线性系统将使用 mldivide 运算符求解。请注意，正确选择预处理器可以显著提高方法的整体性能。除了传入矩阵 m1 和 m2，用户也可以传入两个函数，分别返回将 m1 和 m2 的逆应用于向量的结果。如果 m1 被省略或为空 []，则不应用任何预处理。如果 m 没有可用的分解，可以将 m2 省略或留空，并使用输入变量 m1 来传递预处理器 m。
x0 是初始猜测值。如果 x0 省略或为空，函数默认将 x0 设为零向量。

在 x0 之后的参数被视为附加参数，并以适当方式传递给任何已提供给 pcg 的函数（A 或 m1 或 m2）。详见下面的示例。

输出参数如下：

x 是对 A * x = b 解的近似计算结果。如果算法未收敛，则 x 是具有最小残差的那次迭代结果。
flag 报告收敛情况：
- 0：算法在规定的容差内收敛。
- 1：算法未收敛并达到了最大迭代次数。
- 2：预处理矩阵是奇异的。
- 3：算法停滞，即当前迭代 x 与前一次迭代的差的绝对值小于 eps * norm (x,2)。
- 4：算法检测到输入（经过预处理的）矩阵不是 HPD。
relres 是最终残差与其初始值的比值，以欧几里得范数衡量。
iter 表示计算出 x 时所处的迭代次数。由于输出 x 对应于最小残差解，因此该方法执行的总迭代次数由 length(resvec) - 1 给出。
resvec 描述方法的收敛历史。resvec (i, 1) 是残差的欧几里得范数，resvec (i, 2) 是预处理后的残差范数，均在第 (i-1) 次迭代之后，其中 i = 1, 2, …, iter+1。预处理残差范数定义为 r' * (m \ r)，其中 r = b - A * x，另请参见 m 的描述。如果不需要 eigest，则仅返回 resvec (:, 1)。
eigest 返回预处理矩阵 P = m \ A 的最小特征值 eigest(1) 和最大特征值 eigest(2) 的估计值。特别地，如果未使用预处理，则返回 A 的极端特征值估计。eigest(1) 是高估值，eigest(2) 是低估值的，因此 eigest(2) / eigest(1) 是 cond (P, 2) 的下界，理论上极限情况下应等于条件数的实际值。

让我们考虑一个三对角矩阵的简单问题

n = 10;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
b = A * ones (n, 1);
M1 = ichol (A); # 在此三对角情形下等同于 chol (A)'
M2 = M1';
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

示例 1：pcg 的最简单用法

x = pcg (A, b)

示例 2：使用计算 A * x 的函数调用 pcg

x = pcg (Afcn, b)

示例 3：使用预处理矩阵 M 调用 pcg

x = pcg (A, b, 1e-06, 100, M)

示例 4：使用函数作为预处理器调用 pcg

x = pcg (Afcn, b, 1e-6, 100, Mfcn)

示例 5：使用预处理矩阵 M1 和 M2 调用 pcg

x = pcg (A, b, 1e-6, 100, M1, M2)

示例 6：使用函数作为预处理器调用 pcg

x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn)

示例 7：输入需要参数的函数调用 pcg

  function y = Ap (A, x, p) # 计算 A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = pcg (Apfcn, b, [], [], [], [], [], 2);

示例 8：展示 pcg 使用分裂预处理器的显式示例

M1 = ichol (A + 0.1 * eye (n)); # 扰动后的 A 的分解
M2 = M1';
M = M1 * M2;

## pcg 两次迭代后计算的参考解
[x_ref, fl] = pcg (A, b, [], 2, M)

## 分裂预处理
[y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
x = M2 \ y # 比较 x 和 x_ref

参考文献：

C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995.（基础 PCG 算法）
Y. Saad, Iterative Methods for Sparse Linear Systems, PWS, 1996.（来自 PCG 的条件数估计）本书的修订版可在线获取： https://www-users.cs.umn.edu/~saad/books.html

另请参阅： sparse, pcr, gmres, bicg, bicgstab, cgs.

x = pcr (A, b, tol, maxit, m, x0, …) ¶

[x, flag, relres, iter, resvec] = pcr (…) ¶

通过预处理共轭残差（Preconditioned Conjugate Residuals）迭代法求解线性方程组 A * x = b。

输入参数如下：

A 可以是方阵（最好是稀疏矩阵）、函数句柄、内联函数或包含计算 A * x 的函数名称的字符串。原则上 A 应该是对称且非奇异的；如果 pcr 发现 A 在数值上是奇异的，您将收到一条警告消息，并且 flag 输出参数将被设置。
b 是右侧向量。
tol 是残差误差所需的相对容差，即 b - A * x。当 norm (b - A * x) <= tol * norm (b - A * x0) 时迭代停止。如果 tol 为空或被省略，函数默认设置 tol = 1e-6。
maxit 是允许的最大迭代次数；如果为 maxit 提供了 []，或者 pcr 的参数较少，则使用默认值 20。
m 是（左）预处理矩阵，因此迭代（理论上）等价于通过 pcr 求解 P * x = m \ b，其中 P = m \ A。请注意，正确选择预处理器可以显著提高方法的整体性能。除了矩阵 m，用户也可以传入一个函数，该函数返回将 m 的逆应用于向量的结果（这通常是使用预处理器的首选方式）。如果为 m 提供了 [] 或省略了 m，则不应用预处理。
x0 是初始猜测值。如果 x0 为空或被省略，函数默认设置 x0 为零向量。

在 x0 之后的参数被视为附加参数，并以适当方式传递给任何已提供给 pcr 的函数（A 或 m）。详见下面的示例。

输出参数如下：

x 是对 A * x = b 解的近似计算结果。
flag 报告收敛情况。flag = 0 表示解已收敛且满足由 tol 给出的容差。flag = 1 表示达到了 maxit 的迭代次数限制。flag = 3 表示 pcr 分解，详见 [1]。
relres 是最终残差与其初始值的比值，以欧几里得范数衡量。
iter 是实际执行的迭代次数。
resvec 描述方法的收敛历史，resvec (i) 包含第 (i-1) 次迭代后残差的欧几里得范数，i = 1,2, …, iter+1。

让我们考虑一个对角矩阵的简单问题（我们利用 A 的稀疏性）

n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);

示例 1：pcr 的最简单用法

x = pcr (A, b)

示例 2：使用计算 A * x 的函数调用 pcr

function y = apply_a (x)
  y = [1:10]' .* x;
endfunction

x = pcr ("apply_a", b)

示例 3：带完整诊断的预处理迭代。预处理器（有些特别，因为即使原始矩阵 A 很简单）被定义为一个函数

function y = apply_m (x)
  k = floor (length (x) - 2);
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m")
semilogy ([1:iter+1], resvec);

示例 4：最后，一个依赖参数 k 的预处理器

function y = apply_m (x, varargin)
  k = varargin{1};
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m"', [], 3)

参考文献：

W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, section 9.5.4, Springer, 1994.

另请参阅： sparse, pcg.

使用预处理矩阵 M 可以加快迭代求解器收敛到解的速度。在这种情况下，转而求解线性方程 M^-1 * x = M^-1 * A \ b。典型的预处理矩阵是原始矩阵的部分分解。

L = ichol (A) ¶

L = ichol (A, opts) ¶

计算稀疏方阵 A 的不完全 Cholesky 分解。

默认情况下，ichol 仅使用 A 的下三角部分，并返回一个下三角因子 L，使得 L*L' 近似于 A。

该例程给出的因子可用作通过迭代方法（如 PCG（预处理共轭梯度法））求解的线性方程组的预处理器。

可以通过在结构体 opts 中传递选项来修改分解。选项名称是结构体的字段，设置是字段的值。名称和取值区分大小写。

type

分解的类型。

"nofill"（默认）

无填充的不完全 Cholesky 分解（IC(0)）。

额外支持的选项：michol，shape。

"ict"

带阈值丢弃的不完全 Cholesky 分解（ICT）。

额外支持的选项：droptol，michol，shape，udiag。

diagcomp

一个非负标量 alpha，用于对 A + alpha * diag (diag (A)) 而非 A 进行不完全 Cholesky 分解。当 A 不是正定矩阵时，这很有用。默认值为 0。

droptol

一个非负标量，指定 ICT 分解的丢弃容差。默认值为 0，将产生完整的 Cholesky 分解。

L 的非对角元设置为 0，除非

abs (L(i,j)) >= droptol * norm (A(j:end, j), 1)。

michol

修正的不完全 Cholesky 分解：

"off"（默认）: 行和与列和不一定保留。
"on": L 的对角元被修改，以便即使在分解过程中丢弃了元素，也能保留行（和列）和。保留的关系是：A * e = L * L' * e，其中 e 是全 1 向量。

shape

"lower"（默认）: 仅使用 A 的下三角部分，并返回一个下三角因子 L，使得 L*L' 近似于 A。
"upper": 仅使用 A 的上三角部分，并返回一个上三角因子 U，使得 U'*U 近似于 A。

示例

下面的问题演示了如何使用完全 Cholesky 分解和不完全 Cholesky 分解来分解一个对称正定样本矩阵。

A = [ 0.37, -0.05,  -0.05,  -0.07;
     -0.05,  0.116,  0.0,   -0.05;
     -0.05,  0.0,    0.116, -0.05;
     -0.07, -0.05,  -0.05,   0.202];
A = sparse (A);
nnz (tril (A))
ans =  9
L = chol (A, "lower");
nnz (L)
ans =  10
norm (A - L * L', "fro") / norm (A, "fro")
ans =  1.1993e-16
opts.type = "nofill";
L = ichol (A, opts);
nnz (L)
ans =  9
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.019736

另一个分解示例是用于求解单位正方形上边值问题的有限差分矩阵。

nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
               [-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
               [-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans =  239400
opts.type = "nofill";
L = ichol (A, opts);
nnz (tril (A))
ans =  239400
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.062327

已实现算法的参考文献：

[1] Y. Saad, "Preconditioning Techniques", Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996.

[2] M. Jones, P. Plassmann, An Improved Incomplete Cholesky Factorization, 1992.

另请参阅： chol, ilu, pcg.

LUA = ilu (A) ¶

LUA = ilu (A, opts) ¶

[L, U] = ilu (…) ¶

[L, U, P] = ilu (…) ¶

计算稀疏方阵 A 的不完全 LU 分解。

ilu 返回一个单位下三角矩阵 L、一个上三角矩阵 U，以及可选的置换矩阵 P，使得 L*U 近似于 P*A。

该例程给出的因子可用作通过迭代方法（如 BICG（双共轭梯度法）或 GMRES（广义最小残差法））求解的线性方程组的预处理器。

可以通过在结构体 opts 中传递选项来修改分解。选项名称是结构体的字段，设置是字段的值。名称和取值区分大小写。

type

分解的类型。

"nofill"（默认）

无填充的 ILU 分解（ILU(0)）。

额外支持的选项：milu。

"crout"

Crout 版本的 ILU 分解（ILUC）。

额外支持的选项：milu，droptol。

"ilutp"

带阈值和主元选择的 ILU 分解。

额外支持的选项：milu，droptol，udiag，thresh。

droptol

一个非负标量，指定分解的丢弃容差。默认值为 0，将产生完整的 LU 分解。

U 的非对角元设置为 0，除非

abs (U(i,j)) >= droptol * norm (A(:,j))。

L 的非对角元设置为 0，除非

abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j)。

milu

修正的不完全 LU 分解：

"row": 行和修正的不完全 LU 分解。分解保留行和：A * e = L * U * e，其中 e 是全 1 向量。
"col": 列和修正的不完全 LU 分解。分解保留列和：e' * A = e' * L * U。
"off"（默认）: 行和与列和不一定保留。

udiag

如果为 true，则上三角因子中对角线上的任何零将被替换为局部丢弃容差 droptol * norm (A(:,j))/U(j,j)。默认值为 false。

thresh

分解的主元阈值。其范围可在 0（对角主元）到 1（默认值）之间，其中选择列中幅度最大的条目作为主元。

如果 ilu 仅以一个输出参数调用，则返回的矩阵是 L + U - speye (size (A))，其中 L 是单位下三角矩阵，U 是上三角矩阵。

对于两个输出参数，ilu 返回一个单位下三角矩阵 L 和一个上三角矩阵 U。对于 opts.type == "ilutp"，其中一个因子会根据 opts.milu 的值进行置换。当 opts.milu == "row" 时，U 是列置换的上三角因子。否则，L 是行置换的单位下三角因子。

如果有三个命名输出且 opts.milu != "row"，则返回 P，使得 L 和 U 是 P*A 的不完全因子。当 opts.milu == "row" 时，返回 P，使得 L 和 U 是 A*P 的不完全因子。

示例

A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840

nnz (lu (A))
ans = 126478

nnz (ilu (A, opts))
ans = 7840

这表明 A 有 7,840 个非零元素，完全 LU 分解有 126,478 个非零元素，而不完全 LU 分解（填充级别为 0）有 7,840 个非零元素，与 A 相同。摘自： https://www.mathworks.com/help/matlab/ref/ilu.html

A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)

此示例使用 ILU 作为随机 FEM 矩阵的预处理器，该矩阵具有很大的条件数。如果没有 L 和 U，BICG 将无法收敛。

另请参阅： lu, ichol, bicg, gmres.