求解一个矩形线性方程组lsqr使用默认设置，然后调整解决过程中使用的公差和迭代次数。

创建一个随机稀疏矩阵一个50%的密度。也可以创建一个随机向量b的右边 $\mathrm{斧头}＝\mathit{b}$．

rng默认的A = sprand(400,300，.5);B = rand(400,1);

解决 $\mathrm{斧头}＝\mathit{b}$使用lsqr．输出显示包括相对残余误差的值 $\frac{为\mathit{b}-\mathrm{斧头}为}{为\mathit{b}为}$．

x = lsqr(A,b);

LSQR在迭代20处停止，没有收敛到期望的容差1e-06，因为达到了最大迭代次数。返回的迭代(编号20)的相对残差为0.26。

默认情况下lsqr使用20次迭代，公差为1 e-6，但算法无法在这个矩阵的这20次迭代中收敛。由于残差仍然很大，这是一个很好的指标，表明需要更多的迭代(或预处理矩阵)。您还可以使用更大的公差，使算法更容易收敛。

再用一个容差解系统1的军医70次迭代。指定六个输出以返回相对残差relres的计算解，以及残差历史resvec以及最小二乘残差历史lsvec．

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

Flag = 0

自国旗为0，表示该算法在指定迭代次数内能够满足所期望的容错。您通常可以一起调整公差和迭代次数，以这种方式在速度和精度之间进行权衡。

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

relres

Relres = 0.2625

Lsres = lsvec(end)

Lsres = 2.7640e-04

这些剩余规范表明x是最小二乘解，因为relres是否小于规定的公差1的军医．由于线性系统不存在一致解，求解器能做的最好的是使最小二乘残差满足公差。

绘制剩余历史。相对残差resvec迅速达到最小值，不能再前进，而最小二乘残差lsvec在后续迭代中继续最小化。

