大型稀疏线性系统求解，重新排序和预处理器更糟？

Question

我有一个线性系统，我想求解一个 60000x60000 矩阵，其中大约有 6,000,000 个非零项。

我目前的做法是用reverse cuthill mckee对矩阵重新排序，对矩阵进行因式分解，然后用预条件共轭梯度求解，但我没有得到很好的结果，我不明白为什么。重新排序看起来很合理。

下面我附上了一个简单的例子，我只使用了我试图求解的矩阵的一个子系统。

import matplotlib
matplotlib.use('TkAgg') #TkAgg for vizual
from matplotlib import pyplot as plt
import time
import numpy as np
import scipy
from scipy import sparse
from scipy.sparse.linalg import LinearOperator, spilu, cg
from numpy.linalg import norm
L = sparse.load_npz("L_Matrix.npz")
n = 20000
b = np.random.randn((n))
L2 = L[0:n,0:n].copy()

t00 = time.time()
perm = scipy.sparse.csgraph.reverse_cuthill_mckee(L2, symmetric_mode=True)
I,J = np.ix_(perm,perm)
bp = b[perm]
L2p = L2[I, J]
t01 = time.time()

fig = plt.figure(0, figsize=[20, 10])
plt.subplot(1, 2, 1)
plt.spy(L2)
plt.subplot(1, 2, 2)
plt.spy(L2p)
plt.pause(1)
# plt.pause(1)

t0 = time.time()
print("reordering took {}".format(t0-t00))
ilu = spilu(L2p)
t1 = time.time()
print("Factorization took {}".format(t1-t0))

Mx = lambda x: ilu.solve(x)
M = LinearOperator((n, n), Mx)
x,stat = cg(L2p, bp, tol=1e-12, maxiter=500, M=M)
t2 = time.time()
print("pcg took {} s, and had status {}".format(t2-t1,stat))
print("reorder+pcg+factor = {} s".format(t2-t00))
bsol = L2p @ x
R = norm(bsol - bp)
print("pcg residual = {}".format(R))

x,stat = cg(L2, b, tol=1e-12, maxiter=500)
t3 = time.time()
print("cg took {} s, and had status {}".format(t3-t2,stat))
bsol = L2 @ x
R = norm(bsol - b)
print("pcg residual = {}".format(R))

我得到的结果是：

reordering took 66.32699060440063
Factorization took 64.96741151809692
pcg took 12.732918739318848 s, and had status 500
reorder+pcg+factor = 144.0273208618164 s
pcg residual = 29.10655954230801
cg took 1.2132720947265625 s, and had status 500
pcg residual = 2.5236861383747353

因此，不仅重新排序和因式分解需要很长时间，而且使用 cg 求解并没有更快，也没有给出正确的解决方案，事实上，残差更糟糕！

谁能告诉我我到底做错了什么？

Answer 1

您当前的方法很可能没有做错（至少，我没能发现明显的错误）。

一些注意事项：

29.10655954230801 和 2.5236861383747353 的残差在 500 次迭代后实际上是相同的：您的迭代解没有收敛。
您似乎要求 1E-12 的非常高的迭代求解器容差。这在这里无关紧要，因为你有一个根本不收敛的问题。
（ILU 的）因式分解大约需要这个时间。对于这样一个系统，看到这个数字我并不感到惊讶。不太熟悉 Cuthill-McKee 的这个实现。

如果不知道你的系统来自哪里，就很难说什么。然而：

检查您的小版本矩阵的条件编号（如果它在某种程度上代表您的原始问题）。高条件数表示矩阵条件有问题；因此，迭代解决方案（或极端情况下的任何类型的解决方案）可能收敛不良或收敛不良。
共轭梯度适用于 symmetric and positive-definite 的系统。它可以收敛于其他情况；但是，对于非正定的条件良好的问题，它可能会失败。

大型稀疏线性系统求解，重新排序和预处理器更糟？

large sparse linear system solving, worse with reordering and preconditioner?

python

numpy

linear-algebra

scipy

sparse-matrix