bicg
求解线性系统 - 双共轭梯度法
语法
说明
示例
线性系统的迭代解
使用采用默认设置的 bicg
求解系数矩阵为方阵的线性系统,然后在求解过程中调整使用的容差和迭代次数。
创建密度为 50% 的随机稀疏矩阵 A
。另为 的右侧创建随机向量 b
。
rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);
使用 bicg
求解 。输出显示包括相对残差 的值。
x = bicg(A,b);
bicg stopped at iteration 20 without converging to the desired tolerance 1e-06 because the maximum number of iterations was reached. The iterate returned (number 7) has relative residual 0.45.
默认情况下,bicg
使用 20 次迭代和容差 1e-6
,对于此矩阵,算法无法在 20 次迭代后收敛。由于残差仍然很大,这说明需要更多的迭代(或预条件子矩阵)。您也可以使用更大的容差,使算法更容易收敛。
使用容差 1e-4
和 100 次迭代再次求解方程组。
x = bicg(A,b,1e-4,100);
bicg stopped at iteration 100 without converging to the desired tolerance 0.0001 because the maximum number of iterations was reached. The iterate returned (number 7) has relative residual 0.45.
即使采用更宽松的容差和更多迭代,残差也并未改进多少。当迭代算法以这种方式停滞时,显然需要预条件子矩阵。
计算 A
的不完全乔列斯基分解,并使用 L'
因子作为 bicg
的预条件子输入。
L = ichol(A); x = bicg(A,b,1e-4,100,L');
bicg converged at iteration 55 to a solution with relative residual 9.9e-05.
使用预条件子可以充分改进问题的数值属性,使 bicg
能够收敛。
使用指定了预条件子的 bicg
检查使用指定了预条件子矩阵的 bicg
来求解线性系统的效果。
加载 west0479,它是一个非对称的 479×479 实稀疏矩阵。
load west0479
A = west0479;
定义 b
以使 的实际解是全为 1 的向量。
b = sum(A,2);
设置容差和最大迭代次数。
tol = 1e-12; maxit = 20;
使用 bicg
根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息:
x
是计算A*x = b
所得的解。fl0
是指示算法是否收敛的标志。rr0
是计算的解x
的相对残差。it0
是计算x
时所用的迭代序号。rv0
是 的残差历史记录组成的向量。
[x,fl0,rr0,it0,rv0] = bicg(A,b,tol,maxit); fl0
fl0 = 1
rr0
rr0 = 1
it0
it0 = 0
bicg
未在请求的 20 次迭代内收敛至请求的容差 1e-12
,因此 fl0
为 1。实际上,bicg
的行为太差,因此初始估计值 x0 = zeros(size(A,2),1)
是最佳解,并会返回最佳解(如 it0 = 0
所示)。
为了有助于缓慢收敛,您可以指定预条件子矩阵。由于 A
是非对称的,请使用 ilu
生成预条件子 。指定调降容差,以忽略值小于 1e-6
的非对角线元。通过指定 L
和 U
作为 bicg
的输入,求解预条件方程组 。
setup = struct('type','ilutp','droptol',1e-6); [L,U] = ilu(A,setup); [x1,fl1,rr1,it1,rv1] = bicg(A,b,tol,maxit,L,U); fl1
fl1 = 0
rr1
rr1 = 4.1410e-14
it1
it1 = 6
在第六次迭代中,使用 ilu
预条件子产生的相对残差小于规定的容差 1e-12
。输出 rv1(1)
为 norm(b)
,输出 rv1(end)
为 norm(b-A*x1)
。
您可以通过绘制每次迭代的相对残差来跟踪 bicg
的进度。绘制每个解的残差历史记录图,并添加一条表示指定容差的线。
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')
提供初始估计值
检查向 bicg
提供解的初始估计值的效果。
创建一个三对角稀疏矩阵。使用每行的总和作为 右侧的向量,使 的预期解是由 1 组成的向量。
n = 900; e = ones(n,1); A = spdiags([e 2*e e],-1:1,n,n); b = sum(A,2);
使用 bicg
求解 两次:一次是使用默认的初始估计值,一次是使用解的良好初始估计值。对两次求解均使用 200 次迭代和默认容差。将第二种求解中的初始估计值指定为所有元素都等于 0.99
的向量。
maxit = 200; x1 = bicg(A,b,[],maxit);
bicg converged at iteration 35 to a solution with relative residual 9.5e-07.
x0 = 0.99*e; x2 = bicg(A,b,[],maxit,[],[],x0);
bicg converged at iteration 7 to a solution with relative residual 8.7e-07.
在这种情况下,提供初始估计值可以使 bicg
更快地收敛。
返回中间结果
您还可以通过在 for 循环中调用 bicg
来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代,并存储计算出的解。然后,将该解用作下一批迭代的初始向量。
例如,以下代码会循环执行四次,每次执行 100 次迭代,并在 for 循环中每通过一次后均存储解向量:
x0 = zeros(size(A,2),1); tol = 1e-8; maxit = 100; for k = 1:4 [x,flag,relres] = bicg(A,b,tol,maxit,[],[],x0); X(:,k) = x; R(k) = relres; x0 = x; end
X(:,k)
是在 for 循环的第 k
次迭代时计算的解向量,R(k)
是该解的相对残差。
使用函数句柄代替数值矩阵
通过为 bicg
提供用来计算 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
的非零三对角元素。
表达式 变为:
.
结果向量可以写为三个向量的和:
=。
同样, 的表达式变为:
.
.
在 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
(该函数作为局部函数保存在示例的末尾。)
现在,通过为 bicg
提供用于计算 A*x
和 A'*x
的函数句柄,求解线性系统 。使用容差 1e-6
和 25 次迭代。指定 为 的行总和,使得 的实际解是由 1 组成的向量。
b = full(sum(A,2)); tol = 1e-6; maxit = 25; x1 = bicg(@afun,b,tol,maxit)
bicg converged at iteration 19 to a solution with relative residual 4.8e-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
输入参数
A
— 系数矩阵
矩阵 | 函数句柄
系数矩阵,指定为方阵或函数句柄。该矩阵是线性系统 A*x = b
中的系数矩阵。通常,A
是大型稀疏矩阵或函数句柄,它返回大型稀疏矩阵和列向量的乘积。
将 A
指定为函数句柄
您可以选择将系数矩阵指定为函数句柄而不是矩阵。函数句柄返回矩阵向量乘积,而不是构建整个系数矩阵,从而使计算更加高效。
要使用函数句柄,请使用函数签名 function y = afun(x,opt)
。参数化函数说明如何在必要时为函数 afun
提供附加参数。函数 afun
必须满足下列条件:
afun(x,'notransp')
返回乘积A*x
。afun(x,'transp')
返回乘积A'*x
。
可接受的函数的一个示例是:
function y = afun(x,opt,B,C,n) if strcmp(opt,'notransp') y = [B*x(n+1:end); C*x(1:n)]; else y = [C'*x(n+1:end); B'*x(1:n)]; end
afun
使用 B
和 C
中的值来计算 A*x
或 A'*x
(具体取决于指定的标志),而不用真正建立整个矩阵。
数据类型: double
| function_handle
复数支持: 是
b
— 线性方程的右侧
列向量
线性方程的右侧,指定为列向量。b
的长度必须等于 size(A,1)
。
数据类型: double
复数支持: 是
tol
— 方法容差
[]
或 1e-6
(默认) | 正标量
方法容差,指定为正标量。计算中使用此输入可在准确度和运行时间之间进行权衡。bicg
必须在允许的迭代次数内满足容差才能成功。较小的 tol
值意味着解必须更精确才能成功完成计算。
数据类型: double
maxit
— 最大迭代次数
[]
或 min(size(A,1),20)
(默认) | 正整数标量
最大迭代次数,指定为正整数标量。增加 maxit
的值,以允许 bicg
进行更多迭代,从而满足容差 tol
。通常,较小的 tol
值意味着需要更多迭代才能成功完成计算。
M
, M1
, M2
— 预条件子矩阵(以单独参量指定)
eye(size(A))
(默认) | 矩阵 | 函数句柄
预条件子矩阵,指定为由矩阵或函数句柄组成的单独参量。您可以指定预条件子矩阵 M
或其矩阵因子 M = M1*M2
来改进线性系统的数值方面,使 bicg
更容易快速收敛。对于系数方阵,您可以使用不完全矩阵分解函数 ilu
和 ichol
来生成预条件子矩阵。您还可以在分解之前使用 equilibrate
来改进系数矩阵的条件数。有关预条件子的详细信息,请参阅线性系统的迭代方法。
bicg
将未指定的预条件子视为单位矩阵。
将 M
指定为函数句柄
您可以选择将 M
、M1
或 M2
中的任一个指定为函数句柄而不是矩阵。函数句柄执行矩阵向量运算,而不是构建整个预条件子矩阵,从而使计算更加高效。
要使用函数句柄,首先创建签名为 function y = mfun(x,opt)
的函数。参数化函数说明如何在必要时为函数 mfun
提供附加参数。函数 mfun
必须满足下列条件:
mfun(x,'notransp')
返回M\x
或M2\(M1\x)
的值。mfun(x,'transp')
返回M'\x
或M1'\(M2'\x)
的值。
可接受的函数的一个示例是:
function y = mfun(x,opt,a,b) if strcmp(opt,'notransp') y = x.*a; else y = x.*b; end end
mfun
使用 a
和 b
来计算 M\x = x*a
或 M'\x = x*b
(具体取决于指定的标志),而不用真正建立整个矩阵 M
。
数据类型: double
| function_handle
复数支持: 是
x0
— 初始估计值
[]
或由零组成的列向量 (默认) | 列向量
初始估计值,指定为长度等于 size(A,2)
的列向量。如果您能为 bicg
提供比默认的零向量更合理的初始估计值 x0
,则它可以节省计算时间并帮助算法更快地收敛。
数据类型: double
复数支持: 是
输出参量
flag
— 收敛标志
标量
收敛标志,返回下表中的标量值之一。收敛标志指示计算是否成功,并区分几种不同形式的失败。
标志值 | 收敛 |
---|---|
| 成功 - |
| 失败 - |
| 失败 - 预条件子矩阵 |
| 失败 - |
| 失败 - 由 |
iter
— 迭代序号
标量
迭代序号,以标量形式返回。此输出指示计算出 x
的解时所用的迭代序号。
数据类型: double
详细信息
提示
大多数迭代方法的收敛取决于系数矩阵的条件数
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.
扩展功能
基于线程的环境
使用 MATLAB® backgroundPool
在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool
加快代码运行速度。
此函数完全支持基于线程的环境。有关详细信息,请参阅在基于线程的环境中运行 MATLAB 函数。
GPU 数组
通过使用 Parallel Computing Toolbox™ 在图形处理单元 (GPU) 上运行来加快代码执行。
用法说明和限制:
当输入
A
是稀疏矩阵时:仅支持一个稀疏矩阵预条件子
M
。如果使用两个预条件子
M1
和M2
,则它们都必须是函数。对于 GPU 数组,
bicg
不会检测停滞(标志 3)。此时,它会报告收敛失败(标志 1)。
有关详细信息,请参阅在 GPU 上运行 MATLAB 函数 (Parallel Computing Toolbox)。
分布式数组
使用 Parallel Computing Toolbox™ 在集群的组合内存中对大型数组进行分区。
版本历史记录
在 R2006a 之前推出
MATLAB 命令
您点击的链接对应于以下 MATLAB 命令:
请在 MATLAB 命令行窗口中直接输入以执行命令。Web 浏览器不支持 MATLAB 命令。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)