Main Content

本页翻译不是最新的。点击此处可查看最新英文版本。

lsqr

求解线性系统 - 最小二乘法

说明

示例

x = lsqr(A,b) 尝试使用最小二乘法求解关于 x 的线性系统 A*x = blsqr 求最小化 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 的方程组 AM1y=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 的残差。如果 flag0,则 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。输出显示包括相对残差 b-Axb 的值。

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")

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Least-squares residual, Relative residual.

检查使用指定了预条件子矩阵的 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 时所用的迭代序号。

  • rvb-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 的非对角线元。通过将 LU 指定为 lsqrM1M2 输入,求解 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')

Figure contains an axes object. The axes object with xlabel Iteration number, ylabel Relative residual contains 3 objects of type line, constantline. These objects represent No preconditioner, ILU preconditioner, Tolerance.

检查向 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*xA'*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 x=[1020019200120010010200110][x1x2x3x21]=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21].

结果向量可以写为三个向量的和:

A x=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21]=[0x1x2x20]+[10x19x29x2010x21]+2[x2x3x210]

同样,AT x 的表达式变为:

AT x=[1010029100210020010100210][x1x2x3x21]=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21].

AT x=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21]=2[0x1x2x20]+[10x19x29x2010x21]+[x2x3x210].

在 MATLAB® 中,编写一个函数来创建这些向量并将它们相加,从而根据标志输入给出 A*xA'*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*xA'*x 的函数句柄,求解线性系统 Ax=b。使用容差 1e-6 和 25 次迭代。指定 bA 的行总和,使得 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

输入参数

全部折叠

系数矩阵,指定为矩阵或函数句柄。该矩阵是线性系统 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 使用 BC 中的值来计算 A*xA'*x(具体取决于指定的标志),而不用真正建立整个矩阵。

数据类型: double | function_handle
复数支持:

线性方程的右侧,指定为列向量。b 的长度必须等于 size(A,1)

数据类型: double
复数支持:

方法容差,指定为正标量。计算中使用此输入可在准确度和运行时间之间进行权衡。lsqr 必须在允许的迭代次数内满足容差才能成功。较小的 tol 值意味着解必须更精确才能成功完成计算。

数据类型: double

最大迭代次数,指定为正整数标量。增加 maxit 的值,以允许 lsqr 进行更多迭代,从而满足容差 tol。通常,较小的 tol 值意味着需要更多迭代才能成功完成计算。

预条件子矩阵,指定为由矩阵或函数句柄组成的单独参数。您可以指定预条件子矩阵 M 或其矩阵因子 M = M1*M2 来改进线性系统的数值方面,使 lsqr 更容易快速收敛。对于系数方阵,您可以使用不完全矩阵分解函数 iluichol 来生成预条件子矩阵。您还可以在分解之前使用 equilibrate 来改进系数矩阵的条件数。有关预条件子的详细信息,请参阅线性方程组的迭代方法

lsqr 将未指定的预条件子视为单位矩阵。

M 指定为函数句柄

您可以选择将 MM1M2 中的任一个指定为函数句柄而不是矩阵。函数句柄执行矩阵向量运算,而不是构建整个预条件子矩阵,从而使计算更加高效。

要使用函数句柄,首先创建签名为 function y = mfun(x,opt) 的函数。参数化函数说明如何在必要时为函数 mfun 提供附加参数。函数 mfun 必须满足下列条件:

  • mfun(x,'notransp') 返回 M\xM2\(M1\x) 的值。

  • mfun(x,'transp') 返回 M'\xM1'\(M2'\x) 的值。

可接受的函数的一个示例是:

function y = mfun(x,opt,a,b)  
if strcmp(opt,'notransp')
    y = x.*a;
else
    y = x.*b;
end
end
在此示例中,函数 mfun 使用 ab 来计算 M\x = x*aM'\x = x*b(具体取决于指定的标志),而不用真正建立整个矩阵 M

数据类型: double | function_handle
复数支持:

初始估计值,指定为长度等于 size(A,2) 的列向量。如果您能为 lsqr 提供比默认的零向量更合理的初始估计值 x0,则它可以节省计算时间并帮助算法更快地收敛。

数据类型: double
复数支持:

输出参数

全部折叠

线性系统的解,以列向量形式返回。该输出给出线性系统 A*x = b 的近似解。

  • 如果 flag0relres <= tol,则 xA*x = b 的相容解。

  • 如果 flag0relres > tol,则 x 是最小化 norm(b-A*x) 的最小二乘解。在这种情况下,lsvec 输出包含 x 的缩放标准方程误差。

每当计算不成功时 (flag ~= 0),lsqr 返回的解 x 是在所有迭代中计算出的最小范数残差。

收敛标志,返回下表中的标量值之一。收敛标志指示计算是否成功,并区分几种不同形式的失败。

标志值

收敛

0

成功 - lsqrmaxit 次迭代内收敛至所需容差 tol

1

失败 - lsqr 执行了 maxit 次迭代,但未收敛。

2

失败 - 预条件子矩阵 MM = M1*M2 为病态。

3

失败 - lsqr 在经过两次相同的连续迭代后已停滞。

4

失败 - 由 lsqr 算法计算的标量数量之一变得太小或太大,无法继续计算。

相对残差,以标量形式返回。相对残差表明返回的答案 x 的准确度。lsqr 跟踪求解过程中每次迭代的相对残差和最小二乘残差,当任一残差满足指定的容差 tol 时,算法收敛。relres 输出包含收敛的残差的值,即相对残差或最小二乘残差:

  • 相对残差等于 norm(b-A*x)/norm(b),通常是在 lsqr 收敛时满足容差 tol 的残差。resvec 输出会跟踪此残差在所有迭代上的历史记录。

  • 最小二乘残差等于 norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro')。与相对残差相比,此残差会使 lsqr 较少收敛。lsvec 输出会跟踪此残差在所有迭代上的历史记录。

迭代序号,以标量形式返回。此输出指示计算出 x 的解时所用的迭代序号。

数据类型: double

残差,以向量形式返回。残差 norm(b-A*x) 揭示对于给定的 x 值,算法接近收敛的程度。resvec 中元素的数量等于迭代次数。您可以检查 resvec 的内容,以帮助决定是否更改 tolmaxit 的值。

数据类型: double

缩放标准方程误差,以向量形式返回。对于每次迭代,lsvec 包含缩放标准方程残差 norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro') 的一个估计值。lsvec 中元素的数量等于迭代次数。

详细信息

全部折叠

最小二乘法

最小二乘 (LSQR) 算法是对矩形矩阵的共轭梯度 (CG) 方法的改进。分析表明,A*x = b 的 LSQR 产生的残差与标准方程 A'*A*x = A'*b 的 CG 相同,但 LSQR 具有更好的数值属性,因此通常更可靠 [1]

最小二乘法是唯一可处理矩形和不相容系数矩阵的迭代线性系统求解器。

提示

  • 大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)。当 A 是方阵时,您可以使用 equilibrate 来改进其条件数,它本身就能使大多数迭代求解器更容易收敛。但如果您随后会对经平衡处理的矩阵 B = R*P*A*C 进行因式分解,使用 equilibrate 还可以获得质量更好的预条件子矩阵。

  • 您可以使用矩阵重新排序函数(如 dissectsymrcm)来置换系数矩阵的行和列,并在系数矩阵被分解以生成预条件子时最小化非零值的数量。这可以减少后续求解预条件线性系统所需的内存和时间。

参考

[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.

扩展功能

版本历史记录

在 R2006a 之前推出