MATLAB中lsqr函数用法

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 lsqr

提供初始估计值

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


        lsqr函数的功能是求解线性系统 - 最小二乘法。

语法

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

说明

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

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

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

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

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

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

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

        [x,flag,relres] = lsqr(___) 还返回计算的解 x 的残差。如果 flag 是 0,则 x 是最小化 norm(b-A*x) 的最小二乘解。如果 relres 很小,则 x 也是相容的解,因为 relres 表示 norm(b-A*x)/norm(b)。

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

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

        [x,flag,relres,iter,resvec,lsvec] = lsqr(___) 还会返回 lsvec,即每次迭代的缩放标准方程误差的估计值。

示例

线性系统的迭代解

        使用采用默认设置的 lsqr 求解矩形线性系统,然后调整求解过程中使用的容差和迭代次数。

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

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

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

x = lsqr(A,b);
lsqr 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.26.

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

        使用容差 1e-4 和 70 次迭代再次求解方程组。指定六个输出,以返回计算的解的相对残差 relres,以及残差历史记录 resvec 和最小二乘残差历史记录 lsvec。

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,1e-4,70);
flag
flag = 0

        由于 flag 为 0,因此该算法能够在指定的迭代次数内满足所需的误差容限。通常可以一起调整容差和迭代次数,以这种方式在速度和精度之间进行权衡。

        检查计算的解的相对残差和最小二乘残差。

relres
relres = 0.2625
lsres = lsvec(end)
lsres = 2.5375e-04

        这些残差范数表明 x 是最小二乘解,因为 relres 不小于 1e-4 的指定容差。由于线性系统不存在相容的解,求解器的最佳做法是使最小二乘残差满足容差。

        绘制残差历史记录。相对残差 resvec 很快达到最小值且无法取得进一步进展,而最小二乘残差 lsvec 在后续迭代中继续最小化。

N = length(resvec);
semilogy(0:N-1,lsvec,'--o',0:N-1,resvec,'-o')
legend("Least-squares residual","Relative residual")

如图所示:

使用指定了预条件子的 lsqr

        检查使用指定了预条件子矩阵的 lsqr 来求解线性系统的效果。

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

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

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

  • it 是计算 x 时所用的迭代次数。

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

  • lsrv 是最小二乘残差历史记录组成的向量。

[x,fl,rr,it,rv,lsrv] = lsqr(A,b,tol,maxit);
fl
fl = 1
rr
rr = 0.0017
it
it = 20

由于 fl = 1,算法未在最大迭代次数内收敛于指定的容差。

        为了有助于缓慢收敛,可以指定预条件子矩阵。由于 A 是非对称的,请使用 ilu 生成分解形式的预条件子 M=L U。指定调降容差,以忽略值小于 1e-6 的非对角线元。通过将 L 和 U 指定为 lsqr 的 M1 和 M2 输入,求解 y=Mx 的预条件方程组 AM^−1(M x)=b。

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

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

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

