大型稀疏线性系统的求解，重排序和预条件更差？

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

在当前的方法中，很有可能您没有做错任何事情（至少，我没有发现明显的bug）

几点注意：

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

如果不知道您的系统来自何处，则很难说出任何信息。然而：

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