大型稀疏线性系统的求解，使用重排序和预处理更糟

Question

我有一个线性系统，我想求解一个60000x60000矩阵，其中有大约6000000个非零条目。

我目前的方法是用逆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

在当前方法中很有可能没有做错任何事情(至少，我没有发现明显的错误(。

注意事项：

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

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

1. 检查矩阵的小版本的条件号(如果它在某种程度上代表了您的原始问题(。高条件数表示矩阵的条件存在问题；从而导致迭代解(或极端情况下的任何类型的解(的潜在收敛性差或收敛性差
2. 共轭梯度适用于对称和正定的系统。它可以在其他情况下收敛；然而，它可能会失败的条件良好的问题，不是肯定的