ARPACK utilities (structsolve.arpackutils)#

Utilities for ARPACK, against the Intel MKL bug of dsteqr

With Intel MKL 2024.2.0 to 2025.0.0, scipy.sparse.linalg.eigsh() returns wrong eigenpairs, without any error, or raises ArpackError -8, whenever its number of Lanczos vectors is ncv > 32, e.g. with SciPy’s default ncv = 2 k + 1 for k >= 16 eigenpairs. MKL 2025.0.1 or newer is not affected, neither is scipy.sparse.linalg.eigs().

capped_eigsh() keeps ncv <= 32 when arpack_max_ncv() says so, computing more than 15 eigenpairs in several calls with deflation, whose eigenvectors are combined with a Rayleigh-Ritz projection such that they are orthonormal also across the calls.

The cap is decided by arpack_max_ncv() and can be overridden with the module attribute ARPACK_MAX_NCV or the environment variable STRUCTSOLVE_ARPACK_MAX_NCV, e.g. to lift it with a fixed MKL that cannot be identified as such:

import structsolve.arpackutils
structsolve.arpackutils.ARPACK_MAX_NCV = 0  # no cap

Notes#

The symmetric drivers of ARPACK, used by scipy.sparse.linalg.eigsh(), compute the Ritz vectors with the LAPACK routine dsteqr applied to the ncv x ncv tridiagonal matrix of the Lanczos process. The dsteqr of Intel MKL 2024.2.0 to 2025.0.0 returns wrong eigenvectors for matrices larger than 32 x 32, hence the wrong eigenpairs whenever ncv > 32. The non-symmetric drivers of scipy.sparse.linalg.eigs() do not call dsteqr.

structsolve.arpackutils.ARPACK_MAX_NCV = None#

Maximum ncv of capped_eigsh(). None, the default, reads the environment variable STRUCTSOLVE_ARPACK_MAX_NCV and, when it is not defined, uses 'auto': the cap DSTEQR_MAX_SIZE is applied when the LAPACK of SciPy is detected as faulty, see arpack_max_ncv(). 0 (or 'none') disables the cap, an integer of at least 3 sets it.

structsolve.arpackutils.DSTEQR_MAX_SIZE = 32#

Largest tridiagonal matrix for which the faulty dsteqr is correct

structsolve.arpackutils.MKL_FAULTY_VERSIONS = ((2024, 2, 0), (2025, 0, 0))#

First and last Intel MKL versions with the faulty dsteqr

structsolve.arpackutils.PARDISO_ACCURACY_RTOL = 1e-06#

linear_solver='pardiso': largest relative residual of the solution of a factorization of PARDISO, without refinement, for a random right-hand side, above which SuperLU is used instead, see pardiso_solve(), and the count of negative pivots of PARDISO is not trusted, see linear_buckling._negative_pivots(); that of static._factorize() for the factorizations of the non-linear analyses

structsolve.arpackutils.arpack_max_ncv()#

Maximum ncv of capped_eigsh(), None when not capped

The value of ARPACK_MAX_NCV is used, or, when it is None, the environment variable STRUCTSOLVE_ARPACK_MAX_NCV. Both accept 'auto', 0 or 'none' (no cap), or an integer of at least 3. In the 'auto' mode, the default, ncv is capped at DSTEQR_MAX_SIZE when:

The detection runs once per process.

Returns:
max_ncvint or None

The cap of ncv, at least 3, or None when ncv is not capped.

Raises:
ValueError

When the value given is an integer lower than 3 and not 0.

structsolve.arpackutils.arpack_ncv(k, n)#

ncv of scipy.sparse.linalg.eigsh() for k eigenpairs of a problem of size n

SciPy’s default min(n, max(2 k + 1, 20)), capped at arpack_max_ncv(). With a cap, k should not exceed (max_ncv - 1)//2, see capped_eigsh().

Parameters:
kint

Number of eigenpairs.

nint

Size of the problem, the number of rows of the matrices.

Returns:
ncvint

Number of Lanczos vectors.

structsolve.arpackutils.capped_eigsh(A, k, M, sigma, which, mode='normal', tol=0, linear_solver='splu')#

scipy.sparse.linalg.eigsh() with ncv from arpack_ncv()

Solves \(\boldsymbol{A} \boldsymbol{x} = w \boldsymbol{M} \boldsymbol{x}\) with a positive definite M about the shift sigma, in the shift-invert (mode='normal') or Cayley (mode='cayley') mode, selecting the eigenvalues by increasing \(w\), e.g. the Cayley mode with which='SM' and sigma larger than all \(|w|\), or the shift-invert mode with which='LM' and sigma lower than all \(w\).

