迭代技术(GNU Octave(版本11.1.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, SIAM, 1995. (the base PCG algorithm)
Y. Saad, Iterative Methods for Sparse Linear Systems, PWS, 1996. (condition number estimate from PCG) Revised version of this book is available online at 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 (…) ¶

求解线性方程组A * x = b利用预处理共轭残差迭代法。

输入参数为

A可以是平方(最好是稀疏)矩阵,也可以是函数句柄、内联函数或包含计算函数名称的字符串A * x原则上A应该是对称的和非奇异的;如果pcr查找A为了在数字上保持棱角分明,您将get一条警告消息和flag将设置outputparameter。
b是右侧向量。
tol是残余误差所需的相对误差范围,b - A * x。如果标准b - A * x) <=tol标准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被视为参数,并以适当的方式传递给任何函数(A或m)传递给pcr。有关更多详细信息,详见下面的示例。

输出参数为

x是的解的计算近似值A * x = b.
flag关于趋同的返回。flag = 0表示收敛的解和下式给出的误差范围标准tol令人满意。flag = 1意味着maxit达到了盈利计数的极限。flag = 3返回apcr细分,详见[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: pcr具有一个计算A * x.

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

x = pcr ("apply_a", b)

示例3:预处理迭代,具有完整diag函数。预处理器(很奇怪,因为即使是原始矩阵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) ¶

计算稀疏平方矩阵的不完全Cholesky因子分解A.

默认情况下,ichol只使用的下三角矩阵A并返回一个较低的三角因子L使得L*L'近似A.

该子程序给出的因子可以用作通过迭代方法(如PCG(预处理共轭梯度))求解的线性方程系统的预处理器。

因子分解可以通过传递结构体中的参数来修改opts。参数名称是结构体的一个字段,设置是字段的值。键和值区分大小写。

type

因子分解的类型。

"nofill" (default): 没有填充的不完全Cholesky因子分解(IC(0))。
"ict": 具有阈值下降的不完全Cholesky因子分解(ICT)。

diagcomp

非负标量alpha的不完全Cholesky因子分解A + alpha * diag (diag (A))而不是A。当A不是肯定的。默认值为0。

droptol

一种非负标量,指定在执行ICT时因式分解的下降容限。默认值为0,生成完整的Cholesky因子分解。

的非对角矩阵分量L设置为0,除非

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

michol

改进的不完全Cholesky因子分解:

"off" (default): 行和列的总和不一定保留。
"on": 的对角矩阵L被修改为即使在因子分解过程中删除了元素,也可以保留行(和列)和。保留的关系是:A * e = L * L' * e,其中e是一个1的向量。

shape

"lower" (default): 只使用的下三角矩阵A并返回一个较低的三角因子L使得L*L'近似A.
"upper": 仅使用的上三角矩阵A并返回一个上三角因子U使得U'*U近似A.

示例

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

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 (…) ¶

计算稀疏平方矩阵的不完全LU因子分解A.

ilu返回单位下三角矩阵L,上三角矩阵U,以及可选的置换矩阵P,使得L*U近似P*A.

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

因子分解可以通过传递结构体中的参数来修改opts。参数名称是结构体的一个字段,设置是字段的值。键和值区分大小写。

type

因子分解的类型。

"nofill" (default)

没有填充的ILU因子分解(ILU(0))。

其他支持的参数:milu.

"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" (default): 行和列的总和不一定保留。

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类型"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具有7840个非零,完全LU因子分解具有126478个非零;不完全LU因子因子分解具有0级填充,具有7840非零,与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和UBICG不会进行任何转换。

详见: lu, ichol, bicg, gmres.