tessl/pypi-pyamg

Algebraic Multigrid (AMG) solvers for large sparse linear systems with Python interface

—

Pending

Overview

Eval results

Files

Aggregation Methods

Name: tessl/pypi-pyamg
Author: tessl

Algorithms for grouping fine grid points into coarse grid aggregates for Smoothed Aggregation AMG. Aggregation determines the coarse grid structure and significantly impacts solver performance.

Capabilities

Standard Aggregation Methods

Core aggregation algorithms for general problems.

def standard_aggregation(C, **kwargs):
    """
    Standard aggregation algorithm.
    
    Greedy algorithm that creates aggregates by selecting seed points
    and adding strongly connected neighbors. Most commonly used
    aggregation method for Smoothed Aggregation AMG.
    
    Parameters:
    - C: sparse matrix, strength of connection matrix
    - **kwargs: additional parameters:
        * symmetrize: bool, symmetrize strength matrix (default True)
        * aggregate_mat: bool, return aggregation matrix (default False)
    
    Returns:
    array: aggregation array where AggOp[i] = j means
        fine point i belongs to aggregate j
    
    Notes:
    - Creates reasonably sized aggregates with good connectivity
    - Balances aggregate size and quality
    - Default choice for most problems
    """

def naive_aggregation(C, **kwargs):
    """
    Naive aggregation algorithm.
    
    Simple aggregation that pairs neighboring points without
    sophisticated connectivity analysis. Fast but may create
    poor aggregate quality.
    
    Parameters:
    - C: sparse matrix, strength of connection matrix
    - **kwargs: additional aggregation parameters
    
    Returns:
    array: aggregation array
    
    Notes:
    - Fastest aggregation method
    - May create disconnected or poor-quality aggregates
    - Useful for testing and comparison purposes
    """

Usage Examples:

import pyamg
import numpy as np

# Standard aggregation for Poisson problem
A = pyamg.gallery.poisson((50, 50))
C = pyamg.strength.symmetric_strength_of_connection(A)
AggOp = pyamg.aggregation.standard_aggregation(C)
print(f"Created {AggOp.max() + 1} aggregates from {len(AggOp)} points")

# Naive aggregation for comparison
AggOp_naive = pyamg.aggregation.naive_aggregation(C) 
print(f"Naive: {AggOp_naive.max() + 1} aggregates")

# Use in solver construction
ml = pyamg.smoothed_aggregation_solver(A, aggregate='standard')

Advanced Aggregation Methods

Sophisticated aggregation algorithms for challenging problems.

def lloyd_aggregation(C, ratio=0.03, distance='unit', maxiter=10, **kwargs):
    """
    Lloyd clustering-based aggregation.
    
    Uses Lloyd clustering algorithm (k-means variant) to create
    aggregates with optimal geometric properties. Particularly
    effective for problems with geometric structure.
    
    Parameters:
    - C: sparse matrix, strength of connection matrix
    - ratio: float, desired coarsening ratio (0 < ratio < 1):
        * 0.03: aggressive coarsening (default)
        * 0.1: moderate coarsening
        * 0.2: mild coarsening
    - distance: str, distance measure for clustering:
        * 'unit': unit distance (default)
        * 'abs': absolute value distance  
        * 'min': minimum distance
    - maxiter: int, maximum Lloyd iterations
    - **kwargs: additional clustering parameters
    
    Returns:
    array: aggregation array with improved geometric properties
    
    Notes:
    - More expensive than standard aggregation
    - Can create more uniform aggregate sizes
    - Requires coordinate information for best results
    """

def balanced_lloyd_aggregation(C, ratio=0.03, **kwargs):
    """
    Balanced Lloyd clustering aggregation.
    
    Variant of Lloyd aggregation that enforces more uniform
    aggregate sizes while maintaining clustering quality.
    
    Parameters:
    - C: sparse matrix, strength of connection matrix
    - ratio: float, desired coarsening ratio
    - **kwargs: Lloyd clustering parameters
    
    Returns:
    array: aggregation array with balanced aggregate sizes
    
    Notes:
    - Attempts to create aggregates of similar size
    - May sacrifice some geometric optimality for balance
    - Good for problems requiring uniform coarsening
    """