When ncv is capped at max_ncv, see arpack_max_ncv(), and k is larger than (max_ncv - 1)//2, i.e. 15, the eigenpairs are computed in passes of at most 15 eigenpairs, each with ncv=max_ncv. The first pass is a call in the requested mode. The following passes use the shift-invert mode about the shift sigma_d given by _deflation_shift(), with the eigenvectors of the previous passes deflated from the operator, and each of them returns the next eigenvalues by increasing \(w\). Each pass factorizes A - sigma_d*M with SuperLU, and every call starts from start_kwargs(). The eigenvectors of all the passes are finally combined with a Rayleigh-Ritz projection, see _rayleigh_ritz(), such that they are M-orthonormal also across the passes.

With linear_solver='pardiso', A - sigma*M of the first call and A - sigma_d*M of the passes are factorized by PARDISO instead, see pardiso_solve(), and passed to scipy.sparse.linalg.eigsh() as OPinv, \((\boldsymbol{A} - \sigma \boldsymbol{M})^{-1}\) in both the shift-invert and the Cayley modes; SuperLU is used for a matrix whose factorization by PARDISO is inaccurate.

Parameters:
Asparse matrix

Square matrix of the left-hand side, of size n x n.

kint

Number of eigenpairs.

Msparse matrix

Positive definite matrix of the right-hand side, of size n x n.

sigmafloat

Shift of the shift-invert or Cayley mode, in the units of \(w\).

whichstr

Eigenvalues to find, passed to scipy.sparse.linalg.eigsh(), e.g. 'LM' or 'SM'.

modestr, optional

'normal' (shift-invert, default) or 'cayley', passed to scipy.sparse.linalg.eigsh().

tolfloat, optional

Relative accuracy of the eigenvalues, passed to scipy.sparse.linalg.eigsh(); 0 (default) means the machine precision.

linear_solverstr, optional

'splu' (default), the factorizations of SuperLU, or 'pardiso', those of PARDISO, see above.

Returns:
w, xndarray

The eigenvalues and the \(\boldsymbol{M}\)-orthonormal eigenvectors, by increasing \(w\) when computed in passes. w is a 1-D array of size k and x an array of shape (n, k).

See also

arpack_max_ncv, arpack_ncv, start_kwargs, start_vector

Notes

In the shift-invert mode ARPACK computes the eigenvalues \(\nu = 1/(w - \sigma)\) of largest modulus of the operator \((\boldsymbol{A} - \sigma \boldsymbol{M})^{-1} \boldsymbol{M}\), i.e. the \(w\) nearest \(\sigma\) ([1], Chapter 3, Section 3.2, Eq. (3.2.2)). The passes keep k <= (max_ncv - 1)//2, i.e. ncv >= 2 k + 1. With k close to ncv ARPACK converges very slowly, not at all, or misses eigenvalues, e.g. k=31 and ncv=32.

The passes after the first one use the shift-invert mode about \(\sigma_d\), the middle of the largest gap between the last five eigenvalues found, with the eigenvectors \(\boldsymbol{V}\) of the previous passes deflated from the operator:

\[\boldsymbol{y} = (\boldsymbol{I} - \boldsymbol{V} \boldsymbol{V}^\top \boldsymbol{M})(\boldsymbol{A} - \sigma_d \boldsymbol{M})^{-1} \boldsymbol{x}\]

such that the eigenvalues of \(\boldsymbol{V}\), already found, become null and are not found again. Each pass therefore returns the eigenvalues nearest \(\sigma_d\) among those not found yet, i.e. the next ones by increasing \(w\). A shift close to them, instead of sigma, separates them better: in the study of the panels cylinders it was 30 to 45% faster and gave eigenvalues 5 to 10 times more accurate. ARPACK applies the operator to \(\boldsymbol{M} \boldsymbol{x}\); for an eigenvector \(\boldsymbol{v}_i\) of \(\boldsymbol{V}\), \((\boldsymbol{A} - \sigma_d \boldsymbol{M})^{-1} \boldsymbol{M} \boldsymbol{v}_i = \boldsymbol{v}_i/(w_i - \sigma_d)\) and, \(\boldsymbol{V}\) being \(\boldsymbol{M}\)-orthonormal, \((\boldsymbol{I} - \boldsymbol{V} \boldsymbol{V}^\top \boldsymbol{M}) \boldsymbol{v}_i = \boldsymbol{0}\).

References

[1]

Lehoucq, R.B., Sorensen, D.C., Yang, C., ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998. https://doi.org/10.1137/1.9780898719628

structsolve.arpackutils.dsteqr_is_faulty()#

Tell whether the LAPACK dsteqr of SciPy returns wrong eigenvectors

scipy.linalg.lapack.dstev computes the eigenvectors of a symmetric tridiagonal matrix with dsteqr. It is called for three random matrices of sizes 33, 40 and 64, whose eigenpairs must have relative residuals and an orthogonality error below 1e-10. With Intel MKL 2025.0.0 every size above 32 fails, with rare exceptions for sizes close to 32, hence the three sizes.