N = length(resvec);lsvec semilogy (0: n - 1,“- o”0: n - 1、resvec“o”)传说(“最小二乘残差”，“相对剩余”）

检查使用预处理矩阵的效果lsqr解一个线性方程组。

加载west0479，一个真实的479 × 479非对称稀疏矩阵。

负载west0479A = west0479;

定义b这才是真正的解 $\mathrm{斧头}＝\mathit{b}$都是1的向量。

b = sum(A,2);

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

Tol = 1e-12;Maxit = 20;

使用lsqr在要求的公差和迭代次数上找到解决方案。指定六个输出以返回关于解决方案过程的信息:

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

Fl = 1

rr

Rr = 0.0017

它

It = 20

自Fl = 1时，算法在最大迭代次数内没有收敛到指定的公差。

为了帮助处理缓慢的收敛，您可以指定一个预处理矩阵。自一个是非对称的，用ilu生成前置条件 $\mathit{米}＝\mathit{l}\mathit{U}$因式分解。指定一个落差公差以忽略值小于的非对角线项1 e-6．解预条件方程组 ${\mathrm{我}}^{-1}（\mathit{米}\mathit{x}）＝\mathit{b}$为 $\mathit{y}＝\mathrm{Mx}$通过指定l而且U随着M1而且平方米输入lsqr．

设置= struct(“类型”，“ilutp”，“droptol”1 e-6);[L,U] = ilu(A，设置);[x1,fl1,rr1,it1,rv1,lsrv1] = lsqr(A,b,tol,maxit,L,U);fl1

Fl1 = 0

rr1

Rr1 = 7.0954e-14

it1

It1 = 13

使用ilu预处理剂产生的相对残留小于规定的容差1 e-12在第13次迭代中。输出rv1 (1)是规范(b)，输出rv1(结束)是规范(b * x1)．

你可以跟随的进度lsqr通过绘制每次迭代的相对残差。用指定公差的一条线绘制每个溶液的残留历史。

semilogy(0:长度(rv) 1、rv /规范(b),“o”)举行在semilogy(0:长度(rv1) 1, rv1 /规范(b),“o”) yline(托尔,“r——”）;传奇(“没有预调节器”，ILU预处理的，“宽容”，“位置”，“东”)包含(的迭代次数) ylabel (的相对剩余的）

检查供给的效果lsqr对解的初步猜测。

创建一个随机的矩形稀疏矩阵。用每一行的和作为右边的向量 $\mathrm{斧头}＝\mathit{b}$这就是期望的解 $\mathit{x}$是1的向量。

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

使用lsqr来解决 $\mathrm{斧头}＝\mathit{b}$两次:一次是默认的初始猜测，另一次是正确的初始猜测。对两个解决方案都使用75次迭代和默认容差。金宝搏官方网站将第二个解决方案中的初始猜测指定为所有元素都等于0.99的向量。

Maxit = 75;x1 = lsqr(A,b，[]，maxit);

LSQR在迭代64时收敛到相对残差为8.7e-07的解。

x0 = 0.99*ones(size(A,2)，1);x2 = lsqr(A,b，[]，maxit，[]，[]，x0);

LSQR在迭代26时收敛到相对残差为9.6e-07的解。

有了接近预期解的初始猜测，lsqr能够在更少的迭代中收敛。

x0 = 0 (size(A,2)，1);Tol = 1e-8;Maxit = 100;为k = 1:4 [x,国旗,relres] = lsqr (A, b,托尔,麦克斯特[],[],x0);X(:，k) = X;R(k) = relres;X0 = x;结束

求解一个线性方程组lsqr使用一个函数句柄进行计算* x而且‘* x代替系数矩阵一个．

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

画廊(“wilk”+ diag(ones(20,1)，1)

一个=21日×2110 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 9 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 2 5 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⋮

因为这个三对角矩阵有一个特殊的结构，你可以表示这个操作* x使用函数句柄。当一个乘以一个向量，结果向量中的大部分元素都是零。结果中的非零元素对应于的非零三对角元素一个．

表达式 $\mathit{一个}\mathit{x}$就变成:

$\mathit{一个}\mathit{x}＝［\begin{array}{ccccccc}10& 2& 0& \cdots & & \cdots & 0\\ 1& 9& 2& 0& & & ⋮\\ 0& 1& \ddots & 2& 0& & \\ ⋮& 0& 1& 0& \ddots & \ddots & ⋮\\ & & 0& \ddots & 1& \ddots & 0\\ ⋮& & & \ddots & \ddots & \ddots & 2\\ 0& \cdots & & \cdots & 0& 1& 10\end{array}］［\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {\mathit{x}}_{21}\end{array}］＝［\begin{array}{c}{10\mathit{x}}_1+2{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+2{\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20.}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20.}+10{\mathit{x}}_{21}\end{array}］$．

得到的向量可以写成三个向量的和:

$\mathit{一个}\mathit{x}＝［\begin{array}{c}{10\mathit{x}}_1+2{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+2{\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20.}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20.}+10{\mathit{x}}_{21}\end{array}］$＝ $［\begin{array}{c}0\\ {\mathit{x}}_1\\ {\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20.}\end{array}］+［\begin{array}{c}{10\mathit{x}}_1\\ {9\mathit{x}}_2\\ ⋮\\ 9{\mathit{x}}_{20.}\\ 10{\mathit{x}}_{21}\end{array}］+2\cdot ［\begin{array}{c}{\mathit{x}}_2\\ {\mathit{x}}_{3.}\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}］$．

同样，for的表达式 ${\mathit{一个}}^{\mathit{T}}\mathit{x}$就变成:

${\mathit{一个}}^{\mathit{T}}\mathit{x}＝［\begin{array}{ccccccc}10& 1& 0& \cdots & & \cdots & 0\\ 2& 9& 1& 0& & & ⋮\\ 0& 2& \ddots & 1& 0& & \\ ⋮& 0& 2& 0& \ddots & \ddots & ⋮\\ & & 0& \ddots & 1& \ddots & 0\\ ⋮& & & \ddots & \ddots & \ddots & 1\\ 0& \cdots & & \cdots & 0& 2& 10\end{array}］［\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {\mathit{x}}_{21}\end{array}］＝［\begin{array}{c}{10\mathit{x}}_1+{\mathit{x}}_2\\ {2\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20.}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20.}+10{\mathit{x}}_{21}\end{array}］$．

${\mathit{一个}}^{\mathit{T}}\mathit{x}＝［\begin{array}{c}{10\mathit{x}}_1+{\mathit{x}}_2\\ {2\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_{3.}\\ ⋮\\ ⋮\\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20.}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20.}+10{\mathit{x}}_{21}\end{array}］＝2\cdot ［\begin{array}{c}0\\ {\mathit{x}}_1\\ {\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20.}\end{array}］+［\begin{array}{c}{10\mathit{x}}_1\\ {9\mathit{x}}_2\\ ⋮\\ 9{\mathit{x}}_{20.}\\ 10{\mathit{x}}_{21}\end{array}］+［\begin{array}{c}{\mathit{x}}_2\\ {\mathit{x}}_{3.}\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}］$．

在MATLAB®中，编写一个函数来创建这些向量并将它们相加，从而给出的值* x或‘* x，取决于标志输入:

函数Y = fun(x,flag)如果比较字符串(国旗,“notransp”）%计算A*xY = [0;x (1:20)).．.+ ((10: 1:0) ';(1:10)”)。* x.．.+ 2 * (x(2:结束);0);elseif比较字符串(国旗,“透明”）计算A'*xY = 2*[0;x (1:20)).．.+ ((10: 1:0) ';(1:10)”)。* x.．.+ (x(2:结束);0);结束结束

(该函数在示例末尾保存为本地函数。)

现在求解线性方程组 $\mathrm{斧头}＝\mathit{b}$通过提供lsqr使用函数句柄进行计算* x而且‘* x．使用的公差1 e-625次迭代。指定 $\mathit{b}$的行和 $\mathit{一个}$这就是方程的解 $\mathit{x}$是1的向量。

b = full(sum(A,2));Tol = 1e-6;Maxit = 25;X1 = lsqr(@afun,b,tol,maxit)

LSQR在迭代21收敛到相对残差为5.4e-13的解。

x1 =21日×11.0000 - 1.0000 - 1.0000 - 1.0000 - 1.0000 - 1.0000 - 1.0000 - 1.0000

函数Y = fun(x,flag)如果比较字符串(国旗,“notransp”）%计算A*xY = [0;x (1:20)).．.+ ((10: 1:0) ';(1:10)”)。* x.．.+ 2 * (x(2:结束);0);elseif比较字符串(国旗,“透明”）计算A'*xY = 2*[0;x (1:20)).．.+ ((10: 1:0) ';(1:10)”)。* x.．.+ (x(2:结束);0);结束结束