def metis_aggregation(C, ratio=0.03, **kwargs):
    """
    METIS-based graph partitioning aggregation.
    
    Uses METIS graph partitioning library to create aggregates
    by partitioning the strength graph. Produces high-quality
    aggregates but requires METIS installation.
    
    Parameters:
    - C: sparse matrix, strength of connection matrix
    - ratio: float, desired coarsening ratio
    - **kwargs: METIS-specific parameters
    
    Returns:
    array: aggregation array from graph partitioning
    
    Notes:
    - Requires pymetis package installation
    - Creates high-quality graph partitions
    - More expensive than simpler methods
    - Excellent for unstructured problems
    
    Raises:
    ImportError: if METIS is not available
    """

Usage Examples:

# Lloyd aggregation with custom ratio
A = pyamg.gallery.linear_elasticity((30, 30))
C = pyamg.strength.symmetric_strength_of_connection(A)
AggOp_lloyd = pyamg.aggregation.lloyd_aggregation(C, ratio=0.05)

# Balanced aggregation
AggOp_balanced = pyamg.aggregation.balanced_lloyd_aggregation(C, ratio=0.05)

# METIS aggregation (if available)
try:
    AggOp_metis = pyamg.aggregation.metis_aggregation(C, ratio=0.05)
    print(f"METIS: {AggOp_metis.max() + 1} aggregates")
except ImportError:
    print("METIS not available")

# Compare aggregate counts
print(f"Lloyd: {AggOp_lloyd.max() + 1} aggregates")
print(f"Balanced: {AggOp_balanced.max() + 1} aggregates")

Prolongation Smoothing

Methods for smoothing tentative prolongation operators created from aggregation.

def jacobi_prolongation_smoother(A, T, C=None, omega=1.0):
    """
    Jacobi prolongation smoother.
    
    Applies weighted Jacobi smoothing to tentative prolongation
    operator to improve interpolation quality. Most common
    prolongation smoother.
    
    Parameters:
    - A: sparse matrix, fine level operator
    - T: sparse matrix, tentative prolongation operator
    - C: sparse matrix, strength of connection (optional)
    - omega: float, relaxation parameter (0 < omega <= 1):
        * 1.0: full Jacobi smoothing (default)
        * 4/3: common choice for SPD problems
        * 0.67: conservative smoothing
    
    Returns:
    sparse matrix: smoothed prolongation operator P
    
    Notes:
    - P = (I - omega * D^{-1} * A) * T where D = diag(A)
    - Improves interpolation quality at modest cost
    - Default choice for most problems
    """

def richardson_prolongation_smoother(A, T, omega=1.0, **kwargs):
    """
    Richardson prolongation smoother.
    
    Applies Richardson iteration to smooth tentative prolongation.
    Alternative to Jacobi with different stability properties.
    
    Parameters:
    - A: sparse matrix, fine level operator
    - T: sparse matrix, tentative prolongation
    - omega: float, relaxation parameter
    - **kwargs: additional smoothing parameters
    
    Returns:
    sparse matrix: Richardson-smoothed prolongation operator
    
    Notes:
    - Less common than Jacobi smoothing
    - May have different convergence properties
    - Experimental alternative smoother
    """

def energy_prolongation_smoother(A, T, **kwargs):
    """
    Energy-minimizing prolongation smoother.
    
    Computes prolongation that minimizes energy norm of
    interpolation error. Most sophisticated smoother but
    more expensive to compute.
    
    Parameters:
    - A: sparse matrix, fine level operator
    - T: sparse matrix, tentative prolongation
    - **kwargs: energy minimization parameters
    
    Returns:
    sparse matrix: energy-minimized prolongation operator
    
    Notes:
    - Optimal in energy norm sense
    - Higher computational cost than Jacobi
    - May provide better convergence for difficult problems
    """

Usage Examples:

# Manual prolongation smoothing
A = pyamg.gallery.poisson((40, 40))
C = pyamg.strength.symmetric_strength_of_connection(A)
AggOp = pyamg.aggregation.standard_aggregation(C)

# Create tentative prolongation (usually done internally)
n_fine = A.shape[0]
n_coarse = AggOp.max() + 1
T = sp.csr_matrix((np.ones(n_fine), (np.arange(n_fine), AggOp)), 
                  shape=(n_fine, n_coarse))

