Arnoldi

This is an attempt to write an eigensolver for sparse matrices that only relies on NumPy and BLAS/LAPACK, without depending on ARPACK. Ultimately, the hope is to be a viable replacement for scipy.sparse.linalg.eigs, and remove the fortran dependency.

How to install

The project is using uv. Assuming you have uv installed, simply do

uv sync

Example of usage

import numpy as np

from scipy.sparse.linalg import eigs

from arnoldi.krylov_schur import partial_schur
from arnoldi.matrices import mark
from arnoldi.utils import arg_largest_real


# Markov walk matrix
A = mark(50)

TOL = 1e-8
MAX_RESTARTS = 1_000
MAX_DIMS = 20
K = 5
WHICH = "LR"
SORT_FUNCTION = arg_largest_real

Q, T, history = partial_schur(
    A,
    K,
    max_dim=MAX_DIMS,
    stopping_criterion=TOL,
    sort_function=SORT_FUNCTION,
    max_restarts=MAX_RESTARTS,
)
# Convert Schur transform into eigenpairs
vals, S = np.linalg.eig(T)
vecs = Q @ S

print("========== true residuals, our implementation =============")
for k in range(K):
    val, vec = vals[k], vecs[:, k]
    residual = np.linalg.norm(A @ vec - val * vec)
    norm_residual = residual / np.abs(val)
    print(
        f"Residual {k}, value is {np.real(val):.7}, res is {residual:.3e}, norm res is {norm_residual:.3e}"
    )

# Compare to ARPACK
vals, vecs = eigs(A, K, ncv=MAX_DIMS, maxiter=MAX_RESTARTS, which=WHICH,
                  tol=TOL)

print("========== true residuals, ARPACK (scipy) implementation =============")
for k in range(K):
    val, vec = vals[k], vecs[:, k]
    residual = np.linalg.norm(A @ vec - val * vec)
    norm_residual = residual / np.abs(val)
    print(
        f"Residual {k}, value is {np.real(val):.7}, res is {residual:.3e}, norm res is {norm_residual:.3e}"
    )

It is not recommended for real use yet. I expect the API to change, especially in terms of convergence tracking and history.

For more complex examples, look at the scripts SLEPC script and ARPACK script. It can compare ARPACK and SLEPc (to be installed separately) against this implementation for matrices in the sparse suite.

Why ?

ARPACK-NG is a fortran library for sparse eigen solvers. It has the following issues:

• the fortran code is not thread-safe. In particular, it is not re-entrant
• it does not incorporate some of the more recent improvements discovered for large problems, e.g. A Krylov-Schur Algorithm for Large Eigenproblems, G. W. Stewart, SIAM J. Matrix Anal. Appl., Vol. 23, No. 3, pp. 601–614

AI policy

AI tools such as Claude are extensively used for code review, literature review and understanding alternative, license-compatible implementations.

However, as the code is a candidate to be incorporated into SciPy codebase, the code itself is not generated using AI, and all the code has been manually written from the literature. Exceptions are trivial scripts in the scripts/ directory. Files AI-generated will contain a header mentioning this.

TODO

For a first 1.0 release:

• Fundamental support for arbitrary matrices, largest eigen values only
  • basic arnoldi decomp w/ tests
  • add a key set of test matrices, using sparse matrix suite + synthetic (Markov, Laplace, etc.)
  • convergence tracking on Ritz values
  • explicit restart support with deflation
  • Krylov-schur method
  • customizable orthonormalization
• Try to reach parity w/ ARPACK within 2x speed/memory usage
  • compare MGS vs double GS w/ DGKS vs others in terms of precision
  • implement scenario loading from csv
  • implement dynamic p
  • implement locking
• Tracking capabilities
  • callback setup to enable 1) tracking key operations and 2) track invariants
  • implement time spend tracking in our implementation / SLEPCs and profiling support for scipy eigs (ARPACK)
  • try alternative methods (e.g. FEAST) for cases SLEPC/ARPACK fail on (e.g. olm5000 LR nev >= 22)
• Simple to add
  • LinearOperator support -> enable SHIFT/Inverse mode (S* and sigma support)
  • handle happy breakdown in Krylov-Schur
  • Customizable convergence support + 2 implementations (currnet + mix relative/absolute ala ArnoldiMethod)
  • Single precision ?

Post 1.0:

• optimize for the case of Hermitian/symmetric matrices (Lanczos)
• add support for shift-invert (arbitrary eigen values)
• add support for general inverses problems
• single precision support
• optimization:
  • optimize orthonormalization
  • ensure memory space is mostly V and nothing else as function of input size
  • block Krylov-Schur ?

Existing alternative implementations

• julia
  • The Arnoldi method with Krylov-Schur restart, natively in pure Julia
    • According to https://discourse.julialang.org/t/ann-arnoldimethod-jl-v0-4/110604, it is more stable than ARPACK.
    • implementation ~ 2k julia (not including tests)
    • license compatible with SciPy (MIT-like)

References

• Mathematical references:
  • Lecture Notes on Solving Large Scale Eigenvalue Problems by Prof. Dr. Peter Arbenz from the Computer Science Department at ETH. See in particular the latest version of the notes
  • Templates for the Solution of Algebraic Eigenvalue Problems: a Practical Guide. Overview on numerical methods for eigen values, including dense and sparse
  • A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems.
  • Applied Numerical Linear Algebra
  • Implicit application of polynomial filters in a k-step Arnoldi method
    • An Implicitly Restarted Lanczos Method for Large Symmetric Eigenvalue Problems: a specialization of the implicit ARNOLDI to the hermitian operator case.
  • Thick-Restart Lanczos Method For Large Symmetric Eigen Values Problems: a special case of Krylov-Schur for Hermitian operators
  • NUMERICAL METHODS FOR LARGE EIGENVALUE PROBLEMS: 2nd edition, only covers up to early 1990s techniques (explicit/implicit restarted Arnoldi). It explains clearly deflation, locking and gives some numerical examples that can be used as reference.