semilogy(0:length(rv)-1,rv/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')

如图所示:

提供初始估计值

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

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

A = sprand(700,900,0.1);
b = sum(A,2);

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

maxit = 75;
x1 = lsqr(A,b,[],maxit);
lsqr converged at iteration 64 to a solution with relative residual 8.7e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = lsqr(A,b,[],maxit,[],[],x0);
lsqr converged at iteration 26 to a solution with relative residual 9.6e-07.

随着初始估计值接近预期解,lsqr 能够在更少的迭代次数后收敛。

返回中间结果

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

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

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

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

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

        通过为 lsqr 提供用来计算 A*x 和 A'*x 的函数句柄(而非系数矩阵 A)来求解线性系统。

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

A = gallery('wilk',21) + diag(ones(20,1),1)
A = 21×2110     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     01     9     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     00     1     8     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     00     0     1     7     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     00     0     0     1     6     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     00     0     0     0     1     5     2     0     0     0     0     0     0     0     0     0     0     0     0     0     00     0     0     0     0     1     4     2     0     0     0     0     0     0     0     0     0     0     0     0     00     0     0     0     0     0     1     3     2     0     0     0     0     0     0     0     0     0     0     0     00     0     0     0     0     0     0     1     2     2     0     0     0     0     0     0     0     0     0     0     00     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*xy = [0; x(1:20)] ...+ [(10:-1:0)'; (1:10)'].*x ...+ 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*xy = 2*[0; x(1:20)] ...+ [(10:-1:0)'; (1:10)'].*x ...+ [x(2:end); 0];
end
end

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

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

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = lsqr(@afun,b,tol,maxit)
lsqr converged at iteration 21 to a solution with relative residual 4.5e-12.
x1 = 21×11.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000⋮

局部函数

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

最小二乘法

        ​最小二乘 (LSQR) 算法是对矩形矩阵的共轭梯度 (CG) 方法的改进。分析表明,A*x = b 的 LSQR 产生的残差与标准方程 A'*A*x = A'*b 的 CG 相同,但 LSQR 具有更好的数值属性,因此通常更可靠 [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, "LSQR: An Algorithm for Sparse Linear Equations And Sparse Least Squares," ACM Trans. Math. Soft., Vol.8, 1982, pp. 43-71.

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/400848.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【Qt】QWidget的windowIcon属性

QWidget的windowIcon属性 windowIcon表示窗口的图标 当我们使用默认的windowIcon的时候,其窗口的图标如下: API说明 windowIcon() 获取到控件的窗⼝图标. 返回 QIcon 对象. setWindowIcon(const QIcon& icon) 设置控件的窗⼝图标. 在Qt中&…

【STM32 FreeRTOS】任务

使用 RTOS 的实时应用程序可以被构建为一组独立的任务。每个任务在自己的上下文中执行,不依赖于系统内的其他任务或 RTOS 调度器本身。在任何时间点,应用程序中只能执行一个任务,实时 RTOS 调度器负责决定所要执行的任务。因此, R…

react中使用nextjs框架,前端调后端接口跨域解决方式

前端在项目目录中next.config.js文件中添加以下代码 async rewrites() {return [{source: "/api/:path*",destination: ${process.env.NEXT_PUBLIC_API_DOMAIN}/api/:path*,basePath: false}]} 截图: source: "/api/:path*": 定义了一个 URL …

C++ 面试题常用总结 详解(满足c++ 岗位必备,不定时更新)

📚 本文主要总结了一些常见的C面试题,主要涉及到语法基础、STL标准库、内存相关、类相关和其他辅助技能,掌握这些内容,基本上就满足C的岗位技能(红色标记为重点内容),欢迎大家前来学习指正&…

【MySQL 03】库的操作 (带思维导图)

文章目录 🌈 一、创建数据库🌈 二、查看数据库🌈 三、使用数据库🌈 四、修改数据库🌈 五、删除数据库🌈 六、备份数据库🌈 七、恢复数据库🌈 八、字符集和校验规则⭐ 1. 查看系统默认…

HAProxy 效能飞跃先锋队

目录 一 负载均衡 1.1 四层负载 1.2 七层负载 1.3 四层负载和七层负载的区别 二 Haproxy简介 2.1 概念和内容 2.2 haproxy的基本配置信息 2.2.1 global 配置 2.2.2 proxies 配置 三 Haproxy的算法 3.1 静态算法 3.2 动态算法 3.3 其他算法 四 高级功能及配置 4.…

Mysql原理与调优-索引原理及使用

目录 1.绪论 2.索引原理 2.1 索引采用的数据结构 2.1.1 B树 1.什么是B树 2.B树的优缺点 2.1.2 B树 1.什么是B树 3.B树的优缺点 2.2.3 hash 2.2.4 总结 2.2 InnoDB数据存储结构(聚簇索引) 2.2.1 自底向上探寻索引 3.2.2 查询一条数据的完整流程 3.2.3 聚簇索引 2…

奥运科技观察:AI PC,如何成为当代体育精神的数字捍卫者?

作者 | 曾响铃 文 | 响铃说 数字孪生帮助体育馆建设、超高清直播……这届奥运会科技感拉满,几乎所有前沿技术都能在奥运的赛事运营中发现。 而AI大时代,AI如何帮助帮助奥运会顺利举办、如何帮助运动员拥有更好的表现,同样值得业界关注&…

haproxy最强攻略

1、负载均衡 负载均衡(Load Balance,简称 LB)是高并发、高可用系统必不可少的关键组件,目标是 尽力将网络流量平均分发到多个服务器上,以提高系统整体的响应速度和可用性。 负载均衡的主要作用如下: 高并发…

# Spring Cloud Alibaba Nacos_配置中心与服务发现(四)

Spring Cloud Alibaba Nacos_配置中心与服务发现(四) 一、Nacos 配置管理-集群部署 1、 把 nacos 应用程序包,复制3份,分别命名为 nacos1, nacos2, nacos3 分别在 conf 目录下,修改 application.properties 配置文件…

数据结构——循环队列

目录 循环队列的基本知识 循环队列的实现 定义 各个接口的实现 循环队列的基本知识 循环队列的定义 循环队列(Circular Queue)是一种使用固定大小的数组实现的队列,它将数组的首尾相连,形成环形,以充分利用空间并实…

Spring Boot的配置文件

目录 一、配置文件 1.properties为后缀的配置文件 1.1基本语法 1.2读取配置文件 1.3properties的优缺点 1.4加中文注释出现乱码 2.yml格式的配置文件 2.1基础语法 2.2读取配置文件 2.2.1对象存储到配置文件中 2.3yml的优缺点 2.4用不用加单引号或者双引号呢&#xf…

【C语言篇】编译和链接以及预处理介绍(上篇)

文章目录 前言翻译环境和运行环境翻译环境编译预处理(预编译)编译词法分析语法分析语义分析 汇编 链接 运行环境预处理(预编译)详解预定义符号#define定义常量#define定义宏带有副作用的宏参数宏替换的规则宏和函数的对比 写在最后…

opencv基础的图像操作

1.读取图像,显示图像,保存图像 #图像读取、显示与保存 import numpy as np import cv2 imgcv2.imread(./src/1.jpg) #读取 cv2.imshow("img",img) #显示 cv2.imwrite("./src/2.jpg",img) #保存 cv2.waitKey(0) #让程序进入主循环(让…

RAG系列之四:深入浅出 Embedding

在 RAG 系列之三:文本切分中介绍了如何将文本切分成更小的语义单元,接下来便是将拆分的文本块进行向量化。 什么是文本向量化? 文本向量化就是将文本数据转成数字数据,例如:将文本 It was the best of times, it was…

Android全面解析之context机制(二): 从源码角度分析context创建流程(上)

前言 这篇文章从源码角度分析context创建流程。 在上一篇Android全面解析之Context机制(一) :初识context一文中讲解了context的相关实现类。经过前面的讨论,读者对于context在心中有了一定的理解。但始终觉得少点什么:activity是什么时候被创建的&…

Python数据可视化案例——地图

目录 简单案例: 进阶案例: 继上文数据可视化案例,今天学习用pyecharts练习数据可视化案例2-构建地图。 简单案例: 首先构建一个简单的地图。 代码: import json from pyecharts.charts import MapmapMap() data[…

培训学校课程管理系统-计算机毕设Java|springboot实战项目

🍊作者:计算机毕设残哥 🍊简介:毕业后就一直专业从事计算机软件程序开发,至今也有8年工作经验。擅长Java、Python、微信小程序、安卓、大数据、PHP、.NET|C#、Golang等。 擅长:按照需求定制化开发项目、 源…

大数据面试SQL(八):求连续段的起始位置和结束位置

文章目录 求连续段的起始位置和结束位置 一、题目 二、分析 三、SQL实战 四、样例数据参考 求连续段的起始位置和结束位置 一、题目 有一张表t2_id记录了id,id不重复,但是会存在间断,求出连续段的起始位置和结束位置。 样例数据&…

结构体structure、共用体union

目录 结构体 结构体类型的定义形式 结构体类型的大小 内存计算例子 共用体union 用共用体判断大小端 结构体和共用体对比 qsort() 结构体 结构体类型——用来描述复杂数据的一种数据类型 构造类型(用户自定义类型) struc…