Skip to content

Instantly share code, notes, and snippets.

@denis-bz
Last active June 11, 2024 15:27
Show Gist options
  • Save denis-bz/502a22ee092f0758ee11c9b731828602 to your computer and use it in GitHub Desktop.
Save denis-bz/502a22ee092f0758ee11c9b731828602 to your computer and use it in GitHub Desktop.
Poisson2 - 4 I: nice test matrices for linear solvers and eigensolvers 25 Jun 2022

Poisson2 - 4 I: nice test matrices for linear solvers and eigensolvers

Purpose: a simple generator of test matrices for linear solvers and eigensolvers

Keywords: sparse-matrix, linear-solver, eigensolver, arpack, nullspace, python

poisson2( n ) below generates sparse $n^2 \ x \ n^2$ matrices $A$ with $A$ symmetric (aka Hermitian) and positive-definite; solving $A \ x = b$ is then relatively easy. $A_4 = A - 4 I$ is more interesting: there are $n$ vectors $z_i$ in the hyperplane $A_4 \ z_i = 0 $ to within floating-point eps, 2e-16. If $A_4 \ x = b$, $A_4 \ (x + z_i) = b$ too; how should a solver pick one $x$ in the whole plane ? This is an old problem, discussed in thousands of papers, programs and books. (If anyone knows of an introduction / overview for practitioners, please comment.)

Here's a plot of 10 nullspace vectors of $poisson2( 10 ) - 4 \ I$ (there are more, any rotation of these 10 also span the nullspace):

22jun-plot-nullspace

Which $b$ make $A \ x = b$ hard to solve ?

How hard it is to solve $A x = b$ depends of course on $b$. For poisson2, $\ b = $ a point source, 1 in the middle and the rest 0, is tough. $(poisson1 - 2I)^{-1}$ shows why this is so:

tri101-inv

Finding eigenvalues near 0: shift-invert

To find eigenvalues near 0, scipy.sparse.linalg.eigsh( A, ... sigma=0 ... ) calls Arpack in "shift-invert" mode, which solves $A \ x = b$ many times with various $b$. $A$ cannot have eigenvalues exactly 0, or 0 within floating-point precision. $poisson2 - 4 \ I$ is singular, so that won't work -- shift-invert will (should) then raise an error and quit. Instead, nudge it a little away from 0: $poisson2 - 4.0001 I$ is non-singular, with a cloud of $n$ eigenvalues around -0.0001. Such clouds can break eigensolvers.

Arpack is (too me) a very black box, with labyrinths of parameters, solvers and solver options. Scipy Arpack uses

  • a direct solver splu if $A$ is a matrix (fast if umfpack is installed), or
  • an iterative solver gmres if $A$ is a LinearOperator (very slow -- Lgmres is much faster, at least for $A$ non-posdef).

See also

scicomp.stack arpack or gmres

Dhillon, "Current inverse iteration software can fail", 1998 ff.

"Invop() inverse operator for scipy.sparse eigenvalues with ARPACK" under my gists, to track and plot convergence of Arpack with various solvers

Comments welcome

cheers
-- denis 25 June 2022

