MATLAB中gmres函数用法

目录

语法

说明

示例

线性系统的迭代解

使用 gmres 预条件子而不重启

使用 gmres 预条件子并重启

提供初始估计值

使用函数句柄代替数值矩阵

---

        gmres函数的功能是求解线性系统 - 广义最小残差法。

语法

x = gmres(A,b)
x = gmres(A,b,restart)
x = gmres(A,b,restart,tol)
x = gmres(A,b,restart,tol,maxit)
x = gmres(A,b,restart,tol,maxit,M)
x = gmres(A,b,restart,tol,maxit,M1,M2)
x = gmres(A,b,restart,tol,maxit,M1,M2,x0)
[x,flag] = gmres(___)
[x,flag,relres] = gmres(___)
[x,flag,relres,iter] = gmres(___)
[x,flag,relres,iter,resvec] = gmres(___)

说明

        x = gmres(A,b) 尝试使用广义最小残差法求解关于 x 的线性系统 A*x = b。如果尝试成功，gmres 会显示一条消息来确认收敛。如果 gmres 无法在达到最大迭代次数后收敛或出于任何原因暂停，则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。对于此语法，gmres 不会重新启动；最大迭代次数为 min(size(A,1),10)。

        x = gmres(A,b,restart) 每 restart 次内迭代重新启动该方法一次。最大外迭代次数为 outer = min(size(A,1)/restart,10)。最大总迭代次数为 restart*outer，因为 gmres 对每次外迭代执行 restart 次内迭代。如果 restart 为 size(A,1) 或 []，则 gmres 不重新启动，最大总迭代次数为 min(size(A,1),10)。

        x = gmres(A,b,restart,tol) 指定该方法的容差。默认容差是 1e-6。

        x = gmres(A,b,restart,tol,maxit) 指定最大外迭代次数，使总迭代次数不超过 restart*maxit。如果 maxit 为 []，gmres 使用默认值 min(size(A,1)/restart,10)。如果 restart 为 size(A,1) 或 []，则最大总迭代次数为 maxit（而不是 restart*maxit）。如果 gmres 未能在最大迭代次数内收敛，将显示诊断消息。

        x = gmres(A,b,restart,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解方程组 M^−1Ax=M^−1b 来计算 x。使用预条件子矩阵可以改善问题的数值属性和计算的效率。

        x = gmres(A,b,restart,tol,maxit,M1,M2) 指定预条件子矩阵 M 的因子，使得 M = M1*M2。

        x = gmres(A,b,restart,tol,maxit,M1,M2,x0) 指定解向量 x 的初始估计值。默认值为由零组成的向量。

        [x,flag] = gmres(___) 返回一个标志，指示算法是否成功收敛。当 flag = 0 时，收敛成功。您可以将此输出语法用于之前的任何输入参数组合。如果指定了 flag 输出，gmres 将不会显示任何诊断消息。

        [x,flag,relres] = gmres(___) 还返回相对残差 norm(M\(b-A*x))/norm(M\b)，其中包括预条件子矩阵 M。如果 flag 为 0，则 relres <= tol。

        [x,flag,relres,iter] = gmres(___) 还以向量 [outer inner] 形式返回计算出 x 时的内迭代和外迭代次数。外迭代数在 0 <= iter(1) <= maxit 范围内，内迭代数在 0 <= iter(2) <= restart 范围内。

        [x,flag,relres,iter,resvec] = gmres(___) 还会返回每次内迭代中的残差范数向量（包括第一个残差 norm(M\(b-A*x0))）。这些是预条件方程组的残差范数。

示例

线性系统的迭代解

        使用采用默认设置的 gmres 求解系数矩阵为方阵的线性系统，然后在求解过程中调整使用的容差和迭代次数。

        创建一个随机稀疏矩阵 A，其密度为 50%，主对角线上的元素不为零。另为 Ax=b 的右侧创建随机向量 b。

rng default
A = sprandn(400,400,0.5) + 12*speye(400);
b = rand(400,1);

        使用 gmres 求解 Ax=b。输出显示包括相对残差 的值。

x = gmres(A,b);

gmres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.35.

        默认情况下，gmres 使用 10 次迭代和容差 1e-6，对于此矩阵，算法无法在 10 次迭代后收敛。由于残差仍然很大，这说明需要更多的迭代（或预条件子矩阵）。也可以使用更大的容差，使算法更容易收敛。

        使用容差 1e-4 和 100 次迭代再次求解方程组。

tol = 1e-4;
maxit = 100;

x = gmres(A,b,[],tol,maxit);
gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.0045.

        即使采用更宽松的容差和更多迭代，残差的改进也不足以收敛。当迭代算法以这种方式停滞时，显然需要预条件子矩阵。不过，gmres 也有控制内迭代次数的输入。通过为内迭代指定值，gmres 将针对每次外迭代做更多的工作。

        使用 restart 值 100 和 maxit 值 20 再次求解方程组。gmres 不是一次执行 100 次迭代，而是在两次重启之间执行 100 次迭代，并重复 20 次。

restart = 100;
maxit = 20;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 2 (inner iteration 75) to a solution with relative residual 9.3e-05.

        在这种情况下，为 gmres 指定较大的重启值会使其能够在允许的迭代次数内收敛于某个解。但是，当 A 也很大时，大的重启值会消耗大量内存。

使用 gmres 预条件子而不重启

        检查使用预条件子矩阵和非重启 gmres 求解线性系统的效果。加载 west0479，它是一个非对称的 479×479 实稀疏矩阵。

load west0479
A = west0479;

        定义 b 以使 Ax=b 的实际解是全为 1 的向量。

b = sum(A,2);

设置容差和最大迭代次数。

tol = 1e-12;
maxit = 20;

使用 gmres 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息：

• x 是计算 A*x = b 所得的解。

• fl0 是指示算法是否收敛的标志。

• rr0 是计算的解 x 的相对残差。

• it0 是二元素向量 [outer inner]，它表示计算 x 时的内迭代和外迭代次数。

• rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量。

[x,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.7603
it0
it0 = 1×2

     1    20

        gmres 未在请求的 20 次迭代内收敛至请求的容差 1e-12，因此 fl0 为 1。gmres 返回的最佳近似解是最后一个（如 it0(2) = 20 所示）。MATLAB® 将残差历史记录存储于 rv0 中。

        为了有助于缓慢收敛，可以指定预条件子矩阵。由于 A 是非对称的，请使用 ilu 生成预条件子 M=L U。指定调降容差，以忽略值小于 1e-6 的非对角线元。通过指定 L 和 U 作为 gmres 的输入，求解预条件方程组 M^−1A x=M^−1b。请注意，当指定预条件子时，gmres 会针对输出 rr1 和 rv1 计算预条件方程组的残差范数。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1
fl1 = 0
rr1
rr1 = 6.8998e-14
it1
it1 = 1×2

     1     6

        在第六次迭代中，使用 ilu 预条件子产生的相对残差小于规定的容差 1e-12。第一个残差 rv1(1) 是 norm(U\(L\b))，其中 M = L*U。最后一个残差 rv1(end) 是 norm(U\(L\(b-A*x1)))。

        可以通过绘制每次迭代的相对残差来跟踪 gmres 的进度。绘制每个解的残差历史记录图，并添加一条表示指定容差的线。

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

如图所示：

使用 gmres 预条件子并重启

        使用指定了预条件子的重启的 gmres。

        加载 west0479，它是一个非对称的 479×479 实稀疏矩阵。

load west0479
A = west0479;

        定义 b 以使 Ax=b 的实际解是全为 1 的向量。

b = sum(A,2);

        构造一个调降容差为 1e-6 的不完全 LU 预条件子。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

        使用重新启动的 gmres 的好处是限制执行该方法所需的内存量。如果不重新启动，gmres 需要存储 maxit 向量来保存基本 Krylov 子空间。此外，gmres 必须在每一步与以前的所有向量正交。重新启动可限制使用的工作区量以及每次外迭代执行的工作量。

        使用不完全 LU 因子作为预条件子执行 gmres(3)、gmres(4) 和 gmres(5)。使用容差 1e-12 和最多 20 次外迭代。

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3
fl3 = 0
fl4
fl4 = 0
fl5
fl5 = 0

        fl3、fl4 和 fl5 均为 0，因为在每种情况下重新启动的 gmres 使相对残差趋向小于 1e-12 的规定公差。

        以下绘图显示了每个重启的 gmres 方法的收敛历史记录。gmres(3) 收敛于外迭代 5、内迭代 3 (it3 = [5, 3])，这将与外迭代 6、内迭代 0 相同，因此最终刻度线上的标记是 6。

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

如图所示：

        一般情况下，内迭代次数越多，gmres 针对每个外迭代所做的工作会越多，收敛速度也会越快。

提供初始估计值

        检查向 gmres 提供解的初始估计值的效果。

        创建一个三对角稀疏矩阵。使用每行的总和作为 Ax=b 右侧的向量，使 x 的预期解是由 1 组成的向量。

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

        使用 gmres 求解 Ax=b 两次：一次是使用默认的初始估计值，一次是使用解的良好初始估计值。对两次求解均使用 200 次迭代和默认容差。将第二种求解中的初始估计值指定为所有元素都等于 0.99 的向量。

maxit = 200;
x1 = gmres(A,b,[],[],maxit);
gmres converged at iteration 27 to a solution with relative residual 9.5e-07.
x0 = 0.99*e;
x2 = gmres(A,b,[],[],maxit,[],[],x0);

gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

        在这种情况下，提供初始估计值可以使 gmres 更快地收敛。

返回中间结果

        还可以通过在 for 循环中调用 gmres 来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代，并存储计算出的解。然后，将该解用作下一批迭代的初始向量。

        例如，以下代码会循环执行四次，每次执行 100 次迭代，并在 for 循环中每通过一次后均存储解向量：

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = gmres(A,b,[],tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

        X(:,k) 是在 for 循环的第 k 次迭代时计算的解向量，R(k) 是该解的相对残差。

使用函数句柄代替数值矩阵

        通过为 gmres 提供用来计算 A*x 的函数句柄（而非系数矩阵 A）来求解线性系统。

        gallery 生成的 Wilkinson 测试矩阵之一是 21×21 三对角矩阵。预览该矩阵。

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

        Wilkinson 矩阵有特殊的结构，因此您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。此外，只有主对角线具有不等于 1 的非零值。

        表达式 Ax 变为：

结果向量可以写为三个向量的和：

        在 MATLAB® 中，编写一个函数来创建这些向量并将它们相加，从而给出 A*x 的值：

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

（该函数作为局部函数保存在示例的末尾。）

        现在，通过为 gmres 提供用于计算 A*x 的函数句柄，求解线性系统 Ax=b。重启前，使用容差 1e-12、15 次外迭代和 10 次内迭代。

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)
gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

        检查 afun(x1) 是否产生由 1 组成的向量。

afun(x1)
ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

局部函数

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

提示

• ​大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)。当 A 是方阵时，可以使用 equilibrate 来改进其条件数，它本身就能使大多数迭代求解器更容易收敛。但如果您随后会对经平衡处理的矩阵 B = R*P*A*C 进行因式分解，使用 equilibrate 还可以获得质量更好的预条件子矩阵。​

• 可以使用矩阵重新排序函数（如 dissect 和 symrcm）来置换系数矩阵的行和列，并在系数矩阵被分解以生成预条件子时最小化非零值的数量。这可以减少后续求解预条件线性系统所需的内存和时间。

参考

        [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 and Martin H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput., July 1986, Vol. 7, No. 3, pp. 856-869.