Inverse operators for scipy.sparse eigs, arpack 15 Nov 2020

0-invop.md

Invop() inverse operator for scipy.sparse eigenvalues with ARPACK

Problem: find a few eigenvalues λ and eigenvectors v of
A v = λ v, \ A a scipy.sparse matrix, largest or smallest |λ|
or K v = λ M v with stiffness matrix K and mass matrix M.

Purpose: experiment with different linear solvers, with verbose to track solver calls.
Keywords: scipy.sparse, linear-solvers, eigenvalue, generalized-eigenvalue, Arpack, logging

In the following examples, eigs( A [M] [Minv] [OPinv] [sigma] ... )
is short for evals, V = eigs or eigsh( A ... ) with parameters e.g. k=10, tol=1e-5, v0=np.ones(n) ), see below.

Examples, the largest eigenvalues and eigenvectors of A v = λ M v:

from scipy.sparse.linalg import eigs, eigsh
from invop import Invop, Aminus
eigs( A, M=M ... )  # or eigsh
eigs( A, M=M, Minv=Invop( M, solver="splu" ) ... )      # the default for M a matrix
eigs( A, M=M, Minv=Invop( M, solver="cholesky" ) ... )  # fast for M posdef
eigs( A, M=M, Minv=Invop( M, solver="gmres" ) ... )     # the default for M a LinOp
eigs( A, M=M, Minv=Invop( M, solver="lgmres" ) ... )    # see below

Shift-invert mode

eigs( ... sigma=0 ... ) tells Arpack to find eigenvalues near 0, in the so-called "shift-invert" mode. For example:

Ainv = Invop( A, solver="splu" )  # b -> b \ A
eigs( A, OPinv=Ainv, sigma=0, [M=M] ... )

For this mode, Arpack needs to solve A x = b many many times. The default Invop( ... verbose=2 ) prints a line per solver call; see the example logfiles below.