Returns:
faultybool or None

None when dstev is not available.

structsolve.arpackutils.lapack_is_mkl()#

Tell whether the LAPACK of SciPy is Intel MKL

scipy.show_config(mode='dicts') is read first, NumPy’s build configuration only when SciPy’s does not name the LAPACK library. Older versions of SciPy without mode='dicts' are checked through scipy.__config__.

Returns:
is_mklbool or None

None when the build configuration does not tell.

structsolve.arpackutils.mkl_version()#

Version of the Intel MKL loaded in the process, as a tuple (major, minor, update)

The build configuration of SciPy only gives the major version, e.g. '2025', therefore the version is read at runtime from the optional packages mkl-service (mkl.get_version()) or threadpoolctl.

Returns:
versiontuple or None

None when neither package is installed or no MKL is loaded.

structsolve.arpackutils.pardiso_factorization(A)#

Factorization of the symmetric sparse matrix A by PARDISO and the relative residual of its solution for a random right-hand side

A is factorized by static._factorize() with linear_solver='pardiso': the symmetric indefinite factorization of PARDISO of the matrix scaled to a unit diagonal, with the fallback on perturbed pivots, see its Notes.

Returns:
factorizationcallable or None

The solve_kT of static._factorize(), with num_negative_eigvals and num_perturbed_pivots; None when A has a null row, i.e. it is exactly singular, PARDISO failed, or A is unsymmetric for static._factorize().

residualfloat or None

\(\|\boldsymbol{A} \boldsymbol{x} - \boldsymbol{b}\|/\|\boldsymbol{b}\|\) of the solution \(\boldsymbol{x}\) of the factorization, without refinement, for a fixed random \(\boldsymbol{b}\); None with factorization.

structsolve.arpackutils.pardiso_solve(A, refine=True)#

Solution function of the symmetric sparse matrix A factorized by PARDISO, for the eigenvalue analyses with linear_solver='pardiso'

A is factorized by pardiso_factorization(). With refine=True every solution whose relative residual is larger than static.PARDISO_REFINE_TOL is improved by GMRES, see static._factorize(); with refine=False the solution of the factorization is returned as is, as by SuperLU. When PARDISO fails, or its factorization is inaccurate, i.e. the relative residual of pardiso_factorization() is larger than PARDISO_ACCURACY_RTOL, A is factorized by SuperLU instead, with a RuntimeWarning.

Returns:
solvecallable

solve(b) returns \(\boldsymbol{A}^{-1} \boldsymbol{b}\) for a real or complex b of shape (n,); its attribute linear_solver is 'pardiso' or 'splu', the solver used, and its attribute factorization the solve_kT of static._factorize() with PARDISO, None with SuperLU, which holds the memory of PARDISO, freed when it is deleted.

Raises:
RuntimeError

When SuperLU is used and A is exactly singular.

ImportError

When oneMKL is not found.

structsolve.arpackutils.release_memory()#

Collect the reference cycles of SciPy’s ARPACK wrappers in the browser

The objects of scipy.sparse.linalg.eigsh() and scipy.sparse.linalg.eigs() keep the SuperLU factorizations of the shift-invert modes in reference cycles, freed only by the cyclic garbage collector. Natively this is harmless, but the WebAssembly heap of Pyodide (sys.platform == 'emscripten') is exhausted before the collector runs, e.g. after about ten eigs calls on a problem with 1,000 dofs, raising MemoryError. gc.collect() is then called after each ARPACK call; it does nothing on other platforms.

structsolve.arpackutils.start_kwargs(n)#

Arguments of scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh() such that the results are reproducible

v0 is start_vector() and, with SciPy 1.17 or newer, rng is a numpy.random.Generator with the seed 42, from which ARPACK draws a new starting vector whenever the Krylov subspace is invariant, e.g. when the operator is nilpotent.

Parameters:
nint

Size of the problem.

Returns:
kwargsdict

{'v0': start_vector(n)}, with the key 'rng' when supported.

Notes

Since SciPy 1.17, ARPACK asks the caller for a random vector when the Arnoldi or Lanczos process breaks down, i.e. finds an invariant subspace, and SciPy draws it from rng. With rng=None it comes from the entropy of the operating system: with the same v0 the results then differ from one call to the next. Older versions of SciPy draw it inside ARPACK, with a seed kept from one call to the next in the same process, which rng cannot fix.

structsolve.arpackutils.start_vector(n)#

Fixed starting vector of ARPACK, uniform in [-1, 1] as ARPACK’s random one, such that the results are reproducible

The vector is drawn from numpy.random.RandomState with the seed 42, so that the same n always gives the same vector. The calls of ARPACK use start_kwargs(), which also fixes the vectors drawn when ARPACK restarts.

Parameters:
nint

Size of the problem.

Returns:
v0ndarray

1-D array of size n.