2026a

# ty_symmlq

求解线性系统 - 对称的 LQ 方法

函数库: TyMath

# 语法

x = ty_symmlq(A,b)
x = ty_symmlq(A,b,tol)
x = ty_symmlq(A,b,tol,maxit)
x = ty_symmlq(A,b,tol,maxit,M)
x = ty_symmlq(A,b,tol,maxit,M1,M2)
x = ty_symmlq(A,b,tol,maxit,M1,M2,x0)
x,flag,relres,iter,resvec,resveccg = ty_symmlq(__;verbose)

# 说明

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

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

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

x = ty_symmlq(A,b,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解关于 y 的方程组来计算 x，其中且。该算法不显式形成 H。使用预条件子矩阵可以改善问题的数值属性和计算的效率。

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

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

x,flag,relres,iter,resvec,resveccg = ty_symmlq(___,verbose) 指定诊断信息显示指示，指定 verbose 为 false 时，将不会返回诊断信息，而是返回诊断数据，这包括标志 flag，指示算法是否成功收敛。当 flag = 0 时，收敛成功；解中的残差 relres。如果 flag 为 0，则 relres <= tol；计算出 x 时的迭代次数 iter；残差范数向量（包括第一个残差 norm(b-A*x0)）resvec 以及共轭梯度残差范数向量 resveccg。示例

提示

不同电脑可能结果有差异。

# 示例

线性系统的迭代解

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

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

using TyMath
n = 400
on = ones(n)
A = spdiagm(n, n, -1 => -2 * on[1:(end - 1)], 0 => 4 * on, 1 => -2 * on[2:end])
b = full(sum(A; dims=2));

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

x = ty_symmlq(A,b);

symmlq 在迭代 20 停止，而没有收敛到所需容差 1.0e-6。因为已达到最大迭代数。迭代返回的 (数目 20) 的相对残差为 0.04545454545454513。

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

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

x = ty_symmlq(A,b,1e-4,250);

symmlq 在解的迭代 199 处收敛，并且相对残差为 1.456530523957716e-14。

使用 ty_symmlq 返回求解过程信息

通过指定 verbose 为 false 以关闭诊断信息显示并返回有关求解过程的信息：

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

using TyMath
A = delsq(numgrid("S",102));

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

b = sum(A,dims=2);

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

tol = 1e-12;
maxit = 100;

使用 ty_symmlq 根据请求的容差和迭代次数求解。指定 verbose 为 false 以关闭诊断信息显示并返回有关求解过程的信息：

x 是计算 A*x = b 所得的解；
fl0 是指示算法是否收敛的标志；
rr0 是计算的解 x 的相对残差；
it0 是计算 x 时所用的迭代次数；
rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量；
rvcg0 是的共轭梯度残差历史记录组成的向量。

x,fl0,rr0,it0,rv0,rvcg0 = ty_symmlq(A,b,tol,maxit,verbose=false);

fl0

fl0 = 1

rr0

rr0 = 0.003131755925661587

it0

it0 = 100

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

提供初始估计值

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

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

using TyMath
n = 900
e = ones(n)
A = spdiagm(n, n, -1 => e[2:end], 0 => 2 * e, 1 => e[2:end])
b = sum(A; dims=2);

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

maxit = 200
x = ty_symmlq(A, b, [], maxit);

symmlq 在解的迭代 34 处收敛，并且相对残差为 9.481042389826163e-7。

x0 = 0.99 * e
x = ty_symmlq(A, b, [], maxit, [], [], x0);

symmlq 在解的迭代 6 处收敛，并且相对残差为 8.719318069818795e-7。

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

返回中间结果

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

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

R = Vector{Float64}(undef, 4)
X = Matrix{Float64}(undef, size(A, 1), 4)
x0 = zeros(size(A, 2))
tol = 1e-8
maxit = 100
for k in 1:4
    x, flag, relres = ty_symmlq(A, b, tol, maxit, [], [], x0; verbose=false)
    X[:, k] = x
    R[k] = relres
    x0 = x
end

R = 4-element Vector{Float64}:
 4.166699708427977e-8
 9.741482636943538e-9
 9.741482636943538e-9
 9.741482636943538e-9

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

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

通过为 ty_symmlq 提供用来计算 A*x 的函数句柄（而非系数矩阵 A）来求解线性系统。

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

using TyMath
A = TestArrays.wilk(21)

A = 21×21 Matrix{Float64}:

 10.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  1.0  9.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  1.0  8.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  1.0  7.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  1.0  6.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  1.0  5.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  1.0  4.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  1.0  3.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  2.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  2.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  3.0  1.0  0.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  4.0  1.0  0.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  5.0  1.0  0.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  6.0  1.0  0.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  7.0  1.0  0.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  8.0  1.0   0.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  9.0   1.0
  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  10.0

Wilkinson 矩阵有特殊的结构，因此您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。此外，只有主对角线具有不等于 1 的非零值。

表达式 Ax 变为：

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

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

afun = x -> vcat(0, x[1:20]) + vcat(10:-1:0, 1:10) .* x + vcat(x[2:21], 0)

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

b = ones(21)
tol = 1e-12
maxit = 50
x1 = ty_symmlq(afun,b,tol,maxit)

symmlq 在解的迭代 10 处收敛，并且相对残差为 1.974021674086588e-15。
21-element Vector{Float64}:
 0.09101018481309166
 0.08989815186908856
 0.09990644836511509
 0.11085026120999211
 0.12414172316494042
 0.14429939980036552
 0.1543612778332317
 0.23825548886670742
 0.13087225556664606
 0.5000000000000004
 0.369127744433354
 0.5000000000000006
 0.13087225556664614
 0.2382554888667074
 0.15436127783323172
 0.14429939980036552
 0.12414172316494042
 0.11085026120999211
 0.09990644836511506
 0.08989815186908859
 0.09101018481309164

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

afun(x1)

ans = 21-element Vector{Float64}:
 1.000000000000005
 1.0000000000000038
 1.0000000000000013
 1.0000000000000002
 1.0
 0.9999999999999998
 0.9999999999999998
 1.0
 1.0
 1.0000000000000004
 1.0000000000000009
 1.0000000000000007
 1.0000000000000002
 1.0
 0.9999999999999999
 0.9999999999999998
 1.0000000000000002
 1.0000000000000002
 1.000000000000001
 1.000000000000004
 1.000000000000005

# 输入参数

A — 系数矩阵

矩阵 | 函数句柄

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

将 A 指定为函数句柄

您可以选择将系数矩阵指定为函数句柄而不是矩阵。函数句柄返回矩阵向量乘积，而不是构建整个系数矩阵，从而使计算更加高效。

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

复数支持： 是

b — 线性方程的右侧

列向量

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

复数支持： 是

tol — 方法容差

[] 或 1e-6 （默认） | 正标量

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

maxit — 最大迭代次数

[] 或 min(size(A,1),20) （默认） | 正整数标量

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

M, M1, M2 — 预条件子矩阵（以单独参数指定）

eye(size(A)) （默认） | 矩阵 | 函数句柄

预条件子矩阵，指定为由矩阵或函数句柄组成的单独参数。您可以指定预条件子矩阵 M 或其矩阵因子 M = M1*M2 来改进线性系统的数值方面，使 ty_symmlq 更容易快速收敛。

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

将 M 指定为函数句柄

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

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

复数支持： 是

x0 — 初始估计值

[] 或由零组成的列向量（默认） | 列向量

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

复数支持： 是

# 输出参数

x — 线性系统的解

列向量

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

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

flag — 收敛标志

标量

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

标志值	收敛
0	成功 - ty_symmlq 在 maxit 次迭代内收敛至所需容差 tol。
1	失败 - ty_symmlq 执行了 maxit 次迭代，但未收敛。
2	失败 - 预条件子矩阵 M 或 M = M1*M2 为病态。
3	失败 - ty_symmlq 在经过两次相同的连续迭代后已停滞。
4	失败 - 由 ty_symmlq 算法计算的标量数量之一变得太小或太大，无法继续计算。
5	失败 - 预条件子矩阵 M 不是对称正定矩阵。

relres — 相对残差

标量

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

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

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

iter — 迭代编号

标量

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

resvec — 残差

向量

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

resveccg — 共轭梯度残差范数

向量

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

# 详细信息

对称的 LQ 方法

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

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

# 提示

大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)；
您可以使用矩阵重新排序函数来置换系数矩阵的行和列，并在系数矩阵被分解以生成预条件子时最小化非零值的数量。这可以减少后续求解预条件线性系统所需的内存和时间。

# 参考文献

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

# 另请参阅