# Apply different smoothers
P_jacobi = pyamg.aggregation.jacobi_prolongation_smoother(A, T, omega=4.0/3.0)
P_energy = pyamg.aggregation.energy_prolongation_smoother(A, T)

print(f"Tentative nnz: {T.nnz}")
print(f"Jacobi nnz: {P_jacobi.nnz}")  
print(f"Energy nnz: {P_energy.nnz}")

# Use in solver construction
ml = pyamg.smoothed_aggregation_solver(A, smooth='jacobi')
ml_energy = pyamg.smoothed_aggregation_solver(A, smooth='energy')

Near-Nullspace Handling

Functions for incorporating near-nullspace information into aggregation.

def fit_candidates(AggOp, B, **kwargs):
    """
    Fit near-nullspace candidates to aggregation pattern.
    
    Creates tentative prolongation operator that preserves
    given near-nullspace vectors on each aggregate. Essential
    for problems with non-trivial nullspaces.
    
    Parameters:
    - AggOp: array, aggregation array mapping fine to coarse
    - B: array, near-nullspace vectors (n_fine × n_candidates):
        * Constant vector for scalar problems
        * Rigid body modes for elasticity
        * Custom modes for specific problems
    - **kwargs: fitting parameters
    
    Returns:
    sparse matrix: tentative prolongation operator that
        interpolates near-nullspace exactly
    
    Notes:
    - Critical for elasticity and other systems problems
    - B should span the near-nullspace of the operator
    - Quality of B significantly affects AMG performance
    """

Usage Examples:

# Near-nullspace for elasticity problem
A = pyamg.gallery.linear_elasticity((25, 25))
n_dof = A.shape[0]

# Create rigid body modes (3 modes for 2D elasticity)
B = np.zeros((n_dof, 3))
# Translation modes
B[0::2, 0] = 1.0  # x-translation
B[1::2, 1] = 1.0  # y-translation
# Rotation mode (simplified)
B[0::2, 2] = np.arange(0, n_dof//2)  # x-coordinates
B[1::2, 2] = -np.arange(0, n_dof//2)  # -y-coordinates

# Create aggregation and fit candidates
C = pyamg.strength.symmetric_strength_of_connection(A)
AggOp = pyamg.aggregation.standard_aggregation(C)
T = pyamg.aggregation.fit_candidates(AggOp, B)

# Use in solver
ml = pyamg.smoothed_aggregation_solver(A, B=B)

Aggregation Quality Analysis

Aggregate Statistics

def analyze_aggregation(AggOp, C=None):
    """Analyze aggregation quality."""
    n_fine = len(AggOp)
    n_coarse = AggOp.max() + 1
    
    # Basic statistics
    coarsening_ratio = n_coarse / n_fine
    avg_agg_size = n_fine / n_coarse
    
    # Size distribution
    agg_sizes = np.bincount(AggOp)
    min_size, max_size = agg_sizes.min(), agg_sizes.max()
    
    print(f"Coarsening ratio: {coarsening_ratio:.3f}")
    print(f"Average aggregate size: {avg_agg_size:.1f}")
    print(f"Aggregate size range: {min_size} - {max_size}")
    
    return {
        'coarsening_ratio': coarsening_ratio,
        'avg_size': avg_agg_size,
        'size_range': (min_size, max_size)
    }

# Example usage
A = pyamg.gallery.poisson((50, 50))
C = pyamg.strength.symmetric_strength_of_connection(A)
AggOp = pyamg.aggregation.standard_aggregation(C)
stats = analyze_aggregation(AggOp, C)

Selection Guidelines

Method Selection by Problem Type

General problems: Standard aggregation
Geometric problems: Lloyd or balanced Lloyd
Unstructured meshes: METIS aggregation
Testing/debugging: Naive aggregation

Parameter Tuning

# Aggressive coarsening (fewer levels, less memory)
AggOp_aggressive = pyamg.aggregation.lloyd_aggregation(C, ratio=0.02)

# Conservative coarsening (more levels, better convergence)
AggOp_conservative = pyamg.aggregation.lloyd_aggregation(C, ratio=0.1)