MATLAB中qmr函数用法

最新推荐文章于 2024-11-05 10:03:31 发布
最新推荐文章于 2024-11-05 10:03:31 发布
阅读量843 收藏 17
点赞数 13
分类专栏： Matlab 文章标签： matlab 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jk_101/article/details/133174002
版权

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 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,M) 指定预条件子矩阵 M 并通过有效求解方程组 M^−1Ax=M^−1
b 来计算 x。使用预条件子矩阵可以改善问题的数值属性和计算的效率。

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.

jk_101
关注
专栏目录
matlab数学手册
01-11
- 使用`numel`函数计算矩阵中元素的数量。 - 示例：`X = numel(A);` **1.3 矩阵分解** - **1.3.1 Cholesky分解** - 使用`chol`函数进行Cholesky分解。 - 示例：`Y = chol(A);` - **1.3.2 LU分解** - 使用`lu...
Matlab从入门到精通（五）--常用函数
热门推荐
数据之道
11-23 2万+
1、傅立叶(Fourier)变换及拉普拉斯(Laplace)变换 在MATLAB中，进行傅立叶变换的函数是： fourier(f,x,t)：求函数f(x)的傅立叶像函数F(t)。 ifourier(F,t,x)：求傅立叶像函数F(t)的原函数f(x)。   在MATLAB中，进行拉普拉斯变换的函数是： laplace(fx,x,t)：求函数f(x)的拉普拉斯像函数F(t)。 ilap...
matlab常用函数
good_learne的博客
09-07 1万+
一、软件操作函数 1）命令窗口函数： clc：清空命令窗口，使用向上箭头翻看命令。 open:打开文件，文本文件（*.doc），可执行文件（*.exe），图形文件（*.fig），超文本文件（*.html，*.htm），MATLAB数据库文件（*.mat），simulink模型文件*.mdl），MATLAB p文件（*.p），PDF文件（*,pdf），PPT文件（*.ppt），工程文件（*.prj...
算法之QMR--ST算法
Sizaif's 小屋
09-19 1017
RMQ (RangeMinimum/Maximum Query)问题是指：对于长度为n的数列A，回答若干询问RMQ(A,i,j)(i,j，返回数列A中下标在i,j里的最小(大）值，也就是说，RMQ问题是指求区间最值的问题。 解决的方法有:   1、朴素（即搜索），O(n)-O(qn) online。 2、线段树，O(n)-O(qlogn)online。 3、ST（实质是动态规划
MATLAB函数速查手册
u010510134的专栏
03-19 1万+
MATLAB函数速查手册
Matlab数学手册
07-07
- 使用`qmr`函数。 **1.5 特征值与二次型** - **1.5.1 特征值与特征向量的求法** - 使用`eig`函数求解。 - **1.5.2 提高特征值的计算精度** - 使用`balance`函数对矩阵进行平衡。 - **1.5.3 复对角矩阵转化为实...
MATLAB数学手册（pdf）
04-08
- 示例：使用`intersect`函数找出两个矩阵中的公共元素。 - **1.2.4 除法运算** - 左除`\`和右除`/`，分别用于求解线性方程组和比例问题。 - 示例：求解线性方程组`A \ b`；比例问题`A / b`。 - **1.2.5 矩阵...
【风力发电】基于Matlab的双馈风力发电机建模仿真设计
2401_84610415的博客
11-03 437
本文针对双馈风力发电机（DFIG）系统的建模与仿真，提出了一种基于MATLAB/Simulink的设计方法。双馈风力发电机在风能利用和电力输出方面具有较高的灵活性。通过MATLAB建模和仿真，验证了DFIG在不同风速条件下的性能表现。
基于数字图像的水果边缘检测方法研究
2401_84149290的博客
11-01 1267
基于数字图像的水果边缘检测方法研究1绪论1.1课题研究的背景我国水果品种繁多，资源丰富。作为一项新兴产业,在农业和农村经济发展中的地位已经十分重要,在很多地区己成为农村经济的支柱产业。纵观 世界水果园科技的发展, 二十一世纪我国的水果园也将从传统作业转化到以现代科学技术为基础的现代“集约持续水果园”上来。水果园的发展将全面地体现生产、经济、技术、社会和生态的可持续性。在水果分级方面,我国相继引进和开发了按大小和重量分级的分级 设备,近几年国内也出现了基于机器视觉技术的水果分级设备。但是,进口设备价格昂 贵,
基于 MATLAB的混沌序列图像加密算法的研究
2401_84462152的博客
11-02 579
图像信息生动形象，它已成为人类表达信息的重要手段之一，网络上的图像数据很多是要求发送方和接受都要进行加密通信，信息的安全与保密显得尤为重要，因此我想运用异或运算将数据进行隐藏，连续使用同一数据对图像数据两次异或运算图像的数据不发生改变，利用这一特性对图像信息进行加密保护。其内容，实现了图像数据的隐藏，这种混沌序列方法对图像数据的加密是十分有用，实现了信息的保密，在网络传输中能够很好地保护图像数据不被第三方轻易获取其内容，实现数据的隐藏保护。imwrite(e,' 加密后的 ','bmp');
零基础快速入门MATLAB
mailangduoduo的博客
11-02 985
在学习基于matlab仿真的数字图像处理教程之前，需要我们对matlab的基本语法有一定的了解，本篇博客带你快速入门matlab编程语言。本篇博客介绍了matlab的基本语法，适用于没学过该编程语言，快速入门的基于matlab仿真的数字信号处理实验，本篇博客在介绍时，并未对各种函数的各种参数展开讲解，需要时可以查看官方文档，更适合，更高效。希望本篇快速的matlab能够对你有所帮助！！！！！
MATLAB计算朗格朗日函数
corn1949的博客
10-31 551
MATLAB计算朗格朗日函数
034_Structural_Transient_In_Matlab结构动力学问题求解
withstand的博客
11-02 1083
函数进行可视化，这个函数可以把求解的结果映射到网络结构上，对点、线、面进行着色。当时求解这个问题，在最外面的竖直切面加载了一个静态的固定的力。注意，我们这里设置边界条件的函数，由小写字母开头，而设置好的边界条件的变量，首字母是大写的。这里，因为是一个动态问题，我们实际更想要更好地理解在动态载荷下，结构是如何发生变形和位移的。，这个函数提供了更多的可视化功能，可以对结构的变形、位移、应力、应变等进行可视化。加载的设置方法，类似于边界条件的设置。当这里使用一个函数是， 对函数的输入和输出有固定的要求。
【MATLAB源码-第206期】基于matlab的差分进化算法(DE)机器人栅格路径规划，输出做短路径图和适应度曲线。
最新发布
weixin_67624305的博客
11-05 420
差分进化算法的性能在很大程度上依赖于参数的选择，主要包括种群大小NP、缩放因子F和交叉概率CR。对于每个目标向量和对应的试验向量，如果试验向量的适应度（通常是目标函数值）优于目标向量的适应度，那么在下一代种群中替换目标向量；其中，𝑟𝑎𝑛𝑑𝑗(0,1)randj​(0,1) 是一个随机数，𝑟𝑎𝑛𝑑(𝑖)rand(i) 是一个随机选择的维度索引，保证每个试验向量至少有一个来自变异向量的分量。其中，𝑟1,𝑟2,𝑟3r1,r2,r3 是从种群中随机选择的三个不同的索引，𝐹F 是缩放因子，用于控制差分扰动的幅度。
计算机的错误计算（一百四十）
zaim1的博客
11-01 381
探讨MATLAB 中函数log1p(x) = ln(1+x)的计算精度。
【MATLAB源码-第200期】基于matlab的鸡群优化算法(CSO)机器人栅格路径规划，输出做短路径图和适应度曲线。
weixin_67624305的博客
10-31 1170
CSO算法首先将待优化问题中的解候选视为鸡群中的个体，根据其适应度（即解的质量），将这些个体分为不同的社会角色：公鸡、母鸡、小鸡和观察鸡。母鸡的位置更新则是基于它们与公鸡的相对位置，加上随机扰动，模拟母鸡跟随公鸡并在公鸡确定的区域内进行搜索的行为。观察鸡（Observer）：其余的个体，自由搜索并观察其它鸡的行为，以决定未来的行动策略。小鸡（Chick）：适应度较低的个体，它们依赖母鸡获得保护和食物，学习母鸡的行为。小鸡的位置更新是依赖于母鸡的位置，这模仿了小鸡跟随母鸡的行为。评估：计算每只鸡的适应度。
MATLAB车道检测与跟踪
2401_84457897的博客
10-29 617
它使用的本体匹 配算法利用的则是频率特征 ,通过对似然概率的估 计 ,得到最符合条件的车道线.。车道检测是智能车辆发展的智能因素。近年来,自主驾驶技术受到研究者们的普遍关注,车道线的提取是自主驾驶系统的重要组成部分.车道线的提取就是从车载摄像头获得的视频图像中,根据车道线的颜色,形状和纹理等特征,将车道线与背景进行分离,从而获得车道线的走向,车辆相对于车道线的位置等信息.从国内外研究现状来 车辆检测的目标是车道形状，包括宽度、曲率等几何参数，车辆在车道中的位置，包括横向偏移量，车辆与道路的夹角。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值