MATLAB中lsqr函数用法

于 2024-08-15 08:31:43 发布
阅读量628 收藏 17
点赞数 18
分类专栏： Matlab 文章标签： matlab 算法 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jk_101/article/details/133172804
版权

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 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×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

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

        现在，通过为 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×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

最小二乘法

        ​最小二乘 (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.

jk_101
关注
专栏目录
LSMR：最小二乘问题的迭代算法：由 David Fong 和 Michael Saunders 开发的 LSQR 算法的官方继承者。-matlab开发
05-30
提出了求解线性系统和线性最小二乘系统的迭代方法。 该方法基于 Golub-Kahan 双对角化过程。 它在分析上等效于应用于正规方程的 MINRES 标准方法。 与 LSQR 相比，尽早终止 LSMR 更安全。 有关 LSMR 的详细信息可以在http://www.stanford.edu/group/SOL/software/lsmr.html http://www.stanford.edu/~clfong/lsmr.html
matlab中lsqr函数用法
She_Said的博客
11-05 5327
lsqr LSQR 方法 语法 x = lsqr(A,b) lsqr(A,b,tol) lsqr(A,b,tol,maxit) lsqr(A,b,tol,maxit,M) lsqr(A,b,tol,maxit,M1,M2) lsqr(A,b,tol,maxit,M1,M2,x0) [x,flag] = lsqr(A,b,tol,maxit,M1,M2,x0) [x,flag,relres] ...
matlab 空间 lsq,MATLAB 和 Simulink 在空间系统领域的应用
weixin_34392227的博客
03-18 233
请选择其一AlabamaAlaska美属萨摩亚APO/FPO AAAPO/FPO AEAPO/FPO APArizonaArkansasCaliforniaCaroline IslandsColoradoConnecticutDelawareDistrict of ColumbiaFlorida格鲁吉亚关岛HawaiiIdahoIllinoisIndianaIowaKansasKentuckyLo...
Matlab中使用lsqcurvefit进行多变量进行最小二乘计算
weixin_48966416的博客
10-03 3131
在matlab中使用lsqcurvefit进行最小二乘拟合时，遇到多变量的情况，如何进行书写和拟合
LSQR: Sparse Equations and Least Squares
技术、思维
08-16 2954
http://www.stanford.edu/group/SOL/software/lsqr.html LSQR: Sparse Equations and Least Squares AUTHORS: Chris Paige, Michael Saunders.CONTRIBUTORS: James Howse, Michael Friedlander, John Toml
MATLAB有关lsqcurvefit函数对多个未知参数函数拟合的相关问题
热门推荐
m0_56555154的博客
03-10 1万+
在MATLAB中关于函数拟合的函数主要分为两类。一类是线性拟合，另一类则是非线性拟合。非线性拟合相对来说，难度会比线性拟合更大一些。本篇文章主要是关于如何使用lsqcurvefit函数进行非线性拟合的，从而得到函数中的多个未知量。 1、首先，我们通过MATLAB命令行窗口输入help lsqcurvefit，对lsqcurvefit函数的使用作一个大致的了解。 Nonlinear least-squares solver——基于最小二乘法的非线性拟合 以下是使用lsqcurvefit的格式 x = lsq
LSQR 算法
chrisqinxz的博客
09-08 5355
http://www.stanford.edu/group/SOL/software/lsqr.html   LSQR: Sparse Equations and Least Squares AUTHORS: Chris Paige, Michael Saunders. CONTRIBUTORS: James Howse, Michael Friedlander, John Tomlin,...
LSQR.rar_LSQR_lsqr matlab_matlab lsqr_under determined
07-15
LSQR算法是求解大型线性系统的一种高效方法，尤其适用于处理...如果你需要在MATLAB中使用LSQR，可以将FORTRAN代码转换为MATLAB函数，或者直接调用MATLAB内置的`lsqr`函数，该函数已经封装了LSQR算法，方便用户使用。
matlab的lsqr
最新发布
09-22
lsqr是MATLAB中的一个函数，用于求解线性方程组A*x=b的最小二乘解。它采用了迭代的方法，基于正交投影和最小二乘逼近的原理。lsqr函数具有多个输入和输出参数，可以根据不同的需求进行设置。其中，输入参数包括矩阵A...
基于matlab与遗传优化算法的LQR控制器优化系统设计与实现
05-01
通过这个项目，不仅可以掌握LQR控制器的设计原理和MATLAB的使用技巧，还能深入理解遗传优化算法在解决复杂优化问题中的优势。同时，这种结合传统控制理论与现代优化方法的实践，对于提升控制系统性能和应对不确定性...
LSQR C++源码
12-10
This is the C++ version of LSQR as translated from the C version by John Tomlin Tke license conditions are the same as described in the file lsqr_c.README included in this distribution. "make" will produce test program "testit". Running testit should produce output similar to "testit_output"
lsqr最小二乘分解
04-29
地震走时层析成像反演LSQR程序 单元数：9  12 单元边长：3.0  5.0m 参数：速度（km/s) 单边激发单边接收 (左发右收) 数据：12  12
matlab基本函数
freshman_kx的博客
04-05 7880
*matlab基本函数 1.生成矩阵 （1）直接法：a = [1,2,3;4,5,6;7,8,9] （2）冒号一维矩阵： a = 开始：步长：结束，步长为1可省略 b = 1:1:10; b = 1:10; （与上面等价） （3）函数生成： linspace(开始，结束，元素个数)，等差生成指定元素数的一维矩阵，省略个数则生成100个 c = linspace(0,10,5); （4）特殊矩阵： zeros函数:产生零矩阵 ones函数：产生全1矩阵 eye函数：产生对角线为1的矩阵 rand函数：
MATLAB小技巧（30）非线性拟合 lsqcurefit
sinat_34897952的博客
07-07 1408
MATLAB进行图像处理相关的学习是非常友好的，可以从零开始，对基础的图像处理都已经有了封装好的许多可直接调用的函数，这个系列文章的话主要就是介绍一些大家在MATLAB中常用一些概念函数进行例程演示！在科学技术的各领域中，不同变量之间的相互关系通常用函数描述。有些函数关系可由经典理论分析推导得出，不仅为进一步的分析研究工作提供理论基础，也可以方便的解决实际工程问题。但是，很多工程问题难以直接推导出变量之间的函数表达式，或者即使能得出表达式，公式也十分复杂，不利于进一步的分析与计算。由于研究需要，希望能得到这
LSQR算法的程序设计（带阻尼与无阻尼）
weixin_44498266的博客
07-13 4010
LSQR算法是计算大型稀疏线性方程组的算法，由Paige和Saunders于1982年提出。该方法主要为求解以下线性方程组（A为m*n的矩阵，m>n），且保持二阶残差范数最小： Ax=bMin∥Ax−b∥2 Ax=b\\ Min\parallel Ax-b \parallel_{2} Ax=bMin∥Ax−b∥2​ 对于传统的最小二乘问题，可以直接求出方程的解，但是对于大型稀疏矩阵，往往是病态问题（超定方程组或欠定方程组），因此，LSQR方法不是直接求解x，而是遵守最小二乘的原则，在法方程限定的空间中
基于MATLAB的稀疏矩阵与符号求解线性方程组（附完整代码）
forest_LL的博客
04-19 8914
前言 前面文章已经更新了至少三种求解线性方程组的MATLAB指令，分别是： 逆矩阵法：x=inv(A)*b，或者写成 x=A^(-1)*b 伪逆法：x=pinv(A)*b 左除法：x=A\b 本节将更新另外两种方法，符号矩阵法与稀疏矩阵求解法。 一. 符号解法 在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，MATLAB格式如下： linsolve(sym(A),sym(b)) %结果默认小数表示 还有另外一种函数也可以，形式如下： X=s..
Python Numpy包中的linalg.lstsq()函数与Matlab中 \ 的区别
qq_44341554的博客
05-29 993
Python Numpy包中的linalg.lstsq()函数与Matlab中 \ 的区别 起因：将Matlab的代码转换为Python形式的过程中发现这一情况。 经过：根据Python Numpy包中各种操作符与Matlab的对应关系将程序转换完毕后。在测试的时候发现Python结果输出与Matlab的结果并不一致。于是开启调试模式，在分步调试过程中发现是在进行 \ 运算时，Python和Matlab开始出现区别，于是进行分析百度。 已知：A \ B 和 linalg.lstsq(A，B)[0]是对线性方
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值