# tfqmr
求解线性方程组 - 无转置拟最小残差法
函数库: TyMath
# 语法
x, flag, relres, iter, resvec = tfqmr(A, b)
_ = tfqmr(A, b, tol)
_ = tfqmr(A, b, tol, maxit)
_ = tfqmr(A, b, tol, maxit, M)
_ = tfqmr(A, b, tol, maxit, M1, M2)
_ = tfqmr(A, b, tol, maxit, M1, M2, x0)
_ = tfqmr(_; verbose=true)
# 说明
x, flag, relres, iter, resvec = tfqmr(A, b) 尝试使用无转置拟最小残差法求解关于 x 的线性方程组 A*x = b。其中 flag 为一个指示算法是否成功收敛的标志。当 flag = 0 时，收敛成功。relres 为相对残差 norm(b-A*x)/norm(b)。如果 flag 为 0，则 relres <= tol。iter 为计算出 x 时的迭代次数。resvec 为每次迭代中的残差范数向量包括第一个残差 norm(b-A*x0)）。您总可以使用 [] 返回部分输出，例如 x = tfqmr(A,b)[1] 。示例
_ = tfqmr(A,b,tol) 指定该方法的容差。默认容差是 1e-6。示例
_ = tfqmr(A,b,tol,maxit) 指定要使用的最大迭代次数。示例
_ = tfqmr(A,b,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解关于 y 的方程组
_ = tfqmr(A,b,tol,maxit,M1,M2) 指定预条件子矩阵 M 的因子，使得 M = M1*M2。
_ = tfqmr(A,b,tol,maxit,M1,M2,x0) 指定解向量 x 的初始估计值。默认值为由零组成的向量。示例
_ = tfqmr(_; verbose = true) 在上述用法的基础上增加了 verbose 变量。如果该变量值为 true（默认值），则 tfqmr 会显示收敛信息：如果尝试成功，tfqmr 会显示一条消息来确认收敛。如果 tfqmr 无法在达到最大迭代次数后收敛或出于任何原因暂停，则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。否则不会显示。示例
# 示例
线性方程组的迭代解
使用采用默认设置的 tfqmr 求解系数矩阵为方阵的线性方程组，然后在求解过程中调整使用的容差和迭代次数。
创建密度为 50% 的随机稀疏矩阵 A。另为
using TyMath
rng = MersenneTwister(1234)
A = sprand(rng, 400, 400, 0.5)
A = A'*A
b = rand(rng, 400);
使用 tfqmr 求解
x, = tfqmr(A, b);
tfqmr 在迭代 40 停止，而没有收敛到所需容差 1.0e-6。因为已达到最大迭代数。迭代返回的 (数目 15) 的相对残差为 0.27221337088664255。
默认情况下，tfqmr 使用 40 次迭代和容差 1e-6，对于此矩阵，算法无法在 40 次迭代后收敛。由于残差仍然很大，这说明需要更多的迭代（或预条件子矩阵）。您也可以使用更大的容差，使算法更容易收敛。
使用容差 1e-4 和 100 次迭代再次求解方程组。
x, = tfqmr(A,b,1e-4,100);
tfqmr 在迭代 200 停止，而没有收敛到所需容差 0.0001。因为已达到最大迭代数。迭代返回的 (数目 15) 的相对残差为 0.27221337088664255。
即使使用更宽松的容差和更多迭代，残差也并未改进。当迭代算法以这种方式停滞时，显然需要预条件子矩阵。
提供初始初始值
检查向 tfqmr 提供解的初始估计值的效果。
创建一个三对角稀疏矩阵。使用每行的总和作为
using TyMath
n = 900
e = vec(ones(n,1))
A = spdiagm(n,n,-1=>e[1:end-1],0=>2 .*e,1=>e[1:end-1])
b = sum(A,dims=2);
使用 tfqmr 求解
maxit = 200
x1, = tfqmr(A,b,1e-6,maxit);
tfqmr 在解的迭代 19 处收敛，并且相对残差为 9.565651220382304e-7。
x0 = 0.99 .* e
x2, = tfqmr(A,b,1e-6,maxit,[],[],x0);
tfqmr 在解的迭代 4 处收敛，并且相对残差为 7.926006665268287e-7。
在这种情况下，提供初始估计值可以使 tfqmr 更快地收敛。
返回中间结果
您还可以通过在 for 循环中调用 tfqmr 来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代，并存储计算出的解。然后，将该解用作下一批迭代的初始向量。
例如，以下代码会循环执行四次，每次执行 100 次迭代，并在 for 循环中每通过一次后均存储解向量：
x0 = zeros(size(A,2),1)
tol = 1e-8
maxit = 100
X = zeros(size(A,2),4)
R = zeros(4)
for k = 1:4
x, flag, relres = tfqmr(A,b,tol,maxit,[],[],x0,verbose=false)
X[:,k] = x
R[k] = relres
x0 = x
end
X[:,k] 是在 for 循环的第 k 次迭代时计算的解向量，R[k] 是该解的相对误差。
使用函数句柄代替数值矩阵
通过为 tfqmr 提供用来计算 A*x 的函数句柄（而非系数矩阵 A）来求解线性方程组。
使用 wilkinson 函数生成 Wilkinson 21×21 三对角测试矩阵。预览该矩阵。
using TyMath
A = Int.(wilkinson(21))
A = 21×21 Matrix{Int64}:
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
0 0 0 0 0 0 0 0 0 1 0 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 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 3 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 5 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 7 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 9 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 10
Wilkinson 矩阵有特殊的结构，因此您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。此外，只有主对角线具有不等于 1 的非零值。
表达式 Ax 变为：
结果向量可以写为三个向量的和：
在 Syslab 中，编写一个函数来创建这些向量并将它们相加，从而给出 A*x 的值。
function afun(x)
return [0;x[1:20]]+[10:-1:0;1:10].*x+[x[2:21];0]
end
现在，通过为 tfqmr 提供用于计算 A*x 的函数句柄，求解线性方程组 Ax=b。使用容差 1e-12 和 50 次迭代。
b = ones(21)
tol = 1e-12
maxit = 50
x1 = tfqmr(afun,b,tol,maxit)[1]
tfqmr 在解的迭代 10 处收敛，并且相对残差为 1.3876413177894163e-14。
21×1 Matrix{Float64}:
0.09101018481309114
0.08989815186908816
0.09990644836511493
0.11085026120999206
0.12414172316494035
0.14429939980036555
0.15436127783323192
0.2382554888667065
0.1308722555666494
0.49999999999999045
0.36912774443337176
0.49999999999998945
0.13087225556665014
0.23825548886670603
0.15436127783323206
0.14429939980036546
0.12414172316494036
0.11085026120999206
0.09990644836511503
0.08989815186908838
0.09101018481309138
检查 afun(x1) 是否产生由 1 组成的向量。
afun(x1)
ans = 21×1 Matrix{Float64}:
0.9999999999999996
0.9999999999999994
0.9999999999999997
0.9999999999999997
0.9999999999999998
1.0
0.9999999999999998
1.0000000000000009
0.9999999999999958
1.0000000000000115
0.9999999999999799
1.0000000000000113
0.9999999999999957
1.0000000000000004
0.9999999999999997
0.9999999999999997
0.9999999999999997
0.9999999999999998
1.0000000000000007
1.0000000000000018
1.0000000000000022
使用指定了预条件子的 tfqmr
检查使用指定了预条件子矩阵的 tfqmr 来求解线性方程组的效果。
using TyMath
function west0479()
n = 479
A = zeros(n, n)
for i in 1:n
A[i, i] = i
if i > 1
A[i, i-1] = sqrt(i * (i - 1))
A[i-1, i] = sqrt(i * (i - 1))
end
end
return A
end
生成 west0479 矩阵，它是一个非对称的 479×479 实稀疏矩阵。
A = west0479();
对矩阵 A 的每一行进行求和，并转换成向量形式。
b = vec(sum(A,dims=2));
设置容差和最大迭代次数。
tol = 1e-12;
maxit = 20;
使用 tfqmr 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息：
x0 是计算 A*x0 = b 所得的解；
fl0 是指示算法是否收敛的标志；
rr0 是计算的解 x0 的残差；
it0 是计算出 x0 时的迭代次数；
rv0 是 ‖Ax−b‖ 的残差历史记录组成的向量。
x0,fl0,rr0,it0,rv0 = tfqmr(A, b,tol,maxit);
fl0
tfqmr stopped at iteration 40 without converging to the desired tolerance 1.0e-12 because the maximum number of iterations was reached.The iterate returned (number 3) has relative residual 0.024956282808941715.
1
# 输入参数
A - 系数矩阵矩阵 | 函数句柄
系数矩阵，指定为方阵或函数句柄。该矩阵是线性方程组 A*x = b 中的系数矩阵。通常，A 是大型稀疏矩阵或函数句柄，它返回大型稀疏矩阵和列向量的乘积。
将 A 指定为函数句柄
您可以选择将系数矩阵指定为函数句柄而不是矩阵。函数句柄返回矩阵向量乘积，而不是构建整个系数矩阵，从而使计算更加高效。
要使用函数句柄，请使用函数签名 function afun(x)。参数化函数说明如何在必要时为函数 afun 提供附加参数。函数调用 afun(x) 必须返回 A*x 的值。
数据类型： Float64 | Function
复数支持： 是
b - 线性方程的右侧列向量
线性方程的右侧，指定为列向量。b 的长度必须等于 size(A,1)。
数据类型： Float64
复数支持： 是
tol - 方法容差1e-6（默认） | 正标量
方法容差，指定为正标量。计算中使用此输入可在准确度和运行时间之间进行权衡。tfqmr 必须在允许的迭代次数内满足容差才能成功。较小的 tol 值意味着解必须更精确才能成功完成计算。
数据类型： Float64
maxit - 最大迭代次数min(size(A,1),20)（默认） | 正整数标量
最大迭代次数，指定为正整数标量。增加 maxit 的值，以允许 tfqmr 进行更多迭代，从而满足容差 tol。通常，较小的 tol 值意味着需要更多迭代才能成功完成计算。
M, M1, M2 - 预条件子矩阵（以单独参数指定）eye(size(A))（默认） | 矩阵 | 函数句柄
预条件子矩阵，指定为由矩阵或函数句柄组成的单独参数。您可以指定预条件子矩阵 M 或其矩阵因子 M = M1*M2 来改进线性方程组的数值方面，使 tfqmr 更容易快速收敛。
tfqmr 将未指定的预条件子视为单位矩阵。
将 M 指定为函数句柄
您可以选择将 M、M1 或 M2 中的任一个指定为函数句柄而不是矩阵。函数句柄执行矩阵向量运算，而不是构建整个预条件子矩阵，从而使计算更加高效。
要使用函数句柄，请使用函数签名 function mfun(x)。参数化函数说明如何在必要时为函数 mfun 提供附加参数。函数调用 mfun(x) 必须返回 M\x 或 M2(M1)\x 的值。
数据类型： Float64 | Function
复数支持： 是
x0 - 初始估计值由零组成的列向量（默认） | 列向量
初始估计值，指定为长度等于 size(A,2) 的列向量。如果您能为 tfqmr 提供比默认的零向量更合理的初始估计值 x0，则它可以节省计算时间并帮助算法更快地收敛。
数据类型： Float64
复数支持： 是
verbose - 信息显示指示器true（默认） | false
信息显示指示器，指定为逻辑值。如果该变量值为 true（默认值），则 tfqmr 会显示收敛信息：如果尝试成功，tfqmr 会显示一条消息来确认收敛。如果 tfqmr 无法在达到最大迭代次数后收敛或出于任何原因暂停，则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。否则不会显示。
# 输出参数
x - 线性方程组的解列向量
线性方程组的解，以列向量形式返回。该输出给出线性方程组 A*x = b 的近似解。如果计算成功 (flag = 0)，则 relres 小于或等于 tol。
每当计算不成功 (flag != 0) 时，tfqmr 返回的解 x 是在所有迭代中计算出的残差范数最小的解。
flag - 收敛标志标量
收敛标志，返回下表中的标量值之一。收敛标志指示计算是否成功，并区分几种不同形式的失败。
| 标志值 | 收敛 |
|---|---|
| 0 | 成功 - tfqmr 在 maxit 次迭代内收敛至所需容差 tol。 |
| 1 | 失败 - tfqmr 执行了 maxit 次迭代，但未收敛。 |
| 2 | 失败 - 预条件子矩阵 M 或 M = M1*M2 为病态。 |
| 3 | 失败 - tfqmr 在经过两次相同的连续迭代后已停滞。 |
| 4 | 失败 - 由 tfqmr 算法计算的标量数量之一变得太小或太大，无法继续计算。 |
relres - 相对残差标量
相对残差，以标量形式返回。相对残差 relres = norm(b-A*x)/norm(b) 指示解的准确度。如果计算在 maxit 次迭代内收敛于容差 tol，则 relres <= tol。
数据类型： Float64
iter - 迭代编号标量
迭代编号，以标量形式返回。此输出指示计算出 x 的解时所用的迭代次数。tfqmr 的每个外迭代包括两个内迭代，因此 iter 可以以十进制迭代次数形式返回。
数据类型： Int64
resvec - 残差向量
残差，以向量形式返回。残差 norm(b-A*x) 揭示对于给定的 x 值，算法接近收敛的程度。resvec 中元素的数量等于半迭代次数。您可以检查 resvec 的内容，以帮助决定是否更改 tol 或 maxit 的值。
数据类型： Float64
# 详细信息
无转置最小残差法
开发 TFQMR 方法是为了避免使用转置。与“未平方”的方法相比，这些“平方”方法每步都需要一个额外的矩阵向量积，因此效率稍低。
TFQMR 方法具有平滑的收敛性。尽管如此，由于 TFQMR 最终是使用 BiCG 多项式，因此，它也会有一些故障[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.