Solve sparse Laplacian and SDDM systems with reusable approximate-Cholesky factors. NumPy is required. Install the SciPy extra for sparse-matrix helpers.
python -m pip install "apxchol[scipy]"For an unreleased checkout, use the source installation below; pip installs the published version.
Wheels work on CPU without a GPU. Compatible hardware enables the included GPU backend; no CUDA Toolkit or Homebrew installation is needed.
| Platform | CPU requirement | Optional GPU backend |
|---|---|---|
| Linux | x86-64-v2 or newer | NVIDIA CUDA |
| macOS 15+ | Apple silicon (arm64) | Apple Metal |
Wheels support CPython 3.11–3.15 and PyPy 7.3.23 (Python 3.11, pp73).
On PyPy, install apxchol and use the NumPy-only interface.
import apxchol
import numpy as np
from scipy.sparse import diags
n = 1000
# Tridiagonal matrix: 3 on the diagonal, -1 on each neighboring diagonal.
A = diags([-1.0, 3.0, -1.0], [-1, 0, 1], shape=(n, n), format="csc")
b = np.ones(n)
solver = apxchol.factorize(A)
result = solver.solve(b, rtol=1e-8)
assert result.converged
x = result.xFor one solve, use apxchol.solve(A, b). Results include x, iterations,
residual and converged. Convergence is always checked against the original
A and b; inspect converged for success.
See the shared binding capabilities for comparison with MATLAB/Octave.
Pass the matrix A in A @ x = b. For a graph's adjacency matrix, first use
A = apxchol.as_operator(Adj): it returns a valid operator unchanged and
otherwise forms the graph Laplacian. See supported matrices, input
conversion and RHS requirements.
Any SciPy sparse format is converted to sorted, duplicate-free CSC without changing the input. Supply the full symmetric matrix, including both entries of each off-diagonal pair. Complex inputs are rejected, even with zero imaginary parts.
python -m pip install apxchol installs the solver without SciPy. Supply the
three zero-based CSC arrays directly; all usual solver options still apply:
solver = apxchol.Solver.from_csc(data, indices, indptr, backend="auto")
result = solver.solve(b)The matrix is square with len(indptr) - 1 rows. Integer indices may be 32- or
64-bit; the native importer sorts columns and sums duplicate entries. SciPy
is needed for factorize(A), as_operator, aslinearoperator, and sparse factor
exports. NumPy-only solves and preconditioner application do not import it.
Subsequent solves reuse the factor and workspace:
first = solver.solve(b)
b2 = np.arange(n, dtype=float)
B = np.column_stack((b, b2))
batch = solver.solve(B) # native batched solve, one RHS per column
assert batch.converged.all()Vector inputs have shape (n,); batches have shape (n, nrhs). result.x
and apply outputs preserve that shape and use float64. Batch convergence
fields contain one value per column.
solver.setup_seconds is the wall-clock time the factor took to build, and
result.solve_seconds that of one solve call, for all columns of a batch.
Over k separate one-RHS calls each costs setup_seconds / k + solve_seconds
on average, and one batched call of k columns costs
(setup_seconds + solve_seconds) / k per RHS; see judging whether reuse pays off.
A warm start supplies an initial solution estimate:
result = solver.solve(b2, x0=first.x)| Solve argument | Meaning |
|---|---|
rtol |
Positive, finite relative residual tolerance; omitted uses the native default |
maxiter |
Nonnegative iteration budget; omitted uses the native default |
x0 |
Initial guess with the same shape as the RHS; omitted starts from zero |
out |
Writable, C- or F-contiguous float64 array with the RHS shape; returned as result.x |
The native reference shows the compiled tolerance and iteration defaults.
Overlapping inputs are copied before writing out. Factorization, solve
and apply release the GIL while they run native code. Reading the input
matrix holds it. A call made while another thread uses the same solver raises
RuntimeError. Give each thread copy.copy(solver). The copy shares the factor
and allocates only its own workspace. On a GPU it shares the device factor
and solves on its own CUDA stream or Metal command queue.
Save a reusable solver's operator and factor, then choose the solve backend when loading. These methods also work without SciPy:
solver.save("system.apxchol")
restored = apxchol.Solver.load("system.apxchol", backend="cpu")
result = restored.solve(b)Saving replaces an existing file, writing to <path>.tmp first and renaming it
into place. The archive
stores the original operator and the installed factor, the preconditioner the
solves actually use, in its storage dtype; RHSs, initial guesses and device
workspace are not saved. Loading reinstalls that factor exactly without
repeating factorization, so solves match the saved solver bit for bit on the
same backend and thread count.
backend="auto" adopts on CPU; backend="gpu" uploads to a compatible GPU.
The file fixes the storage dtype and dropping, so factor_storage and
drop_tolerance have no effect when loading. Pass lump_max_mass again if the
saved operator requires positive-entry lumping.
Pass these keywords to factorize or the one-shot solve:
| Keyword | Meaning |
|---|---|
backend |
"auto" chooses a compatible CPU/GPU route; "cpu" or "gpu" requests one explicitly |
num_threads |
CPU thread ceiling; 0 inherits the caller's OpenMP setting |
seed |
Sampling seed; None uses the native default |
partitioner |
"block_greedy", "priority_greedy" or "baumann_kyng"; None uses the native default |
factor_storage |
"auto" (default) for the backend default, "float16" or "float32"; NumPy dtypes also work |
Read the chosen route from solver.backend and storage type from
solver.factor_storage, a NumPy dtype. See backend and thread selection
for automatic selection and explicit-request behavior.
Storage defaults to FP32 on CPU/Metal and FP16 on CUDA. This affects triangular-solve storage and the exported factor, not outer solve precision. Unsupported FP16 requests fail; FP64 storage is unavailable. See precision.
These share the native defaults and can usually be left unset. GPU routes support fewer construction choices.
| Keyword | Meaning |
|---|---|
sampler |
Clique sampler: "gks" (default) or "trace_cycle" (CPU and CUDA only) |
degree_quantile |
Degree cutoff quantile in [0, 1); 0 selects the multiplier-based cutoff (Metal requires (0, 1)); None uses the backend default |
degree_multiplier |
Positive degree multiplier, used only when degree_quantile=0 |
drop_tolerance |
Relative factor-entry threshold in [0, 1]; 0 disables dropping |
regularization |
Minimum diagonal excess of the preconditioner, as a fraction of each diagonal after lumping; excess below 1e-12 of the diagonal counts as zero first; A is not modified; 0 disables it |
lump_max_mass |
Allowed positive off-diagonal mass fraction in [0, 1]; 0 rejects such entries |
verbose |
Print setup diagnostics; off by default |
Unknown keywords raise ValueError, like other invalid input or options;
execution failures raise RuntimeError (see
results and errors).
The Python reference
and native option reference describe the controls.
Use the built-in solver.solve for the standard PCG path. To use apxchol as a
preconditioner for another SciPy iterative method:
from scipy.sparse.linalg import cg
M = solver.aslinearoperator()
x, info = cg(A, b, M=M, rtol=1e-8)SciPy controls convergence. Use this when you need its algorithm or interface. With a GPU factor, each application exchanges vectors with the host; the built-in GPU solve keeps the iteration on the device.
factorize retains exportable arrays by default. Exported factor values are
FP32, regardless of backend or internal storage:
P = solver.P # P[original vertex] = elimination position
G = solver.chol() # lower CSC factor, including its diagonal
p = P.argsort()
# A[p][:, p] ≈ G @ G.Tchol(), P and factor_nnz describe the preconditioner actually used: the
factor after dropping and storage rounding. chol() returns a float64 copy, P
an int64 permutation and factor_nnz an integer. For a Laplacian, they
represent the nonsingular subspace and the last row and column
are empty; see the factor reference.
A GPU-backed solver downloads the factor on the first export and caches it.
Matrix import is limited to signed 32-bit dimensions and stored-entry counts. The native library's wide-index build options do not enlarge this binding.
From a repository checkout, install a development build and run its tests:
python -m pip install -e './python[test,scipy]'
python -m pytest python/tests -vThe -e option links the installed package to the local Python source.
Re-run the install command after changing C++ code. Local builds tune for the
build machine; see build instructions.
Set CMAKE_ARGS="-DAPXCHOL_DEVICE=CUDA" or CMAKE_ARGS="-DAPXCHOL_DEVICE=METAL"
for a GPU source build, with its native dependencies installed. Without build
isolation, install pybind11 scikit-build-core first. License.