大型稀疏线性系统求解，用重排序和前置条件器更糟糕？

Question

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

我目前的方法是用反向切希尔麦基对矩阵进行重排序，对矩阵进行因子化，然后用预设共轭梯度来求解，但我没有得到很好的结果，我不明白为什么。"重排序 "看起来很合理。

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

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. 共轭梯度是针对以下系统的。对称 和 正定. 对于其他情况，它可以收敛；但对于非正定型的条件良好的问题，它可能失败。