# from: test-lgmres.py test="poisson2" n=100 incdiag=-4.0001 solver=lgmres tol=1e-5 maxiter=0 restart=500 innerm=30 printevery=100 save="5jul-lgmres-poisson2-n100-inc-4.0001-tol1e-5.npz"
# 5 Jul 2022 10:58z np/eig/gmres Denis 10.14.6
▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄
testmatrix poisson2: A (10000, 10000) incdiag -4.0001
A.diag: [-0.0001 -0.0001 -0.0001 -0.0001 -0.0001 ... -0.0001 -0.0001 -0.0001 -0.0001 -0.0001]
{ test-lgmres.py lgmres poisson2 N 10000 incdiag -4.0001 tol 1e-05 atol 1e-05 maxiter 0 maxtime 300 restart 500 innerm 30 printevery 100
solvex0 |b - A x0|: 1
# b-Ax || dx x
callback 0 0s res 1 dx 0 .. 0 x 0 .. 0
callback 100 1s res 0.04 dx -1.3e+02 .. 1.3e+02 x 44 .. 132
callback 200 3s res 0.0073 dx -16 .. 16 x 49.2 .. 147
callback 300 4s res 0.00043 dx -0.54 .. 0.54 x 49.3 .. 148
callback 400 6s res 0.00023 dx -0.37 .. 0.37 x 49.5 .. 148
lgmres: iter 485 maxiter 10000 |b-Ax| 1e-05 x -1.5e+02 .. 1.5e+02 A∙1 (10000, 10000) -4 .. -2 b 0 .. 1 tol 1e-05 atol 1e-05 inner_m 30
time: 7.0 sec 1.0 cpus 0 iter lgmres niter 0 solvex0 0
x: || 1.22e+03 min -148.5 at (50, 50) max 148.5 at (49, 49)
quantiles: [-149 -0.11 -0.0099 1.7e-05 0.0099 0.11 149]
middle:
[[ 66 . . . . . . . 33 .]
[ . -66 . . . . . -33 . -33]
[ . . 66 . . . 33 . 33 .]
[ . . . -66 . -33 . -33 . .]
[ . . . . 99 . 33 . . .]
[ . . . -33 . -99 . . . .]
[ . . 33 . 33 . 66 . . .]
[ . -33 . -33 . . . -66 . .]
[ 33 . 33 . . . . . 66 .]
[ . -33 . . . . . . . -66]]
res: || 1e-05 min -9.139e-07 at (49, 49) max 9.015e-07 at (50, 50)
quantiles: [-9.1e-07 -9.4e-08 -4.5e-08 -3e-10 4.5e-08 9.3e-08 9e-07]
middle:
[[ -59 2 -2 3 -5 6 -5 13 -34 13]
[ 2 69 -3 6 -4 5 -7 36 -13 33]
[ -2 -3 -73 3 -6 4 -35 7 -34 13]
[ 3 6 3 72 -3 34 -5 33 -7 .]
[ -5 -4 -6 -3 -100 3 -32 4 1 6]
[ 6 5 4 34 3 99 -4 -2 -3 1]
[ -5 -7 -35 -5 -32 -4 -65 4 . 3]
[ 13 36 7 33 4 -2 4 66 -3 5]
[ -34 -13 -34 -7 1 -3 . -3 -72 3]
[ 13 33 13 . 6 1 3 5 3 77]]
params: test-lgmres.py lgmres poisson2 N 10000 incdiag -4.0001 tol 1e-05 atol 1e-05 maxiter 0 maxtime 300 restart 500 innerm 30 printevery 100
}
> 5jul-lgmres-poisson2-n100-inc-4.0001-tol1e-5.npz
9.86 real time 9.69 user 0.78 sys 106% of 4 cores imac 2.7 GHz
#!/usr/bin/env python3
""" poisson_nd( n, d ) -> sparse matrix `A` 1d n x n, 2d n^2 x n^2, 3d n^3 x n^3
poisson2( 3 ):
[[ 4 -1 . -1 . . . . .]
[-1 4 -1 . -1 . . . .]
[ . -1 4 . . -1 . . .]
[-1 . . 4 -1 . -1 . .]
[ . -1 . -1 4 -1 . -1 .]
[ . . -1 . -1 4 . . -1]
[ . . . -1 . . 4 -1 .]
[ . . . . -1 . -1 4 -1]
[ . . . . . -1 . -1 4]]
(A_h = this / h^2)
poisson2( n ) - 4 I has n null vecs,
poisson2( n ) - 4.0001 I has n near-null vecs -- a nice test for arpack and *gmres
inv( poisson1( n ) - 2 I ) is all 1 0 -1 ! plot under my gists
"""
# see also Invop, Cheblinop
from functools import partial
import numpy as np
from scipy import sparse
#...............................................................................
def poisson1( n ) -> "sparse csc tridiagonal [-1 2 1] ...":
# Make Gilbert Strang's favorite matrix
# http://www-math.mit.edu/~gs/PIX/cupcakematrix.jpg
# w Discrete_Laplace_operator
P1d = sparse.diags( [[-1]*(n-1), [2]*n, [-1]*(n-1)],
[-1, 0, 1]) # ~ - f''
return P1d
def poisson_nd( n, d=2 ):
""" -> A sparse csc 1d n x n, 2d n^2 x n^2, 3d n^3 x n^3 """
# w Kronecker_sum_of_discrete_Laplacians
# https://math.stackexchange.com/questions/2629649/understanding-high-and-low-frequency-eigenmodes
# poisson2 evals: 2 (pi / n)^2 .. 8, cond ~ 0.8 n^2
n = int(n)
assert n > 0, n
P1d = poisson1(n)
if d == 1:
return P1d.tocsc()
return sparse.kronsum( P1d, poisson_nd( n, d - 1 )).tocsc()
# kron I A + kron B I
poisson2 = partial( poisson_nd, d=2 )
poisson3 = partial( poisson_nd, d=3 )
poisson2.__name__ = "poisson2"
poisson3.__name__ = "poisson3"
__version__ = "2022-06-24 June" # denis-bz-py t-online.de
#...............................................................................
if __name__ == "__main__": # run various solvers, direct / iterative
import sys
from time import time
from numpy.linalg import norm
from scipy.sparse import linalg as sparselin
print( 80 * "▄" )
# print( "▄ bp" ); breakpoint()
try:
from scikits import umfpack
print( "umfpack:", umfpack.__version__ )
except ImportError:
print( "umfpack:", None )
np.set_printoptions( edgeitems=10, threshold=10, linewidth=120, formatter = dict(
float = lambda x: "%.2g" % x )) # float arrays %.2g
def sparsesolver( name, A ):
# https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems
if name == "dense":
return lambda b: np.linalg.solve( A.A, b )
elif name == "splu": # direct, with umfpack v fast
return sparselin.splu(A).solve
elif name == "spsolve": # direct, v slow
return lambda b: sparselin.spsolve( A, b )
else:
solver = getattr( sparselin, name ) # iter lgmres ... sparsesolvers.py tl;dr
return lambda b: solver( A, b, atol=1e-5 )
def nbytes( A ):
return A.nbytes if hasattr( A, "nbytes" ) \
else A.data.nbytes + A.indices.nbytes + A.indptr.nbytes
#...........................................................................
solvers = "splu lgmres gmres" # spsolve"
# solvers = "dense"
d = 2
n = 100
incdiag = -4.0001 # -4: n null vecs of the n^2
save = 0
# to change these params, run this.py a=1 b=None 'c = expr' ... in sh or ipython --
for arg in sys.argv[1:]:
exec( arg )
#...........................................................................
A = poisson_nd( n, d=d ) # n^d x n^d
N = A.shape[0]
if incdiag != 0:
A += incdiag * sparse.eye( N )
params = "poisson%d n %g incdiag %g A %s %.0f mbytes " % (
d, n, incdiag, A.shape, nbytes(A) / 1e6 )
print( "params:", params )
# print( "A: \n", A[:10,:10].A )
b = np.zeros( d * (n,) )
b[ d * (n // 2,) ] = 1 # midpoint
b = b.reshape( -1 )
x0 = np.zeros( N ) # ?
for solvername in solvers.split():
print( "\n--", solvername )
t0 = time() # wall clock, e.g. 329% of 4 cores
#.......................................................................
solver = sparsesolver( solvername, A ) # , b, x0
x = solver( b )
if isinstance( x, tuple ):
x, info = x # iterative solvers
b_Ax = norm( b - A.dot(x) ) / norm( b )
print(
"time: %6.1f sec %-8s poisson%d n %d incdiag %g |b-Ax|/|b| %.2g x %g .. %g " % (
time() - t0, solvername, d, n, incdiag,
b_Ax, x.min(), x.max() ))
if d == 1:
print( x )
elif d == 2 and n <= 20:
print( x.reshape( n, n )) # incdiag: ripple
if save:
out = "poisson%d-n%d-%s.npz" % (d, n, solvername)
print( "\n>", out, "\n" )
np.savez( out, x=x, n=n, params=params )
# $sslin/_dsolve/linsolve.py with umfpack --
# time: 1.0 sec splu poisson2 n 300 incdiag -4.0001
# time: 30.8 sec spsolve poisson2 n 300 incdiag -4.0001
# time: 6.7 sec splu poisson3 n 30 incdiag -6.0001
# time: 47.2 sec spsolve poisson3 n 30 incdiag -6.0001
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment