Main Content

本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

symmlq

求解线性方程组 - 对称的 LQ 方法

说明

示例

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

示例

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

示例

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

示例

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

示例

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

示例

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

示例

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

示例

[x,flag,relres] = symmlq(___) 还会返回计算的解中的残差。如果 flag0,则 relres <= tol

示例

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

示例

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

示例

[x,flag,relres,iter,resvec,resveccg] = symmlq(___) 还会在每次迭代中返回共轭梯度残差范数向量。

示例

全部折叠

使用采用默认设置的 symmlq 求解系数矩阵为方阵的线性方程组,然后在求解过程中调整使用的容差和迭代次数。

创建一个稀疏三对角矩阵 A 作为系数矩阵。使用 A 的稠密行总和作为 Ax=b 右侧的向量,以便 x 的预期解是由 1 组成的向量。

n = 400; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = full(sum(A,2));

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

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

默认情况下,symmlq 使用 20 次迭代和容差 1e-6,对于此矩阵,算法无法在 20 次迭代后收敛。由于残差的数量级为 1e-2,显然需要更多迭代。您还可以减小容差,使算法更容易收敛。

使用容差 1e-4 和 250 次迭代再次求解方程组。

x = symmlq(A,b,1e-4,250);
symmlq converged at iteration 199 to a solution with relative residual 1.5e-14.

检查使用指定了预条件子矩阵的 symmlq 来求解线性方程组的效果。

创建一个对称正定带状系数矩阵。

A = delsq(numgrid('S',102));

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 100;

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

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

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

  • rr0 是计算的解 x0 的残差。

  • it0 是计算出 x0 时的迭代次数。

  • rv0Ax-b 的残差历史记录组成的向量。

  • rvcg0ATAx-ATb 的共轭梯度残差历史记录组成的向量。

