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

下一个使用稀疏矩阵的真实生活示例前一个稀疏矩阵上的线性代数上级稀疏矩阵[目录][索引]

22.3稀疏矩阵的迭代技术¶

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

: 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, …) ¶

求解线性方程组A * x = b 采用预处理共轭梯度迭代法。

输入参数为:

A是线性系统的矩阵,它必须是平方。A可以作为矩阵、函数句柄或内联函数传递Afcn使得Afcn(x) = A * x。的附加参数Afcn可能在之后通过x0.
A必须是埃尔米特和正定(HPD)。如果pcg检测A不是肯定的,打印一个警告flag设置输出。
b是右侧向量。
tol是残余误差所需的相对误差范围,b - A * x 。如果标准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(分裂预处理)。在实践中,在共轭梯度法的每次迭代中,都使系统与矩阵等距m用解决mldivide.如果一个特定的因子分解m = m1 * m2是可用的(例如,的完全Cholesky因式分解a),两个矩阵m1和m2可以通过,并且相对线性系统用mldivide运算符请注意,正确选择预处理器可以显著提高该方法的整体性能。而不是矩阵m1和m2,用户可以传递两个函数,返回应用的逆运算的结果m1和m2到向量。如果m1被省略或为空[],则不应用任何预处理。如果没有的因子分解m是可用的,m2可以省略或保留[],并且输入变量m1可以用来通过预处理器m.
x0是最初的猜测。如果x0被省略或为空,则函数集x0默认情况下为零向量。

以下参数x0被视为参数,并以适当的方式传递给任何函数(A或m1或m2)已经给予pcg。有关更多详细信息,详见下面的示例。

输出参数为:

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返回最小值的估计值eigest(1)和最大的eigest(2)预处理矩阵的特征值P = m \ A 特别地,如果使用norepresenting,则的极端特征值的估计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); # in this tridiagonal case it corresponds to 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: pcg具有一个计算A * x

x = pcg (Afcn, b)

示例3: pcg具有预处理矩阵M

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

示例4: pcg具有作为预处理器的函数

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

示例5: pcg具有预处理矩阵M1和M2

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) # compute 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使用asplit预处理器

M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
M2 = M1';
M = M1 * M2;

## reference solution computed by pcg after two iterations
[x_ref, fl] = pcg (A, b, [], 2, M)

## split preconditioning
[y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
x = M2 \ y # compare x and x_ref

参考文献:

C.T.Kelley,Iterative Methods for Linear and Nonlinear Equations,暹罗,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.

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);

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)

使用预处理矩阵可以加快迭代解算器收敛到解的速度M.在这种情况下,线性方程

M^-1 * x = M^-1 *
A \ b

而是被解决。典型的预处理矩阵是原始矩阵的部分因子分解。

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