MATLAB中qmr函数用法

大家好，欢迎来到IT知识分享网。

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 qmr

提供初始估计值

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

        qmr函数的功能是求解线性系统 – 拟最小残差法。

语法

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

说明

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

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

x = qmr(A,b,tol,maxit) 指定要使用的最大迭代次数。如果 qmr 无法在 maxit 次迭代内收敛，将显示诊断消息。

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

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

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

[x,flag,relres] = qmr(___) 还会返回相对残差 norm(b-A*x)/norm(b)。如果 flag 为 0，则 relres <= tol。

[x,flag,relres,iter] = qmr(___) 还会返回计算出 x 时的迭代次数 iter。

[x,flag,relres,iter,resvec] = qmr(___) 还会在每次迭代中返回残差范数向量（包括第一个残差 norm(b-A*x0)）。

示例

线性系统的迭代解

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

        创建密度为 50% 的随机稀疏矩阵 A。另为 Ax=b 的右侧创建随机向量 b。

rng default A = sprand(400,400,.5); A = A'*A; b = rand(400,1);

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

x = qmr(A,b); qmr stopped at iteration 20 without converging to the desired tolerance 1e-06 because the maximum number of iterations was reached. The iterate returned (number 20) has relative residual 0.12.

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

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

x = qmr(A,b,1e-4,100); qmr 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.061.

        即使采用更宽松的容差和更多迭代，残差也并未改进多少。当迭代算法以这种方式停滞时，显然需要预条件子矩阵。

        计算 A 的不完全 Cholesky 分解，并使用 L’ 因子作为 qmr 的预条件子输入。

L = ichol(A); x = qmr(A,b,1e-4,100,L'); qmr converged at iteration 58 to a solution with relative residual 8.4e-05.

        使用预条件子可以充分改进问题的数值属性，使 qmr 能够收敛。

使用指定了预条件子的 qmr

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

load west0479 A = west0479;

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

b = sum(A,2);

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

tol = 1e-12; maxit = 20;

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

• x 是计算 A*x = b 所得的解。
• fl0 是指示算法是否收敛的标志。
• rr0 是计算的解 x 的相对残差。
• it0 是计算 x 时所用的迭代次数。
• rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量。

[x,fl0,rr0,it0,rv0] = qmr(A,b,tol,maxit); fl0 fl0 = 1 rr0 rr0 = 0.7984 it0 it0 = 17

        qmr 未在请求的 20 次迭代内收敛至请求的容差 1e-12，因此 fl0 为 1。第 17 次迭代是最佳近似解，也是按 it0 = 17 指示返回的近似解。

        为了有助于缓慢收敛，您可以指定预条件子矩阵。由于 A 是非对称的，请使用 ilu 生成预条件子 M=L U。指定调降容差，以忽略值小于 1e-6 的非对角线元。通过指定 L 和 U 作为 qmr 的输入，求解预条件方程组 

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

        在第六次迭代中，使用 ilu 预条件子产生的相对残差小于规定的容差 1e-12。输出 rv1(1) 为 norm(b)，输出 rv1(end) 为 norm(b-A*x1)。

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

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

如图所示：

提供初始估计值

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

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

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

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

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

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

返回中间结果

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

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

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

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

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

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

        创建一个非对称三对角矩阵。预览该矩阵。

A = gallery('wilk',21) + diag(ones(20,1),1) A = 21×21 10 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 8 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 2 0 0 0 0 0 0 0 0 0 0 ⋮ 

        由于此三对角矩阵有特殊的结构，可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。

表达式 A x 变为：

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

同样，A^T x 的表达式变为：

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

function y = afun(x,flag) if strcmp(flag,'notransp') % Compute A*x y = [0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + 2*[x(2:end); 0]; elseif strcmp(flag,'transp') % Compute A'*x y = 2*[0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + [x(2:end); 0]; end end

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

        现在，通过为 qmr 提供用于计算 A*x 和 A’*x 的函数句柄，求解线性系统 Ax=b。使用容差 1e-6 和 25 次迭代。指定 b 为 A 的行总和，使得 x 的实际解是由 1 组成的向量。

b = full(sum(A,2)); tol = 1e-6; maxit = 25; x1 = qmr(@afun,b,tol,maxit) qmr converged at iteration 19 to a solution with relative residual 4.7e-07. x1 = 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,flag) if strcmp(flag,'notransp') % Compute A*x y = [0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + 2*[x(2:end); 0]; elseif strcmp(flag,'transp') % Compute A'*x y = 2*[0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + [x(2:end); 0]; end end

拟最小残差法

​        开发 QMR 算法是为了改进 BiCG。GMRES 对 Krylov 子空间使用正交基并计算最小残差解，而 QMR 使用双正交基，因此只计算拟最小残差解。

        QMR 通常比 BiCG 收敛得更平滑，它还使用前瞻性的方法，以保证几乎任何情况下都不出故障。QMR 的计算成本仅略高于 BiCG [1]。

​提示

• ​大多数迭代方法的收敛取决于系数矩阵的条件数 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] Freund, Roland W. and Nöel M. Nachtigal, “QMR: A quasi-minimal residual method for non-Hermitian linear systems,” SIAM Journal: Numer. Math. 60, 1991, pp. 315–339.