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

Question

I have a linear system with a 60000x60000 matrix that I wish to solve, with about 6,000,000 nonzero entries in it.

My current approach is to reorder the matrix with reverse cuthill mckee, factorize the matrix, and then solve it with preconditioned conjugate gradient, but I'm not getting very good results and I don't understand why. The reordering looks reasonable.

Below I have attached a simple example where I only use a subsystem of the matrix I'm trying to solve.

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))

The results I get back from this are:

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

So not only does the reordering and factorization take a very look time, but the solving with cg doesn't go any faster and doesn't give anymore of a correct solution, in fact the residual is way worse!

Can anyone tell me what exactly I'm doing wrong here?

Anton Menshov Anton Menshov · Accepted Answer · 2020-05-04T11:49:35

There is a high chance that you are doing nothing wrong within your current approach (at least, I was not able to spot an obvious bug).

A couple of notes:

Residual of 29.10655954230801 and 2.5236861383747353 after 500 iterations are effectively the same: your iterative solution has not converged.
You seem to request a very high iterative solver tolerance of 1E-12. That would not matter here, as you have a problem that does not converge at all.
Factorization (of ILU) should take approximately this time. I am not surprised to see this number for such a system. Not so familiar with this implementation of Cuthill-McKee.

Without knowing where your system comes from, it would be very hard to say anything. However:

Check the condition number for your small version of the matrix (if it is somewhat representative of your original problem). High condition number would indicate a problem with the conditioning of the matrix; thus, potential poor convergence or ill-convergence of the iterative solution (or any type of solution in the extreme case).
Conjugate gradient is intended for systems that are symmetric and positive-definite. It can converge for other cases; however, it can fail for well-conditioned problems that are not positive-definite.

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

1 Answers