Sparse_Numba

A Numba-compatible sparse linear solver: SuperLU, Eigen (SparseLU and SimplicialLDLT) and a user-supplied UMFPACK behind one interface that runs inside @njit(nogil=True) code and keeps its numba cache across processes.

• New to it? The user guide walks through the backends, the orderings, the handle's life, numba use and the cache.
• Scripts to run and read: examples/.
• What 0.2.0 adds and changes over 0.1.11: the release notes.

0.2.0 vs 0.1.11 at a glance (the largest test matrix, n = 20698; one machine; details and every size in Performance):

• Solve: 3969 us in 0.1.11 (SuperLU, NATURAL only) -> 1426 us with the new ordering choice (SuperLU COLAMD, 2.8x) and 755 us with the new Eigen backend (Eigen LU COLAMD, 5.3x); SuperLU NATURAL stays bit-identical to 0.1.11 on the same CPU and BLAS kernel.
• Factorize: 68961 us -> 16450 us (SuperLU COLAMD, 4.2x) and 12585 us (Eigen LU COLAMD, 5.5x); the new refactorize re-factors the same pattern in 10138 us (SuperLU COLAMD).
• Overhead of the numba wrapper against a direct C++ call (sparse_numba / direct, ratio of medians): 0.85 to 1.23 at the largest size and 0.77 to 2.17 over all four sizes, in both directions; in 19 of the 102 measured cells sparse_numba's and the direct call's min - max ranges do not overlap.

Why Sparse_Numba?

Sparse linear solvers in Python (SciPy, CVXOPT, KVXOPT) are called from the interpreter. Sparse_Numba's functions are njit(nogil=True, cache=True): they are called from Numba-compiled code, release the GIL, and so run in parallel threads (prange) without a round trip through Python.

Installation

pip install sparse-numba

One wheel per Python (CPython 3.10 - 3.13) and platform, each built and tested on its own Python. Requires numba >= 0.61 (the first numba with Python 3.13; it needs Python 3.10) and numpy >= 1.24. Python 3.9 users stay on 0.1.11.

Intel Mac users: pip install sparse-numba --no-binary=sparse-numba. SuperLU's version differs per platform: 7.0.0 in the Windows wheels, 5.x in the Linux wheels (AlmaLinux 8's package, bundled by auditwheel into sparse_numba.libs), Homebrew's on macOS. The Linux wheels are manylinux_2_28: they need glibc >= 2.28 (e.g. Ubuntu 18.10+, Debian 10+, RHEL / AlmaLinux / Rocky 8+); on an older glibc pip builds from the sdist (section "Building from source").

UMFPACK: bring your own SuiteSparse

No SuiteSparse binary ships with sparse_numba (UMFPACK and CHOLMOD are GPL / LGPL; see Licences below). The UMFPACK backend loads the user's library at its first call:

• SPARSE_NUMBA_UMFPACK_LIBRARY=/path/to/libumfpack.{dll,so,dylib} names it (only that path is tried); otherwise
• the platform's library names are tried through the platform loader: libumfpack.dll, umfpack.dll (Windows: the default DLL directories, among them those added with os.add_dll_directory; sparse_numba adds SUITESPARSE_BIN when it is set); libumfpack.so.6, .so.5, .so (Linux); libumfpack.6.dylib, .5.dylib, .dylib (macOS).

Not found: every UMFPACK call returns BACKEND_ABSENT (-20), backends() leaves it out, and require(UMFPACK) raises with this hint. Nothing is printed at import.

Building from source

Windows needs MinGW-w64 (x86_64-posix-seh) on PATH; setup.py selects the mingw32 compiler itself. Then:

python -m build --wheel
pip install dist/sparse_numba-<VERSION>-<TAGS>.whl

