目录
minres函数的功能是求解线性系统 - 最小残差法。
语法
x = minres(A,b)
x = minres(A,b,tol)
x = minres(A,b,tol,maxit)
x = minres(A,b,tol,maxit,M)
x = minres(A,b,tol,maxit,M1,M2)
x = minres(A,b,tol,maxit,M1,M2,x0)
[x,flag] = minres(___)
[x,flag,relres] = minres(___)
[x,flag,relres,iter] = minres(___)
[x,flag,relres,iter,resvec] = minres(___)
[x,flag,relres,iter,resvec,resveccg] = minres(___)
说明
x = minres(A,b) 尝试使用最小残差法求解关于 x 的线性系统 A*x = b。如果尝试成功,minres 会显示一条消息来确认收敛。如果 minres 无法在达到最大迭代次数后收敛或出于任何原因暂停,则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。
x = minres(A,b,tol) 指定该方法的容差。默认容差是 1e-6。
x = minres(A,b,tol,maxit) 指定要使用的最大迭代次数。如果 minres 无法在 maxit 次迭代内收敛,将显示诊断消息。
x = minres(A,b,tol,maxit,M) 指定对称正定预条件子矩阵 M 并通过有效求解关于 y 的方程组 来计算 x,其中。该算法不显式形成 H。使用预条件子矩阵可以改善问题的数值属性和计算的效率。
x = minres(A,b,tol,maxit,M1,M2) 指定预条件子矩阵 M 的因子,使得 M = M1*M2。
x = minres(A,b,tol,maxit,M1,M2,x0) 指定解向量 x 的初始估计值。默认值为由零组成的向量。
[x,flag] = minres(___) 返回一个标志,指示算法是否成功收敛。当 flag = 0 时,收敛成功。您可以将此输出语法用于之前的任何输入参数组合。如果指定了 flag 输出,minres 将不会显示任何诊断消息。
[x,flag,relres] = minres(___) 还会返回计算的解中的残差。如果 flag 为 0,则 relres <= tol。
[x,flag,relres,iter] = minres(___) 还会返回计算出 x 时的迭代次数 iter。
[x,flag,relres,iter,resvec] = minres(___) 还会在每次迭代中返回残差范数向量(包括第一个残差 norm(b-A*x0))。
[x,flag,relres,iter,resvec,resveccg] = minres(___) 还会在每次迭代中返回共轭梯度残差范数向量。
示例
线性系统的迭代解
使用采用默认设置的 minres 求解系数矩阵为方阵的线性系统,然后在求解过程中调整使用的容差和迭代次数。
创建一个稀疏对称三对角矩阵 A 作为系数矩阵。使用 A 的行总和作为 Ax=b 右侧的向量 b,以便 x 的预期解是由 1 组成的向量。
n = 400;
on = ones(n,1);
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2);
使用 minres 求解 Ax=b。输出显示包括相对残差 ‖b−Ax‖/‖b‖ 的值。
x = minres(A,b);
minres 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.017.
默认情况下,minres 使用 20 次迭代和容差 1e-6,对于此矩阵,算法无法在 20 次迭代后收敛。由于残差的数量级为 1e-2,显然需要更多迭代。也可以使用更大的容差,使算法更容易收敛。
使用容差 1e-4 和 250 次迭代再次求解方程组。
x = minres(A,b,1e-4,250);
minres converged at iteration 200 to a solution with relative residual 7e-13.
使用指定了预条件子的 minres
检查使用指定了预条件子矩阵的 minres 来求解线性系统的效果。
创建一个对称正定带状系数矩阵。
A = delsq(numgrid('S',102));
定义 b 以使 Ax=b 的实际解是全为 1 的向量。
b = sum(A,2);
设置容差和最大迭代次数。
tol = 1e-12;
maxit = 100;
使用 minres 根据请求的容差和迭代次数求解。指定六个输出以返回有关求解过程的信息:
-
x 是计算 A*x = b 所得的解。
-
fl0 是指示算法是否收敛的标志。
-
rr0 是计算的解 x 的相对残差。
-
it0 是计算 x 时所用的迭代次数。
-
rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量。
-
rvcg0 是的共轭梯度残差历史记录组成的向量。
[x,fl0,rr0,it0,rv0,rvcg0] = minres(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.0013
it0
it0 = 100
minres 未在请求的 100 次迭代内收敛至请求的容差 1e-12,因此 fl0 为 1。
为了帮助收敛,可以指定预条件子矩阵。由于 A 是对称矩阵,请使用 ichol 生成预条件子 。指定 'ict' 选项以使用阈值调降不完全 Cholesky 分解,并指定对角线偏移量值 1e-6 以避免非正主元。通过指定 L 和 L' 作为 minres 的输入,求解预条件方程组。
setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = minres(A,b,tol,maxit,L,L');
fl1
fl1 = 0
rr1
rr1 = 2.6412e-15
it1
it1 = 4
在第四次迭代中,使用 ilu 预条件子产生的相对残差小于规定的容差 1e-12。输出 rv1(1) 为 norm(b),输出 rv1(end) 为 norm(b-A*x1)。
可以通过绘制每次迭代的相对残差来跟踪 minres 的进度。绘制每个解的残差历史记录图,并添加一条表示指定容差的线。
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','ICHOL preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')
如图所示:
提供初始估计值
检查向 minres 提供解的初始估计值的效果。
创建一个三对角稀疏矩阵。使用每行的总和作为 Ax=b 右侧的向量,使 x 的预期解是由 1 组成的向量。
n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);
使用 minres 求解 Ax=b 两次:一次是使用默认的初始估计值,一次是使用解的良好初始估计值。对两次求解均使用 200 次迭代和默认容差。将第二种求解中的初始估计值指定为所有元素都等于 0.99 的向量。
maxit = 200;
x1 = minres(A,b,[],maxit);
minres converged at iteration 27 to a solution with relative residual 9.5e-07.
x0 = 0.99*e;
x2 = minres(A,b,[],maxit,[],[],x0);
minres converged at iteration 7 to a solution with relative residual 6.7e-07.
在这种情况下,提供初始估计值可以使 minres 更快地收敛。
返回中间结果
还可以通过在 for 循环中调用 minres 来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代,并存储计算出的解。然后,将该解用作下一批迭代的初始向量。
例如,以下代码会循环执行四次,每次执行 100 次迭代,并在 for 循环中每通过一次后均存储解向量:
x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
[x,flag,relres] = minres(A,b,tol,maxit,[],[],x0);
X(:,k) = x;
R(k) = relres;
x0 = x;
end
X(:,k) 是在 for 循环的第 k 次迭代时计算的解向量,R(k) 是该解的相对残差。
使用函数句柄代替数值矩阵
通过为 minres 提供用来计算 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
(该函数作为局部函数保存在示例的末尾。)
现在,通过为 minres 提供用于计算 A*x 的函数句柄,求解线性系统 Ax=b。使用容差 1e-12 和 50 次迭代。
b = ones(21,1);
tol = 1e-12;
maxit = 50;
x1 = minres(@afun,b,tol,maxit)
minres converged at iteration 11 to a solution with relative residual 4.1e-16.
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
最小残差法
MINRES 和 SYMMLQ 方法是支撑共轭梯度法 PCG 的 Lanczos 方法的变体。像 PCG 一样,系数矩阵仍需是对称矩阵,但 MINRES 和 SYMMLQ 允许它是不定矩阵(不要求所有特征值都必须为正值)。这是通过避免 Lanczos 方法中通常存在的隐式 LU 分解来实现的,当不定矩阵遇到主元为零的情况时,该方法容易出现故障。
MINRES 最小化 2-范数残差,而 SYMMLQ 使用 LQ 分解求解投影方程组,并使残差与所有先前的残差正交。GMRES 方法旨在将 MINRES 推广到非对称问题 [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] Paige, C. C. and M. A. Saunders, “Solution of Sparse Indefinite Systems of Linear Equations.” SIAM J. Numer. Anal., Vol.12, 1975, pp. 617-629.