Algebraic Multigrid (AMG) solvers for large sparse linear systems with Python interface
npx @tessl/cli install tessl/pypi-pyamg@5.3.0PyAMG is a comprehensive library of Algebraic Multigrid (AMG) solvers with a convenient Python interface. It provides optimal or near-optimal efficiency for solving large sparse linear systems, particularly excelling with problems discretized on unstructured meshes and irregular grids where geometric information is limited or unavailable.
The library features implementations of Ruge-Stuben (Classical AMG), Smoothed Aggregation (SA), Adaptive Smoothed Aggregation (αSA), and Compatible Relaxation (CR) methods, combining Python accessibility with performance-critical C++ extensions.
pip install pyamgimport pyamgCommon usage imports:
import pyamg
import numpy as np
import scipy.sparse as spimport pyamg
import numpy as np
# Create a 2D Poisson problem
A = pyamg.gallery.poisson((100, 100), format='csr')
# Construct multigrid hierarchy using Classical AMG
ml = pyamg.ruge_stuben_solver(A)
# Create random right-hand side vector
b = np.random.rand(A.shape[0])
# Solve the linear system Ax = b
x = ml.solve(b, tol=1e-10)
# Check residual
residual = np.linalg.norm(b - A*x)
print(f"Residual: {residual}")
# Alternative: use high-level solve interface
x = pyamg.solve(A, b)PyAMG is built around the MultilevelSolver class, which manages a hierarchy of progressively coarser grids and operators. The architecture enables efficient solution of large linear systems through:
Key components:
Simple interface for solving linear systems without detailed AMG knowledge. Automatically selects appropriate solver methods and parameters based on problem characteristics.
def solve(A, b, x0=None, tol=1e-5, maxiter=400, return_solver=False,
existing_solver=None, verb=True, residuals=None):
"""
Solve the linear system Ax = b using AMG.
Parameters:
- A: sparse matrix, coefficient matrix
- b: array, right-hand side vector
- x0: array, initial guess
- tol: float, convergence tolerance (default 1e-5)
- maxiter: int, maximum iterations (default 400)
- return_solver: bool, return solver if True
- existing_solver: MultilevelSolver, reuse existing solver
- verb: bool, verbose output if True
- residuals: list, store residual norms
Returns:
array or tuple: solution vector x, optionally with solver
"""
def solver(A, config):
"""
Create AMG solver for matrix A using configuration.
Parameters:
- A: sparse matrix, coefficient matrix
- config: dict, solver configuration dictionary
Returns:
MultilevelSolver: configured AMG solver
"""
def solver_configuration(A, B=None, verb=True):
"""
Generate solver configuration dictionary.
Parameters:
- A: sparse matrix, coefficient matrix
- B: array, near null-space modes (default None)
- verb: bool, verbose output if True
Returns:
dict: solver configuration parameters
"""Factory functions for creating specific types of AMG solvers with detailed control over algorithm parameters and behavior.
def ruge_stuben_solver(A, strength=('classical', {'theta': 0.25}),
CF=('RS', {'second_pass': False}),
interpolation='classical',
presmoother=('gauss_seidel', {'sweep': 'symmetric'}),
postsmoother=('gauss_seidel', {'sweep': 'symmetric'}),
max_levels=30, max_coarse=10, keep=False, **kwargs):
"""
Create Classical (Ruge-Stuben) AMG solver.
Parameters:
- A: sparse matrix, coefficient matrix
- strength: str or tuple, strength of connection measure
- CF: str or tuple, coarse/fine splitting method
- interpolation: str, interpolation method
- presmoother: str or tuple, pre-smoothing method
- postsmoother: str or tuple, post-smoothing method
- max_levels: int, maximum number of levels
- max_coarse: int, maximum coarse grid size
- keep: bool, keep extra operators for debugging
Returns:
MultilevelSolver: Classical AMG solver
"""
def smoothed_aggregation_solver(A, B=None, BH=None, symmetry='hermitian',
strength='symmetric', aggregate='standard',
smooth=('jacobi', {'omega': 4.0/3.0}),
presmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
postsmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
improve_candidates=(('block_gauss_seidel',
{'sweep': 'symmetric', 'iterations': 4}), None),
max_levels=10, max_coarse=10, diagonal_dominance=False,
keep=False, **kwargs):
"""
Create Smoothed Aggregation AMG solver.
Parameters:
- A: sparse matrix, coefficient matrix
- B: array, near null-space modes
- BH: array, hermitian transpose of B
- symmetry: str, matrix symmetry assumption
- strength: str, strength of connection measure
- aggregate: str, aggregation method
- smooth: tuple, prolongation smoothing method
- presmoother: tuple, pre-smoothing method
- postsmoother: tuple, post-smoothing method
- improve_candidates: tuple, candidate improvement parameters
- max_levels: int, maximum number of levels
- max_coarse: int, maximum coarse grid size
- diagonal_dominance: bool, diagonal dominance flag
- keep: bool, keep extra operators for debugging
Returns:
MultilevelSolver: Smoothed Aggregation AMG solver
"""
def air_solver(A, strength='classical', presmoother='l1_gauss_seidel',
postsmoother='l1_gauss_seidel', max_levels=10, max_coarse=500,
coarse_solver='pinv2', **kwargs):
"""
Create Approximate Ideal Restriction AMG solver.
Parameters:
- A: sparse matrix, coefficient matrix
- strength: str, strength of connection measure
- presmoother: str, pre-smoothing method
- postsmoother: str, post-smoothing method
- max_levels: int, maximum number of levels
- max_coarse: int, maximum coarse grid size
- coarse_solver: str, coarse grid solver method
Returns:
MultilevelSolver: AIR AMG solver
"""
def adaptive_sa_solver(A, initial_candidates=None, symmetry='hermitian',
num_candidates=1, candidate_iters=5, epsilon=0.1,
max_levels=10, max_coarse=10, **kwargs):
"""
Create Adaptive Smoothed Aggregation AMG solver.
Parameters:
- A: sparse matrix, coefficient matrix
- initial_candidates: array, initial candidate basis
- symmetry: str, matrix symmetry assumption
- num_candidates: int, number of candidates to generate
- candidate_iters: int, smoothing iterations per candidate
- epsilon: float, target convergence factor
- max_levels: int, maximum number of levels
- max_coarse: int, maximum coarse grid size
Returns:
MultilevelSolver: Adaptive SA AMG solver
"""The core solver class that manages AMG hierarchies and provides solve methods with extensive control over solution process.
class MultilevelSolver:
"""
Multilevel AMG solver managing hierarchy of operators and grids.
"""
def solve(self, b, x0=None, tol=1e-5, maxiter=None,
cycle='V', accel=None, **kwargs):
"""
Solve Ax = b using multigrid cycles.
Parameters:
- b: array, right-hand side vector
- x0: array, initial guess
- tol: float, convergence tolerance
- maxiter: int, maximum iterations
- cycle: str, multigrid cycle type ('V', 'W', 'F')
- accel: str, acceleration method
Returns:
array: solution vector
"""
def aspreconditioner(self, cycle='V'):
"""
Return linear operator for use as preconditioner.
Parameters:
- cycle: str, multigrid cycle type
Returns:
LinearOperator: preconditioner operator
"""Comprehensive collection of test matrices and problems for AMG development and evaluation, including PDE discretizations and finite element problems.
def poisson(grid, format='csr', dtype=float):
"""
Generate discrete Poisson operator.
Parameters:
- grid: tuple, grid dimensions (nx,) or (nx, ny) or (nx, ny, nz)
- format: str, sparse matrix format
- dtype: data type
Returns:
sparse matrix: Poisson operator
"""
def linear_elasticity(grid, format='csr', dtype=float):
"""
Generate linear elasticity operator.
Parameters:
- grid: tuple, grid dimensions
- format: str, sparse matrix format
- dtype: data type
Returns:
sparse matrix: elasticity operator
"""
def demo(**kwargs):
"""
Run PyAMG demonstration showing solver capabilities.
"""Complete suite of Krylov subspace methods for iterative solution of linear systems, optimized for use with AMG preconditioning.
def gmres(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
"""
Generalized Minimal Residual method.
Parameters:
- A: sparse matrix or LinearOperator
- b: array, right-hand side
- x0: array, initial guess
- tol: float, convergence tolerance
- maxiter: int, maximum iterations
- M: preconditioner
Returns:
tuple: (solution, info)
"""
def cg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
"""
Conjugate Gradient method.
Parameters:
- A: sparse matrix or LinearOperator
- b: array, right-hand side
- x0: array, initial guess
- tol: float, convergence tolerance
- maxiter: int, maximum iterations
- M: preconditioner
Returns:
tuple: (solution, info)
"""Methods for determining the strength of connections between unknowns, fundamental to AMG coarsening strategies.
def classical_strength_of_connection(A, theta=0.25):
"""
Classical strength of connection measure.
Parameters:
- A: sparse matrix, coefficient matrix
- theta: float, strength threshold
Returns:
sparse matrix: strength matrix
"""
def symmetric_strength_of_connection(A, theta=0.0):
"""
Symmetric strength of connection measure.
Parameters:
- A: sparse matrix, coefficient matrix
- theta: float, strength threshold
Returns:
sparse matrix: strength matrix
"""
def energy_based_strength_of_connection(A, theta=0.0, **kwargs):
"""
Energy-based strength of connection measure.
Parameters:
- A: sparse matrix, coefficient matrix
- theta: float, strength threshold
Returns:
sparse matrix: strength matrix
"""Algorithms for grouping fine grid points into coarse grid aggregates for Smoothed Aggregation AMG.
def standard_aggregation(C, **kwargs):
"""
Standard aggregation algorithm.
Parameters:
- C: sparse matrix, strength of connection matrix
Returns:
array: aggregation array mapping fine to coarse
"""
def naive_aggregation(C, **kwargs):
"""
Naive aggregation algorithm.
Parameters:
- C: sparse matrix, strength of connection matrix
Returns:
array: aggregation array
"""
def lloyd_aggregation(C, ratio=0.03, distance='unit', **kwargs):
"""
Lloyd clustering-based aggregation.
Parameters:
- C: sparse matrix, strength of connection matrix
- ratio: float, desired coarsening ratio
- distance: str, distance measure
Returns:
array: aggregation array
"""Smoothing and relaxation methods for error correction in multigrid cycles.
def gauss_seidel(A, x, b, iterations=1, sweep='forward'):
"""
Gauss-Seidel relaxation.
Parameters:
- A: sparse matrix, coefficient matrix
- x: array, solution vector (modified in-place)
- b: array, right-hand side
- iterations: int, number of relaxation sweeps
- sweep: str, sweep direction ('forward', 'backward', 'symmetric')
"""
def jacobi(A, x, b, iterations=1, omega=1.0):
"""
Jacobi relaxation.
Parameters:
- A: sparse matrix, coefficient matrix
- x: array, solution vector (modified in-place)
- b: array, right-hand side
- iterations: int, number of relaxation sweeps
- omega: float, relaxation parameter
"""Helper functions for matrix operations, norms, and linear algebra utilities supporting AMG algorithms.
def norm(x, ord=None):
"""
Compute vector or matrix norm.
Parameters:
- x: array-like, input vector or matrix
- ord: norm type
Returns:
float: computed norm
"""
def approximate_spectral_radius(A, tol=0.1, maxiter=10):
"""
Estimate spectral radius of matrix.
Parameters:
- A: sparse matrix
- tol: float, tolerance for estimation
- maxiter: int, maximum iterations
Returns:
float: estimated spectral radius
"""
def make_system(A, b, x0=None, formats=None):
"""
Create consistent linear system with proper formats.
Parameters:
- A: matrix, coefficient matrix
- b: array, right-hand side
- x0: array, initial guess
- formats: dict, desired matrix formats
Returns:
tuple: (A, b, x0) in consistent formats
"""# Type aliases for common PyAMG objects
MultilevelSolver = pyamg.multilevel.MultilevelSolverPyAMG raises standard Python exceptions:
Common error patterns:
.tocsr() for CSR format)