Linux and macOS need the SuperLU and OpenBLAS development packages (and SuiteSparse's headers for the UMFPACK interface). Build switches (setup.py, read at build time):

The release flags: -O3 -DNDEBUG on every extension and platform (plus -DMS_WIN64 on 64-bit Windows, -fPIC on Linux / macOS); the Eigen extension adds -std=c++17 -DEIGEN_MPL2_ONLY -DEIGEN_DONT_PARALLELIZE; MinGW links with -Wl,--no-insert-timestamp (a rebuild of unchanged source gives the same extension bytes). The flags compiled in are part of cache_key().

More: the installation guide (wheels, UMFPACK, building from source, checking the install) and the user guide (its parts 1, Install, and 8, UMFPACK's library).

The API

Every function below is njit(nogil=True, cache=True) unless marked python-level. Input: CSC of an n x n matrix, data float64, indices and indptr int32; b / out float64, 1-D, length n, C-contiguous; B / OUT (n, k) C-contiguous.

ADDED in 0.2.0: the backend-neutral interface

Backends: SUPERLU, UMFPACK, EIGEN_LU, EIGEN_LDLT. Orderings: NATURAL, COLAMD, AMD, MMD_AT_PLUS_A. There is no default backend and no default ordering: both are arguments of factorize.

import numpy as np
from numba import njit
import sparse_numba as sn

# a 3 x 3 CSC matrix
data = np.array([4.0, 1.0, 1.0, 4.0, 1.0, 1.0, 4.0])
indices = np.array([0, 1, 0, 1, 2, 1, 2], dtype=np.int32)
indptr = np.array([0, 2, 5, 7], dtype=np.int32)
n = 3

handle, info = sn.factorize(data, indices, indptr, n, sn.SUPERLU, sn.COLAMD)
assert info == 0, sn.explain(info)

@njit(nogil=True)
def many_solves(handle, rhs, out):
    status = 0
    for i in range(rhs.shape[0]):
        status |= sn.solve_into(handle, rhs[i], out[i])
    return status

rhs = np.random.rand(1000, n)
out = np.empty_like(rhs)
assert many_solves(handle, rhs, out) == 0

data2 = data * 1.5                         # new values, the same pattern
assert sn.refactorize(handle, data2) == 0
assert sn.free(handle) == 0

The same with Eigen's SparseLU at AMD, or SimplicialLDLT on a symmetric positive definite matrix: sn.factorize(..., sn.EIGEN_LU, sn.AMD), sn.factorize(..., sn.EIGEN_LDLT, sn.AMD).

What each backend supports

"refused" = the call returns ORDERING_UNSUPPORTED or BAD_SHAPE by name, never a silent fallback.

Notes:

• EIGEN_LDLT checks that the PATTERN is symmetric (the stored entries, explicit zeros among them), not the values; it reads the lower triangle. An indefinite matrix is refused after the factorization: SINGULAR (= 1) when an entry of D is <= 0.
• A row index repeated inside one column: the Eigen backends sum the entries; UMFPACK refuses it (BACKEND_ERROR); SuperLU passes it to dgstrf as given.
• SINGULAR: SuperLU and UMFPACK return k > 0, the zero pivot's 1-based index; the Eigen backends return 1.

The status codes

info == 0 is success; info != 0 is failure (the 0.1.11 contract). sparse_numba.errors names every code and explain(info) says it in words; nothing is printed by any call.

The handle is an int64 registry index: the caller who factorized owns it and frees it; it is process-local (it never crosses processes).

KEPT from 0.1.11 (same names, signatures and return shapes)

superlu_solve_csc / _coo / _csr, superlu_factorize_csc / _coo / _csr, superlu_solve_factored, superlu_free_factors, and the same umfpack_* set; convert_coo_to_csc, convert_csr_to_csc, convert_coo_to_csr, sparse_matvec_csr, validate_sparse_matrix, check_matrix_properties, validate_csc_matrix; is_slu_available, is_umf_available. They are thin shims over the interface above; SuperLU's give bit-identical solutions to 0.1.11 (NATURAL ordering) on the same CPU and BLAS kernel (OpenBLAS picks its kernel by CPU at run time; CI pins OPENBLAS_CORETYPE=Haswell for this check). The values some codes carry changed (the 0.2.0 release notes, docs/RELEASE_NOTES_0.2.0.md).

from sparse_numba import superlu_factorize_csc, superlu_solve_factored, superlu_free_factors

b = np.ones(n)
handle, info = superlu_factorize_csc(data, indices, indptr)
x, info = superlu_solve_factored(handle, b)
superlu_free_factors(handle)

DEPRECATED (kept through 0.2.x, removed at 0.3)

*_solve_factored(handle, b) -> (x, info) (allocates x) in favour of solve_into; the *_factorize_* forms without an ordering in favour of factorize(..., backend, ordering); *_free_factors in favour of free. The deprecation is in the documentation only (njit functions do not warn).

The numba cache across processes

• Every function is cache=True; the C entry points are registered with llvmlite under stable names (sparse_numba.symbols.TABLE) at import, so a cached caller loads in a new process without recompiling.
• No cache file ships in a wheel (no __pycache__, *.nbi, *.nbc).
• A read-only install (site-packages not writable): numba writes its cache where NUMBA_CACHE_DIR points.
• Upgrades: the symbol names are stable across releases; the contract for a caller that caches its own compiled code is to fold sparse_numba.cache_key() into its cache key (it names the version, the symbol table, each extension's sha256, the SuperLU / UMFPACK / Eigen versions and the compile flags). For example, a downstream DAE solver that keeps factorizations across steps records the key in its generated modules' manifest and regenerates them when it changes.

Performance: the comparison record

The numbers below are one run of tools/compare_record.py at tree cb0617e on one machine (its json record, holding every sample, sha256 ff2cc551f5347012b8aa88870be4b88afd36834445dde17f4d923364f51bf1fe; the record is kept with the development tree). Machine: Intel(R) Core(TM) Ultra 7 258V, 8 logical CPUs, Windows power plan Balanced; Python 3.11.9, numpy 2.4.6, numba 0.67.0; the wheel's build (-O3 -DNDEBUG, MinGW gcc 14.2.0). Single-threaded (OPENBLAS / OMP / MKL / NUMEXPR threads = 1; Eigen with EIGEN_DONT_PARALLELIZE). The protocol: per arm the factorization is the second of two calls, 7 repeats; the substitution is solve_into on a held out from a numba loop, 10,000 calls a repeat x 7 repeats; refactorize 100 calls x 7; a repeat's calls cut to fit 2 s (the actual N in the record); an arm whose first factorization takes over 0.5 s is run as a function test only (its solution checked, no timing). Every value is the median of the 7 repeats, in microseconds; min and max are in the record. Every arm's solution is within 1e-10 (relative, max norm) of the 0.1.11 SuperLU solution; SuperLU NATURAL's is bit-identical to it at every size (the same CPU and BLAS kernel).

The matrices: G, the unsymmetric chain matrix of a power-system model (n = 688 / 2068 / 6898 / 20698); S, its symmetric positive definite nodal reduction (n = 207 / 621 / 2070 / 6210). The same right-hand side every call.

Substitution (us per solve_into) against n

Factorization and refactorization (us) against n

Each cell: factorize / refactorize (--: 0.1.11 has no refactorize). UMFPACK is the user's library (6.3.5, SuiteSparse built by MSYS2), loaded at run time.

The wrapper's overhead against n: sparse_numba / a C++ direct call

A C++ driver (tools/cpp_direct/direct_arms.cpp) calls the same SuperLU and Eigen builds directly -- the same DLLs, the extensions' own compile flags, the shims' options and template parameters, no sparse_numba in between -- in the same process, on the same matrix and right-hand side, under the same protocol. Each cell is the sparse_numba median divided by the direct-call median (1.000 = the same time). Single runs: the spread between repeats is in the record's min / max.

Column headers: the vehicle, then the n of G (the .G rows) and of S (the .S rows).

The direct-call arms' own medians, the substitution (us):

Function-only by the 0.5 s rule: EIGEN_LU at AMD on G at n 6898 and 20698 (Eigen's AMDOrdering on the symmetrized pattern of the unsymmetric G; first factorizations 2.7 s and 39.5 s), on both sides.

Architecture

User code (@njit)
    |
    v
sparse_numba.api / the kept interface modules (njit, cache=True; numba
ExternalFunctions on stable LLVM symbol names)
    |
    v
Cython layer (cdef api entry points)
    |
    v
the shims: superlu_wrapper.c, umfpack_wrapper.c (the user's library loaded
at run time), eigen_wrapper.cpp; one handle registry per backend
    |
    v
SuperLU 7.0.0 + OpenBLAS (vendored on Windows) / Eigen 3.4.0 (compiled
in) / the user's UMFPACK

Every layer releases the GIL.

Licences

sparse_numba is BSD 3-Clause. The components it ships or compiles in (docs/LICENSE_AUDIT.md holds each one's version, licence file and sha256):

On macOS SuperLU and OpenBLAS are Homebrew's, loaded at run time; the macOS wheels bundle nothing. CI reads every wheel's own listing (tests/test_sn6_wheel.py): each library a Linux wheel bundles must have its licence text in the wheel.

Citation

@software{hong2025sparse_numba,
  author = {Hong, Tianqi},
  title = {Sparse_Numba: A Numba-Compatible Sparse Solver},
  year = {2025},
  publisher = {GitHub},
  url = {https://github.com/th1275/sparse_numba}
}

Release process

1. Set the version in setup.py, setup.cfg and sparse_numba/__init__.py (cache_key() reads sparse_numba.__version__).
2. Open a pull request to main: CI builds one wheel per Python (cp310 - cp313) on Windows, manylinux x86_64 and macOS arm64 and tests each on its own Python (on Windows the full test suite from the wheel, with no UMFPACK library on the runner, and the wheel's contents; on Linux and macOS the import check and the wheel's contents); no wheel may carry a SuiteSparse library; no upload.
3. When CI is green, merge, tag the merge commit vX.Y.Z and push the tag: the upload_pypi job builds from the tagged tree and publishes the wheels and the sdist (trusted publishing).
4. A failed tag push before upload_pypi ran: delete the tag, fix, re-tag. A version uploaded to PyPI is never re-uploaded: bump the version.

Contributing

Contributions are welcome. Please open an issue or pull request on GitHub.