(In fact eigs( A, OPinv=Ainv, sigma=0 ... ) doesn't use A at all, except to check its dtype.)

Notes on some solvers

There are dozens of linear solvers. Their runtimes (multicore), memory, understandability, can vary a lot, depending on A's size, spectrum, sparsity pattern ... One size cannot fit all. The default solvers for scipy Arpack shift-invert mode, splu for sparse matrices and GMRES for LinearOperators, are based on experience which may or may not match your problem.

If the dense matrix A.A fits in memory, running scipy.linalg.eig or eigh( A.A ), even overnight, can save you hours of messing about with ARPACK parameters and different solvers. Sampling the spectrum with subset_by_index or subset_by_value can possibly give some insight into your problem.

It's far easier to find eigenpairs of symmetric or Hermitian matrices, which have real spectrum, than of general matrices with complex spectrum. For posdef matrices, sksparse.cholmod.cholesky is very fast. Use cholinc=1e-10 or so to nudge almost-posdef away from 0 -- X'X can easily have tiny eigenvalues < 0.

Iterative solvers like GMRES are much slower than direct solvers, but direct solvers splu etc. will need more memory. Wikipedia GMRES says "fast convergence when the eigenvalues of A are clustered away from the origin". When this is not the case, I find that LGMRES can be ~ 10 x faster than GMRES, e.g. on poisson3 - 5.9 I; see the example logfiles below.

Notes on ARPACK parameters

The default tol=0, machine precision, can be very slow. I'd start with tol=1e-6 and see if eigcheck A V - V Λ looks reasonable. (For Linops, arpack.py gmres_loose is not loose at all, e.g. 1000 * np.sqrt(b.size 1e6) * eps 2e-16 = 2e-10 .)

The default start vector v0 random is good for the Journal of Irreproducible Results. A BIG and easy improvement is to use v0 from previous runs, or from scaling up a half-size problem.

Lgmres can run much faster with a preconditioner, a fast cheap approximate solver. "Preconditioning is an art, a science, an industry."

The two-matrix aka generalized eigenvalue problem K v = λ M v

This is a common problem in the finite-element method for analyzing structures such as bridges: find eigenvalues λ near 0 for K v = λ M v, with K a positive-definite stiffness matrix and M a positive-definite mass matrix.

eigsh( K, M=M, OPinv=Invop( K, solver=... ), sigma=0, ... )

Similarly, for λ near sigma=2:

inv_A_minus_sigmaM = Invop( Aminus( A, 2, M ), solver=... )
eigsh( K, M=M, OPinv=inv_A_minus_sigmaM, sigma=2 ... )

There are many papers on this problem, perhaps more than real test cases :(

(By the way, the generalized eigenvalue problem and generalized eigenvector problem are different altogether -- the latter is (A - λI)^k x = 0.)

Example logs, splu vs. spsolve

Invop() with verbose=2 (the default) can help understand what a solver is doing. For example, solve = scipy.linalg.splu(A).solve first builds a representation of A-1; this can be slow and can explode memory 💥 , depending on the sparsity structure of A^-1 and on the options. Arpack calls of solve(b) are then pretty fast:

A = testmatrix( "poisson3", n=40, incdiag=-5.9  )  # 40^3, see testmatrix.py
Ainv = Invop( A, solver="splu" )  # default verbose=2 
eigsh( A, sigma=0, OPinv=Ainv ... )

Invop splu  A: (64000, 64000) csc_matrix  6.85 nnz/row ...

# solver  seconds wall clock   |b-Ax|/|b|   dx   x   [lgmres: nr outer iter]
invop splu 142s  |b-Ax|/|b| 1e-13     dx -0.0039 .. 0.0022   x -0.0039 .. 0.0022   
...
invop splu 148s  |b-Ax|/|b| 1e-10     dx -25 .. 25   x -12 .. 11   
eigsh: 147.8 480.9 sec  poisson3   n 64000  58 calls of splu

Spsolve with umfpack is slower per call, but ~ linear in the number of calls:

# solver  seconds wall clock   |b-Ax|/|b|   dx   x   [lgmres: nr outer iter]
invop spsolve   3s  |b-Ax|/|b| 1e-16     dx -0.0039 .. 0.0022    x -0.0039 .. 0.0022   
invop spsolve   7s  |b-Ax|/|b| 3.8e-14   dx -3 .. 4.4    x -3 .. 4.4
...
invop spsolve 214s  |b-Ax|/|b| 9.4e-14   dx -18 .. 19    x -6.7 .. 6.4   
eigsh: 214.2 786.2 sec  poisson3   n 64000  64 calls of spsolve

Another example, GMRES vs. Lgmres:

testmatrix poisson2: A (2500, 2500)  4.92 nnz/row  -1 .. 0.1  incdiag -3.9  max |A - A'|: 0

GMRES:
eigsh: 1755.9 sec  30 calls of gmres, 30 stopped at maxiter 25000
evals [-0.00065 -0.00068 -0.0011 -0.0013 -0.0029 0.004]
eigcheck: |Av - λv|  max 0.00089  [0.0003 4e-05 8e-05 0.0009 4e-05 9e-06]

Lgmres:
eigsh: 107.3 sec  30 calls of lgmres
evals [-0.00054 -0.00054 -0.001 -0.001 -0.0029 0.004]
eigcheck: |Av - λv|  max 5e-07  [1e-07 2e-07 4e-07 2e-07 3e-08 5e-07]

See also

shift-inverse-in-pictures
http://www.netlib.org/lapack/explore-html/modules.html
https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html
https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html
https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html \

Dhillon, "Current inverse iteration software can fail", 1998 \

https://scicomp.stackexchange.com/questions/tagged/linear-solver 300
https://scicomp.stackexchange.com/questions/tagged/gmres 40 \

cheers
-- denis 15 Nov 2020

1-invop-files.md

Invop files outline

-- 108 invop.py
class Invop( object ):  """ Invop( A, M=None, solver="splu" ... ): a function b -> b \ A 
	def __init__( s, 
	def __call__( s, b ):	""" b -> x = b \ A  verbose """ 
	def shape( s ): 
	def dtype( s ): 
def Aminus( A, c, M=None, inplace=False ):	""" A - c I or A - c M  for Invop( Aminus( A, c )), shift-invert """ 


-- 112 test-mikotapair.py
def KM_smalleigs( K, M, sigma=0, k=10, tol=1e-6, **kw ):	""" λ nearest sigma, shift-invert: 
def KM_bigeigs( K, M, k=10, tol=1e-6, **kw ):	""" large λ, K x = λ M x  ->  Minv K x = λ x """ 


--  94 poisson_nd.py
def poisson1(n): 
def poisson_nd( n, d=2 ):	""" poisson 1d / 2d n^2, n^2 / 3d n^3, n^3 ... """ 


--  78 sparserotate.py
def sparserotate( eigvals, p=2 ):	""" :param eigvals: 1d arraylike, len n a multiple of 2^p 
def hadamard( p=2 ):	""" -> H 2^p x 2^p +- 1 


-- 155 sparsesolvers.py
class Sparsesolver( object ):		""" Ainv = Sparsesolver( "splu", A ... ) for arpack eigs( OPinv or Minv ) """
	def __call__( s, b, **kw ):		""" solve(b), verbose: print solver time b-Ax xmin xmax [per iter] """ 


-- 186 testmatrix.py
def testmatrix( name, n=0, incdiag=0, testsymm=True, **kw ):	""" load or generate a sparse csc matrix 
def transforms( trans, A ):		""" :param trans: "flip Ltri LL rotate float32" 
def _getmatrix( name, n=0, **kw ):	""" generate a matrix, poisson* / powerlaw / randomsparse ... 
def try_npz_mtx( dir ): 
def load_npz_mtx( filename ): 
def curve( n,	""" -> spectrum like scikit-learn 20newsgroups ndoc=4000 A A'


-- run-invop.sh

Comments welcome, test cases most welcome

invop.py

#!/usr/bin/env python3
""" please see invop.md """

import sys
import time
import numpy as np
from scipy import sparse
        # $scipy/sparse/linalg/dsolve/linsolve.py
        # $scipy/sparse/linalg/isolve/iterative.py
from scipy.sparse.linalg import LinearOperator, aslinearoperator
        # $scipy/sparse/linalg/interface.py

from sparsesolvers import Sparsesolver  # ( "splu ...", A, tol=1e-5, x0=1 )
from zutil import arraydesc, lohi

__version__ = "2020-10-02 Oct"

#...............................................................................
class Invop( object ):
    """ Invop( A, M=None, solver="spsolve" ... ): a function b -> b \ A
        wraps Sparsesolver
        verbose=1 or 2 tracks iterations
        please see invop.md
    """
    # (object) not (LinearOperator) -- arpack then uses gmres

    def __init__( s,
            A,              # sparse
            M=None,         # b -> A^-1 (M b)
            solver="splu",  # see *-{spsolve,splu}*.log
            tol=1e-5,
            atol=None,
            x0=1,
            verbose=2,
            # maxtime=600,
            **solverkw ):
        s.A = A
        s.solvername = solver
        s.M = M  # b -> (M b) \ A
        s.verbose = verbose
        s.t0 = time.time()
        s.xprev = 0

        #.......................................................................
        s.solver = Sparsesolver( s.solvername, A=A, tol=tol, atol=atol,  # b -> b \ A
                    x0=x0, **solverkw )
        if verbose:
            print( "\nInvop %s  A: %s  M: %s  tol %g  maxiter %d " % (
                    s.solvername, arraydesc( A ), arraydesc( M ),
                    tol, s.solver.maxiter ))
            if verbose >= 2:
                print( "# solver  seconds wall clock   |b-Ax|/|b|   dx   x   [lgmres: nr outer iter]" )

    #...........................................................................
    def __call__( s, b ):
        """ b -> x = b \ A  verbose """
        if s.M is not None:
            b = s.M.dot( b )  # (M b) \ A

        x = s.solver( b )  # b \ A

        if s.verbose >= 2 or s.solver.info != 0:
            t = time.time() - s.t0  # wallclock
            dx = x - s.xprev
            niter = s.solver.info  # lgmres gcrotmk maxiter / patch niter
            niter = "\t%d iter" % niter if niter  else ""  # patch lgmres gcrotmk
            print(
"invop %s %3.0fs  |b-Ax|/|b| %-8.2g  dx %s\t x %s  %s " % (
                s.solvername,
                t,
                s.solver.resbyb, lohi(dx), lohi(x),
                niter ))
            s.xprev = x
        # todo: if totaltime >= s.maxtime:

        return x


    matvec = _matvec = dot = __call__  # ?

    @property
    def shape( s ):
        return s.A.shape

    @property
    def dtype( s ):
        return s.A.dtype


def Aminus( A, c, M=None, inplace=False ):
    """ A - c I or A - c M  for Invop( Aminus( A, c )), shift-invert """
    if c == 0:  # better -1e-6 if there's a cluster very near 0
        return A
    if isinstance( A, LinearOperator ) \
    or isinstance( M, LinearOperator ):
        if M is None:
            M = interface.IdentityOperator( A.shape )
        return (aslinearoperator( A )
            - c * aslinearoperator( M ))  # arpack: gmres
    if M is None:
        n = A.shape[0]
        M = sparse.eye( n ) if sparse.issparse( A ) \
            else np.eye( n )
    if inplace:
        A -= c * M
    else:
        A = A - c * M
    return A

sparserotate.py

""" sparserotate( eigvals ): block_diag, blocks 2^p x 2^p, with the same eigvals
    keywords: eigenvalues, sparse rotation, sparse-matrix, Hadamard-matrix, python
"""
    # is eigsh( sparserotate( diag )) slower than eigsh( diag ) ?
    # arpack powerlaw n=20k p=10 ~ 3* slower

import numpy as np
from scipy import sparse

#...........................................................................
def sparserotate( eigvals, p=2 ):
    """ :param eigvals: 1d arraylike, len n a multiple of 2^p
        :return: n x n sparse.block_diag() with the same eigenvalues as diag( eigvals )
            2^p x 2^p blocks, 2^p nnz per row
    """
        # random subsets, ix_ ?
    n = len(eigvals)
    assert eigvals.ndim == 1 and (n % 4 == 0), \
            [type(eigvals).__name__, eigvals.shape]
    H = hadamard( p=p ) / (2 ** (p / 2))
    blocks = sparse.block_diag([ (H * eig4) @ H
                                for eig4 in eigvals.reshape( -1, 2**p )])
    return blocks.tocsc()


_H2 = np.array( [[1., 1],
                [1, -1]] )

def hadamard( p=2 ):
    """ -> H 2^p x 2^p +- 1
        eigvals( H' Diag H' ) = eigvals( Diag ), H' = H / 2^(p/2)
    """
    return _H2 if p <= 1 \
        else np.kron( _H2, hadamard( p - 1 ))


#...........................................................................
if __name__ == "__main__":
    import sys
    from scipy.linalg import eigvalsh, norm
    from scipy.sparse.linalg import eigsh

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=120,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g

    n = 10000  # 10k 20k 30k: 2 7 13 sec
    p = 2  # blocks 2^p x 2^p
    powerlaw = 10
    do = "ar dense"
    k = 100

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

    n = (n // 2**p) * 2**p
    print( "%s  n %d  p %d  powerlaw %g  do %s  k %d " % (
            sys.argv[0], n, p, powerlaw, do, k ))

    D0 = 10 ** (-10 / powerlaw)
    D = np.linspace( D0, 1, n ) ** powerlaw  # 1e-10 .. 1
    print( "D:", D )

    A = sparserotate( D, p=p )
    if p <= 3:
        print( "A: \n", A[:16,:16].A )  # blocks

    do = do.split()
    if "ar" in do:
        arevals, arV = eigsh( A, sigma=0, k=k, tol=1e-5, v0=np.ones(n) )
        print( "arpack shift-invert %d evals: \n%s " % (k, arevals ))

    if "dense" in do:
        evals = eigvalsh( A.A )  # 10k: 82 sec real * 384% of 4 cores
        print( "all evals dense:", evals )
        print( "max |eigvalsh( sparserotate D ) - D|: %.2g" % 
                norm( evals - D, ord=np.inf ))  # ~ eps

sparsesolvers.py

#!/usr/bin/env python3
""" solver = Sparsesolver( solvername, A, tol= ... )
    x = solver( b ): Ax ~ b
    direct or iterative solvers spsolve splu gmres lgmres ...
        posdef (non-neg def): cholesky lsq, *cg*
    verbose=1 tracks iterations
    A = Invop( A, solver= ) wraps this for Arpack,
    see https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html
"""
# not much error checking
# W gmres: fast convergence when the eigenvalues of A are clustered away from the origin
#   else verrry slow -- see gist gmres vs. Lgmres
# https://scicomp.stackexchange.com/questions/tagged/linear-solver 300
# https://scicomp.stackexchange.com/questions/tagged/gmres 40

import time
import numpy as np
from numpy.linalg import norm
from scipy import sparse
from scipy.sparse import linalg as sparselin
try:
    from scikits import umfpack
    umfpack_version = umfpack.__version__
except ImportError:
    umfpack = umfpack_version = None
try:
    from sksparse.cholmod import cholesky, cholesky_AAt, CholmodNotPositiveDefiniteError
    from sksparse import __version__ as cholesky_version
    # https://scikit-sparse.readthedocs.io/en/latest/cholmod.html
except ImportError:
    cholesky = cholesky_version = None

__version__ = "2020-10-03 Oct"  # 06-02 June  # denis

#...............................................................................
_directsolvers = "splu spsolve "
_cgsolvers = "bicg bicgstab cg cgs minres "  # posdef / symmetric
_itersolvers = "lgmres gcrotmk qmr gmres " + _cgsolvers
_posdefsolvers = "cholesky lsq "
_allsolvers = "%s %s %s" % (_directsolvers, _itersolvers, _posdefsolvers)

#...............................................................................
class Sparsesolver( object ):
    """ Ainv = Sparsesolver( "splu", A ... ) for arpack eigs( OPinv or Minv ) """
    def __init__( s,
            solvername,     # "spsolve" "splu" for matrices, "lgmres" ... LinOps
            A,              # usually sparse csc, or LinOp
            x0=None,        # None -> 0 / scalar * ones
            tol=1e-5,
            atol=None,      # None: tol
            maxiterfac=10,  # maxiter = this * n, default 10 n
            maxouter=1000,  # lgmres maxiter, default 1000
            cholinc=0,      # cholesky (A + cholinc I) nudge posdef
            prepend_outer_v=True,   # lgmres gcrotmk ?
            useumfpack=True,    # spsolve: umfpack if installed, else splu
            verbose=0 ):
        s.solvername = solvername
        s.A = A
        n = A.shape[0]
        if np.isscalar( x0 ):  # most iter solvers None -> 0
            x0 = x0 * np.ones( A.shape[0], dtype=A.dtype )
        s.x0 = x0
        s.tol = tol  # iterative.py |b - Ax| < max( tol |b|, atol )
        if atol is None:
            atol = tol
                # / sqrt n ?  norm(res) < atol  ==  rms(res) < atol * sqrt n
        s.atol = atol
        if solvername in ["gcrotmk", "lgmres"]:
            s.maxiter = maxouter  # default 1000
            s.outeriter = "outer "
        else:
            s.maxiter = maxiterfac * n  # default 10 * n
            s.outeriter = ""
        s.verbose = verbose

        if not sparse.issparse( A ):
            assert solvername in _itersolvers.split(), "%s can't do LinOps" % solvername

        #.......................................................................
        if solvername in _itersolvers.split():
            solver = getattr( sparselin, solvername )
            kw = dict( prepend_outer_v=prepend_outer_v ) if solvername == "lgmres" \
                else {}
            solve = lambda b: \
                    solver( A, b, x0=x0, tol=tol, atol=atol,
                        maxiter=s.maxiter, **kw )

        elif solvername == "spsolve":
            solve = lambda b: \
                    sparselin.spsolve( A, b )
            # except UmfpackWarning: (almost) singular matrix
            sparselin.use_solver( useUmfpack=useumfpack )
            if umfpack is None:
                print( "Warning: umfpack is not installed, splu will be used instead" )
                # splu: slow build, then fast solves

        elif solvername == "splu":
            solve = sparselin.splu(A).solve  # options ...
        elif solvername == "spilu":
            solve = sparselin.spilu(A).solve

        elif solvername == "cholesky":  # sparse posdef, very fast
            A = A.tocsc()
            solve = cholesky( A, beta=cholinc )  # A + beta I
                # except CholmodNotPositiveDefiniteError
            """ If A has tiny eigenvalues, 0 to within machine precision,
                about half of these "zeros" may be negative -- asking for trouble.
                So nudge them > 0, say A + 1e-6 or so;
                this also improves the condition number to ~ rho(A) / beta.
            """
        elif solvername == "lsq":  # cholesky (A'A + beta I) x = A'b
            At = A.T.tocsc()
            solve = lambda b: \
                    cholesky_AAt( At, beta=cholinc )( At.dot( b ))

        elif solvername == "dense":  # for comparison -- 10k^2 * 8 = .8 gb
            Adense = getattr( A, "A", A )  # .astype( np.float32 / complex64 )
            solve = lambda b: \
                    scipy.linalg.solve( Adense, b )
                        # scipy assume_a= "gen" / "sym" / "her" / "pos"

        else:
            assert 0, "no solver " + solvername

        s.solve = solve
        s.ncalls = s.sumtime = s.nstop = 0
        s.res = s.resbyb = np.NaN


    #...........................................................................
    def __call__( s, b, **kw ):
        """ solve(b), verbose: print solver time b-Ax xmin xmax [per iter] """
        t0 = time.time()  # wall clock, some do multicore

        x = s.solve( b )

        if isinstance( x, tuple ):  # iter solvers
            x, s.info = x
            s.nstop += (s.info == s.maxiter)  # lgmres gcrotmk
        else:
            s.info = 0

        t = time.time() - t0
        s.ncalls += 1
        s.sumtime += t - t0
        s.res = b - s.A.dot(x)
        s.resbyb = norm(s.res) / max( norm(b), 1e-20 )
        if s.verbose:
            niter = ("\t%d iter" % s.info) if s.info  else ""
            print( "%6.1f sec  %-8s  |b-Ax|/|b| %-8.2g  tol %-4.2g  x %g .. %g %s " % (
                t, s.solvername, s.resbyb, s.tol,
                x.real.min(), x.real.max(),
                niter ))

        return x

    
    def stopstr( s ):
        return "" if s.nstop == 0 \
            else "%d stopped at maxiter %s " % (
                    s.nstop, s.maxiter )

test-mikotapair.py

#!/usr/bin/env python3
""" test-mikotapair.py  K x = λ M x 
    KM_smalleigs (shift-invert), KM_bigeigs
"""

import sys
import time
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import eigsh

from invop import Invop, Aminus
from eigsutil_ import eigcheck, eigsort
from mikotapair import mikotapair  # K, M eigvals 1 4 9 ... easy

#...............................................................................
def KM_smalleigs( K, M, sigma=0, k=10, tol=1e-6, **kw ):
    """ λ nearest sigma, shift-invert:
        find biggest λ with λ x = Kinv M x
        shift back
    """
    print( "KM_smalleigs, shift-invert: eigsh( OPinv=Invop(K), M=M, sigma=%g )" % sigma )
    assert sigma is not None, sigma
    Kshift = Aminus( K, sigma, M, inplace=True )
    Kinv = Invop( Kshift, tol=tol, **kw )  # kw: solver=solver ...
    t0 = time.time()
    A = sparse.csc_matrix( K.shape )  # arpack with OPinv uses A.dtype only, not A

    evals, V = eigsh( A, OPinv=Kinv, M=M, sigma=sigma, k=k, tol=tol, v0=v0 )
    evals, V = eigsort( evals, V )
    print( "evals:", evals )
    diff, maxdiffs = eigcheck( K, evals, V, M=M )
    t = time.time() - t0
    print( "\neigsh: %.1f sec  %d calls of %s%s" % (
            t, Kinv.solver.ncalls, solver, Kinv.solver.stopstr() ))
    return evals, V

#...............................................................................
def KM_bigeigs( K, M, k=10, tol=1e-6, **kw ):
    """ large λ, K x = λ M x  ->  Minv K x = λ x """
    print( "KM_bigeigs: eigsh( K, M=M, Minv=Invop(M)" )
    Minv = Invop( M, tol=tol, **kw )  # kw: solver=solver ...
    t0 = time.time()

    evals, V = eigsh( K, M=M, Minv=Minv, sigma=None, k=k, tol=tol, v0=v0 )
    evals, V = eigsort( evals, V )
    print( "evals:", evals )
    diff, maxdiffs = eigcheck( K, evals, V, M=M )
    t = time.time() - t0
    print( "\neigsh: %.1f sec  %d calls of %s%s" % (
            t, Minv.solver.ncalls, solver, Minv.solver.stopstr() ))
    return evals, V


#...............................................................................
if __name__ == "__main__":

    np.set_printoptions( threshold=10, edgeitems=5, linewidth=140,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    print( 80 * "▄" )
    # print( nu.versionstr() )

    #...........................................................................
    test = "mikotapair"  # evals 1 4 9 .. n^2, K tridiag M diagonal, easy
    n = 1e6
    incdiag = 0

    solver = "cholesky"
    # solvers = "spsolve splu lgmres gmres "  # arpack shift-invert: gmres / splu
    cholinc = 0  # nudge posdef
    maxiterfac = 10  # 10 n
    maxouter = 5000  # lgmres gcrotmk 1000
    verbose = 2  # 2: each call Invop to solver

        # eigs args --
    sigma = 0  # 0: KM_smalleigs, None: bigeigs
    k = 6
    tol = 1e-6
    v0 = 1  # "laplace"

        # 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 )
    n = int( n )
    k = min( k, n - 1 )

    params = "%s  test %s  n %d  incdiag %.3g  sigma %s  k %d  tol %g  v0 %s  maxiterfac %d  maxouter %d  cholinc %g" % (
            sys.argv[0], test, n, incdiag, sigma, k, tol, v0, maxiterfac, maxouter, cholinc )
    print( "params:", params )

    #...........................................................................
    K, M = mikotapair( n, flip=False )  # evals 1 4 9 .. n^2, K tridiag M diagonal, easy
        # A = testmatrix( test, n=n, incdiag=incdiag, p=p, lo=lo, hi=hi )
    K = K.tocsc()
    n = K.shape[0]
    k = min( k, n - 1 )
    if v0 == 1:
        v0 = np.ones( n )
    else:
        v0 = np.random.default_rng(0) .laplace( size=n )
        # *much* better: from previous runs or smaller problems

    # denseevals, denseV = scipy.linalg.eigh( K.A, M.A )  # K x = λ M x
    # print( "\ndenseevals:", denseevals )
    invopargs = dict( solver=solver, verbose=verbose,
            maxiterfac=maxiterfac, maxouter=maxouter, cholinc=cholinc )

    if sigma is None:
        evals, V = KM_bigeigs( K, M, k=k, tol=tol, **invopargs )
    else:
        evals, V = KM_smalleigs( K, M, sigma=sigma, k=k, tol=tol, **invopargs )

testmatrix.py

#!/usr/bin/env python3
"""testmatrix: load or generate a sparse matrix, with optional transforms """

import numpy as np
from scipy import sparse
from scipy.io import mmread

    # matrix generators --
try:
    from poisson_nd import poisson1, poisson2, poisson3
    from powerlaw import powerlaw_eigs
    from randomsparse import randomsparse
    from Lband import LLband, Ltriband
    from mikotapair import mikotasingle
    from sparserotate import sparserotate
except ImportError:
    pass

from eigsutil import eigcheck, flip_eigmax
from zutil import arraydesc, findfile, globs, maxnorm

#...............................................................................
def testmatrix( name, n=0, incdiag=0, testsymm=True, **kw ):
    """ load or generate a sparse csc matrix 
        :param "name": poisson* powerlaw randomsparse ...
            or filename*.npz filename*.mtx -- * $ ~ are expanded
            or A/B -> $suitesparse/A/B/*.mtx
            or "name-flip-Ltri ...": optional transforms
        :return: a sparse csc matrix
    """
    s = name.replace( "-", " " ).split()
    nm, trans = s[0], s[1:]  # 0 or more transforms

    A = _getmatrix( nm, n=n, **kw )
    n = A.shape[0]
    if trans:
        A = transforms( trans, A )
    if incdiag:
        A += incdiag * sparse.eye( n )
    A = sparse.csc_matrix( A )
    if testsymm:
        A_AT = maxnorm( A - A.conj().T )
        symm = "max |A - A'|: %.2g" % A_AT
    else:
        symm = ""
    print( "\ntestmatrix %s: A %s  incdiag %.3g  %s " % (
            name, arraydesc( A ), incdiag, symm ))
    print( "  A.diag:", A.diagonal() )
    return A

#...............................................................................
def transforms( trans, A ):
    """ :param trans: "flip Ltri LL rotate float32"
        :param A: sparse matrix
        :return: e.g. Ltri( flip( A ))
    """
    if not trans:
        return A
    if isinstance( trans, str ):
        trans = trans.split()
    if "flip" in trans:
            A, eigmax = flip_eigmax( A )  # eigmax I - A, emin <-> emax
    n = A.shape[0]
    if "Ltri" in trans:
            A = Ltriband( n ).tocsc() .dot( A )  # lower-triangular, diags 0 1 2 4 ... ones
    if "LL" in trans:
            L = Ltriband( n ).tocsc()
            A = L .dot( A ) .dot( L.T )
    if "rotate" in trans or "sparserotate" in trans:
            A = sparserotate( A.diagonal() )
    if "float32" in trans:
        A = A.astype( np.float32 )
    return A

#...............................................................................
def _getmatrix( name, n=0, **kw ):
    """ generate a matrix, poisson* / powerlaw / randomsparse ...
        or load findfile( .npz .mtx )
    """
    n = int( n )

    if name == "I":
        return sparse.eye( n ).tocsc()

    if name == "LLband":
        return LLband( n )  # L L' sparse band matrix, may have tiny eigenvalues < 0

    if name == "mikota":
        return mikotasingle( n=n )  # tridiag symmetric with eigvals n^2 - [1 4 9 ...]

    if name == "poisson1":
        return poisson1( n or 100 )  # eigvals 0 .. 4
    if name == "poisson2":
        return poisson2( n or 100 )  # n^2, n^2  eigvals 0 .. 8
    if name == "poisson3":
        return poisson3( n or 20 )  # n^3, n^3

    if name == "powerlaw":
        A, _ = powerlaw_eigs( n or 100, **kw )  # diagonal
            # def powerlaw( n, p=5, lo=1e-3, hi=1, flip=False )
        return A
    for k in "p flip lo hi ".split():
        kw.pop( k, None )

    if name == "randomsparse":
        # def randomsparse( m, n, density=4, symmetric=False, posdef=False,
        kw.pop( "incdiag", None )
        return randomsparse( n, n, **kw )
    for k in "density symmetric posdef ".split():
        kw.pop( k, None )

    if name == "curve":
        return sparse.diags( curve( n ))

    #...........................................................................
    # load files, sparse.load_npz / mmread
    if name.endswith(( ".npz", ".mtx", ".mtx.gz" )):
        filename = findfile( name )  # expand * $ ~, IOError if not found
        return load_npz_mtx( filename )

        # e.g. Norris/fv1 -> $suitesparse/Norris/fv1/*.mtx from https://sparse.tamu.edu
    A = try_npz_mtx( "$suitesparse/" + name )
    assert A is not None, "testmatrix: can't find %s" % name
    return A


def try_npz_mtx( dir ):
    g = globs( dir + "/*.npz" )  # expand * $ ~
    if g:
        return load_npz_mtx( g[0] )  # first only
    g = globs( dir + "/*.mtx" )
    if g:
        return load_npz_mtx( g[0] )
    return None


def load_npz_mtx( filename ):
    if filename.endswith( ".npz" ):
        return sparse.load_npz( filename ).tocsc()
    if filename.endswith(( ".mtx", ".mtx.gz" )):  # matrix market -- text, huge, mem leak
        return mmread( filename ).tocsc()
    return None


def curve( n,
        x = np.r_[ 0, 1, 2, 5, 10, 25, 50, 75, 90, 95, 98, 99, 99.9, 100 ],
        y = np.r_[ 0, 0, 1e-16, 0.13, 0.19, 0.33, 0.61, 1.1, 1.9, 2.6, 3.8, 5, 270, 273 ]
        ):
    """ -> spectrum like scikit-learn 20newsgroups ndoc=4000 A A' """
        # ... 13 14 18 23 273
    print( "curve: %s \n-> %s" % (x, y) )
    xnew = np.linspace( x.min(), x.max(), n )
    return np.interp( xnew, x, y )  # mono spline ?


#...............................................................................
if __name__ == "__main__":
    import sys
    import scipy.linalg

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=120, formatter = dict(
                float = lambda x: "%.2g" % x,  # float arrays %.2g
                complexfloat = " {0.real:.2g} {0.imag:+.2g}j " .format ))

    test = "powerlaw rotate"
    # test = "powerlaw flip"
    # test = "randomsparse"  # 100 A'A: -5e-16 .. -2e-18  8e-33 1e-16 ..
    # test = "I LL"
    n = 100
    incdiag = 0
    neig = 10

    # 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 = testmatrix( test, n=n, incdiag=incdiag, posdef=True )
    n = A.shape[0]
    print( A[:10,:10].A, "..." )
    print( A[-10:,-10:].A )

    if neig > 0:
        evals, V = scipy.linalg.eigh( A.A, subset_by_index=[n - neig, n - 1] )
        print( "eigh dense:", np.sort( evals ))
        eigcheck( A, evals, V )