[x0,fl0,rr0,it0,rv0,rvcg0] = symmlq(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.0031
it0
it0 = 100

symmlq 未在请求的 100 次迭代内收敛至请求的容差 1e-12,因此 fl0 为 1。

为了帮助收敛,您可以指定预条件子矩阵。由于 A 是对称矩阵,请使用 ichol 生成预条件子 M=LLT。指定 'ict' 选项以使用阈值调降不完全 Cholesky 分解,并指定对角线偏移量值 1e-6 以避免非正主元。通过指定 LL' 作为 symmlq 的输入,求解预条件方程组。

setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = symmlq(A,b,tol,maxit,L,L');
fl1
fl1 = 0
rr1
rr1 = 2.5636e-15
it1
it1 = 3

使用 ichol 预条件子极大地改进了问题的数值属性,symmlq 能够快速收敛。输出 rv1(1)norm(b),输出 rv1(end)norm(b-A*x1)

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

semilogy(0:length(rvcg0)-1,rvcg0/norm(b),'-o')
hold on
semilogy(0:length(rvcg1)-1,rvcg1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ICHOL preconditioner','Tolerance','Location','SouthEast')
xlabel('Iteration number')
ylabel('Relative residual')

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

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

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

使用 symmlq 求解 Ax=b 两次:一次是使用默认的初始估计值,一次是使用解的良好初始估计值。对这两个解都使用 200 次迭代,并将初始估计值指定为所有元素均等于 0.99 的向量。

maxit = 200;
x1 = symmlq(A,b,[],maxit);
symmlq converged at iteration 34 to a solution with relative residual 9.5e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = symmlq(A,b,[],maxit,[],[],x0);
symmlq converged at iteration 6 to a solution with relative residual 8.7e-07.

在这种情况下,提供初始估计值可以使 symmlq 更快地收敛。

返回中间结果

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

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

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

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

通过为 symmlq 提供用来计算 A*x 的函数句柄(而非系数矩阵 A)来求解线性方程组

gallery 生成的 Wilkinson 测试矩阵之一是 21×21 三对角矩阵。预览该矩阵。

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

Wilkinson 矩阵有特殊的结构,因此您可以用函数句柄来表示 A*x 运算。由于 A 的每行都乘以 x 中的元素,因此只有少数结果为非零值(对应于三对角线上的非零值)。此外,只有主对角线具有不等于 1 的非零值。表达式 Ax 变为:

[1010001910001810017100161001510014100130001000110][x1x2x3x4x5x21]=[10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21]

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

[0+10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21+0]=[0x1x20]+[10x19x210x21]+[x2x210]

在 MATLAB® 中,编写一个函数来创建这些向量并将它们相加,从而给出 A*x 的值:

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

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

现在,通过为 symmlq 提供用于计算 A*x 的函数句柄,求解线性方程组 Ax=b。使用容差 1e-12 和 50 次迭代。

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = symmlq(@afun,b,tol,maxit)
symmlq converged at iteration 10 to a solution with relative residual 4.5e-16.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

检查 afun(x1) 是否产生由 1 组成的向量。

afun(x1)
ans = 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)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

输入参数

全部折叠

系数矩阵,指定为对称矩阵或函数句柄。该矩阵是线性方程组 A*x = b 中的系数矩阵。通常,A 是大型稀疏矩阵或函数句柄,它返回大型稀疏矩阵和列向量的乘积。您可以使用 issymmetric 来确认 A 是对称矩阵。

A 指定为函数句柄

您可以将系数矩阵指定为函数句柄而不是矩阵,以节省计算中使用的内存。函数句柄返回矩阵向量乘积,而不是构建整个系数矩阵,从而使计算更加高效。

要使用函数句柄,请使用函数签名 function y = afun(x)参数化函数说明如何在必要时为函数 afun 提供附加参数。函数调用 afun(x) 必须返回 A*x 的值。

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

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

数据类型: double
复数支持:

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

数据类型: double

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

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

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

M 指定为函数句柄

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

要使用函数句柄,请使用函数签名 function y = mfun(x)参数化函数说明如何在必要时为函数 mfun 提供附加参数。函数调用 mfun(x) 必须返回 M\xM2\(M1\x) 的值。

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

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

数据类型: double
复数支持:

输出参数

全部折叠

线性方程组解,以向量形式返回。该输出给出线性方程组 A*x = b 的近似解。如果计算成功 (flag = 0),则 relres 小于或等于 tol

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

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

标志值

收敛

0

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

1

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

2

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

3

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

4

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

5

失败 - 预条件子矩阵 M 不是对称正定矩阵。

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

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

  • 共轭梯度残差等于 norm(A'*A*X - A'*b)。与相对残差相比,此残差会使 symmlq 较少收敛。resveccg 输出会跟踪此残差在所有迭代上的历史记录。

数据类型: double

迭代编号,以标量形式返回。此输出指示计算出 x 的解时的迭代次数。

数据类型: double

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

数据类型: double

共轭梯度残差范数,以向量形式返回。resveccg 中元素的数量等于迭代次数。

数据类型: double

详细信息

全部折叠

对称的 LQ 方法

MINRES 和 SYMMLQ 方法是支撑共轭梯度法 PCG 的 Lanczos 方法的变体。像 PCG 一样,系数矩阵仍需是对称矩阵,但 MINRES 和 SYMMLQ 允许它是不定矩阵(不要求所有特征值都必须为正值)。这是通过避免 Lanczos 方法中通常存在的隐式 LU 分解来实现的,当遇到零主元时,该方法容易出现故障。

MINRES 最小化 2-范数残差,而 SYMMLQ 使用 LQ 分解求解投影方程组,并使残差与所有先前的残差正交 [1]

提示

  • 大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)。您可以使用 equilibrate 来改进 A 的条件数,它本身就能使大多数迭代求解器更容易收敛。但如果您随后会对经平衡处理的矩阵 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, "Solution of Sparse Indefinite Systems of Linear Equations." SIAM J. Numer. Anal., Vol.12, 1975, pp. 617-629.

扩展功能

在 R2006a 之前推出