denis-bz/0-Shift-invert-in-pictures.md

## 0-Shift-invert-in-pictures.md

      
    Raw
  

              0-Shift-invert-in-pictures.md
            
          
    How shift-invert finds eigenvalues, in pictures

Keywords: eigenvalues, eigenvectors, shift-invert, ARPACK, scipy, python
Pictures first:

Background

Eigenvalues and eigenvectors
describe a matrix or linear transform A of N-space
in terms of INDEPENDENT (orthogonal, decoupled) vectors V0 V1 ...
which A stretches by "stretch factors" λ0 λ1 ...
A V0 = λ0 V0
A V1 = λ1 V1  # independent, orthogonal to, V0
A V2 = λ2 V2  # independent of V0 and V1 ...
...

Basics of shifted eigenvalues


if A has eigenvalues λ0 λ1 ..., A - shift I has eigenvalues λ0-shift, λ1-shift ...
if A has eigenvalues λ0 λ1 ... all != 0,
the inverse matrix Ainv or inv(A) has eigenvalues 1/λ0 1/λ1 ...

Exercise: what are the eigenvalues of inv( A - shift I ) ?
How are its eigenvectors related to those of A ?
Hint: first do the case A diagonal, e.g. x, y, z  ->  2x, 3y, -z.
In some problem areas, the largest eigenvalues are of interest; in others,
the smallest (for complex eigenvalues, the smallest in magnitude).
In practice, finding large eigenvalues is easier than finding small ones.
The raison d'etre of the shift-invert method is to transform
"small, difficult" to "large, easier", and back again.
Some kinds of matrices


diagonal, x -> [A00*x0  A11*x1 ...] per component.
To understand a matrix algorithm, start with A diagonal.
A real and symmetric, Aij == Aji => λ0 λ1 ... are all real
A positive definite, posdef: real symmetric, and λ0 λ1 ... all > 0
A non-negative-definite: λ0 λ1 ... all >= 0
(Beware: the eigenvalues of matrices X' X are in theory all real >= 0 ,
but in practice you may see a few tiny < 0, or tiny imaginary parts.)
Singular, some eigenvalues very near 0.
Shift-invert for the smallest eigenvalues, shift 0, approaches 1 / 0 -- trouble.
non-Normal matrices
have no orthonormal basis V0 V1 ... at all.

Which matrices are easy, which hard, for iterative eigensolvers like Arpack ?
There's no general answer.
Posdef matrices like X' X are far easier than e.g. sparse random
Circular law
matrices.
Closely spaced eigenvalues, or large null space,
cause slow convegence and inaccurate eigenpairs.
Try scipy.linalg.eigvals first

If your N x N dense matrix fits in memory, try scipy.linalg.eigvalsh / eigh overnight,
first with .astype(np.float32).
For example, on a 29027 x 29027 almost-non-negative-definite matrix from
SO finding-smallest-eigenvectors-of-large-sparse-matrix
, np.linalg.eigh found all eigenvalues and eigenvectors
of the 7 Gb dense matrix in under an hour on my old 4-core, 16 Gb imac.
Messing with Arpack tol etc. for just 10 eigenpairs, of uncertain accuracy at that,
took a lot longer in human time.
Added 2 July:
I see speed ratios ~ 1 : 2 : 5 : 5 for scipy eigvalsh float32, scipy eigvalsh float64,
numpy eigvalsh 32 or 64, on randomsparse().A, dense.
see Test-eigvals-numpy-scipy under
gists.
For dense symmetric matrices, see also
scipy.linalg.eigvalsh
eigvals=[0,9] for the 10 smallest,
eigvals=[n-10,n-1] for the 10 largest eigenpairs.
(Does anyone know of a free cloud computer with numpy, scipy and lots of memory ?)
Scipy Arpack shift-invert

Background:
scipy.sparse.linalg.eigs
.
Shift-invert mode for eigenvalues near 0 (sigma=0) or near a given shift,
must solve the sparse linear system (A - shift I) x = b many many times.
Obviously, slow linear solver => slow Arpack.
But the runtimes of various solvers even on random A
can vary by a factor of 100, depending on where shift puts the eigenvalues --
see Circular law
and gists,
Arpack's default solvers are splu for matrices
and gmres (with curious tol) for LinearOperators.
I'd suggest:


for symmetric matrices with eigenvalues >= 0,
use sksparse.cholmod
with beta=1e-6 or so to nudge near-0 eigenvalues positive.


for scipy.sparse matrices,
use spsolve with scikit-umfpack;
I've seen a 50 times speedup.
(If pip install scikit-umfpack doesn't work,
see the setup.py under
gists,
or try conda.)


for LinearOperators, try
lgmres
.


start with tol=1e-5, not the default tol=0, machine precision,
and a sensible atol.
(IFAQ: what is tol in Arpack ?)


start with v0=np.ones(n), not the default random start vector
(unless you're planning an article in the Journal of Irreproducible Results).
A start vector from a half-size problem, or from previous runs, must work much better.


How to run Arpack with a linear-solver-of-your-choice

Todo: writeup class Invop
Check, check again

Always check the accuracy of calculated eigenvalues and vectors:
evals, V = np.linalg.eig( A )  # or eigh
	# or sparse: evals, V = eigs( A, k=10, tol=1e-6 )  # small: sigma=0
	# evals, V = eigsort( evals, V )
diffs = A.dot( V ) - V * evals  # sparse: N x 10
maxdiffs = np.linalg.norm( diffs, axis=0, ord=np.inf )
print( "|Av - λv|_max:", maxdiffs )

See also

Scipy Arpack tutorial

Eigenvalue algorithm

Numerical Recipes pages 563-599

Scicomp CG GMRES

Test-sparse-solvers under my gists
Comments welcome, test cases welcome

cheers

-- denis-bz-py t-online.de  2 July 2020