Documentation

AMJax rebuilds PyAMG's algebraic multigrid solvers in JAX. The JAX compiler has rules that make some features easy to add, others possible with effort, and some not possible at all. This table shows what is currently available.

JAX constraints

  1. Control flow must be known at compile time: loops and conditionals depending on data cannot be expressed freely.
  2. Data is immutable: every array modification creates a new array; sequential in-place updates cannot be optimized by the compiler.
  3. No side effects inside compiled functions: no external writes, no print statements, no arbitrary Python calls.

Feature Implemented Implementable Justification
Smoothers
JacobiYes
Gauss-SeidelNoNo
  • Update of x_i depends on x_j^{k+1} for j < i: irreducible sequential dependency
  • lax.scan would produce an O(n) XLA graph, not vectorizable on GPU
SORNoNo
  • Gauss-Seidel with omega: same sequential dependency
  • Same lax.scan limitation
Block Gauss-SeidelNoNo
  • Same structure as GS by blocks, sequential inter-block dependency
Block JacobiNoYes
  • Blocks are independent: all can be solved simultaneously
  • Block inverses precomputed at setup
Schwarz (multiplicative)NoNo
  • Each local solve modifies x, creating a sequential intra-sweep dependency
Polynomial (generic)NoYes
  • General polynomial smoother: p(A)r with static degree
  • Richardson and Chebyshev are special cases (see below)
Richardson (via polynomial)NoYes
  • One matrix-vector product and vector additions per iteration: fully parallel
  • omega precomputed at setup
Chebyshev (via polynomial)NoYes
  • Horner's rule with static degree: unrolled at tracing, each matmul parallel
  • Coefficients precomputed from spectrum
Jacobi-NENoYes
  • Same structure as Jacobi, fully parallel
  • D_{AA^H} precomputed at setup
Gauss-Seidel-NENoNo
  • GS on normal equations AA^Hy = b: sequential row by row
Gauss-Seidel-NRNoNo
  • GS on A^HAx = A^Hb: sequential column by column
Gauss-Seidel indexedNoPartial
  • Implementable if the index sets are independent (e.g. graph coloring)
  • Otherwise same sequential dependency as GS
CF-Jacobi / FC-JacobiNoYes
  • Jacobi on C/F-points separately
  • Within each group, fully independent via jnp.where on static mask
CG as smootherNoYes
  • jax.scipy.sparse.linalg.cg with static maxiter is jit-compilable
GMRES as smootherNoYes
  • jax.scipy.sparse.linalg.gmres with static maxiter is jit-compilable
CGNE / CGNR as smootherNoYes (effort)
  • Absent from jax.scipy
  • Requires custom implementation (algorithmically simple: CG on normal equations)
Coarse solvers
JacobiYes
pinvYes
LUYes
QRYes
CholeskyNoYes
  • jax.scipy.linalg.cho_factor/cho_solve exist
  • Same pattern as LU already implemented
  • Requires SPD matrix
pinv2NoYes
  • pinv with rcond parameter to control SVD truncation
  • jnp.linalg.pinv already supports rcond
splu (sparse LU)NoNo
  • Relies on SuperLU (external C library): no native sparse LU in XLA
  • jax.pure_callback possible but breaks vmap and grad
CGNoYes
  • jax.scipy.sparse.linalg.cg with static maxiter
GMRESNoYes
  • jax.scipy.sparse.linalg.gmres with static max_iter
BiCGSTABNoYes
  • jax.scipy.sparse.linalg.bicgstab with static maxiter
FGMRESNoYes (very high effort)
  • Requires variable preconditioner per iteration
  • Not supported by jax.scipy
Cycles
V-cycleYes
W-cycleYes
F-cycleYes
AMLINoPartial
  • Cycle itself is implementable
  • Requires FGMRES as accelerator, not available in JAX
Krylov acceleration
CGYes
  • jax.scipy.sparse.linalg.cg(A, b, M=ml.aspreconditioner())
GMRESYes
  • jax.scipy.sparse.linalg.gmres(A, b, M=ml.aspreconditioner())
BiCGSTABYes
  • jax.scipy.sparse.linalg.bicgstab(A, b, M=ml.aspreconditioner())
FGMRESNoYes (very high effort)
  • Absent from jax.scipy
  • Requires variable preconditioner per Krylov iteration
AMG hierarchies
Ruge-StubenYes
Smoothed AggregationYes
Root NodeYes
Adaptive SANoYes (effort)
  • Adaptive process to find near-null components
  • Implementable with a fixed number of iterations unrolled at trace time
AIR (Approx. Ideal Restriction)NoYes (effort)
  • Requires solving local least-squares problems at setup
  • Solvable with jnp.linalg.lstsq if subdomain sizes are static
Other
Complex matricesNoYes
  • JAX supports complex64/complex128 natively
  • Requires adjustments to the current implementation
Residual tracking under jitNoNo
  • lax.while_loop forbids side effects: list.append impossible inside XLA graph
  • Debug mode only via jax.debug.callback, incompatible with vmap and grad