Pyamg Documentation

This page contains the Pyamg Package documentation.

The __config__ Module

pyamg.__config__.get_info(name)[source]
pyamg.__config__.show()[source]

The graph Module

Algorithms related to graphs

pyamg.graph.maximal_independent_set(G, algo='serial', k=None)[source]

Compute a maximal independent vertex set for a graph

G : sparse matrix
Symmetric matrix, preferably in sparse CSR or CSC format The nonzeros of G represent the edges of an undirected graph.
algo : {‘serial’, ‘parallel’}
Algorithm used to compute the MIS
  • serial : greedy serial algorithm
  • parallel : variant of Luby’s parallel MIS algorithm
An array S where
S[i] = 1 if vertex i is in the MIS S[i] = 0 otherwise

Diagonal entries in the G (self loops) will be ignored.

Luby’s algorithm is significantly more expensive than the greedy serial algorithm.

pyamg.graph.vertex_coloring(G, method='MIS')[source]

Compute a vertex coloring of a graph

G : sparse matrix
Symmetric matrix, preferably in sparse CSR or CSC format The nonzeros of G represent the edges of an undirected graph.
method : {string}
Algorithm used to compute the vertex coloring:
  • ‘MIS’ - Maximal Independent Set
  • ‘JP’ - Jones-Plassmann (parallel)
  • ‘LDF’ - Largest-Degree-First (parallel)

An array of vertex colors (integers beginning at 0)

Diagonal entries in the G (self loops) will be ignored.

pyamg.graph.bellman_ford(G, seeds, maxiter=None)[source]

Bellman-Ford iteration

CLR

pyamg.graph.lloyd_cluster(G, seeds, maxiter=10)[source]

Perform Lloyd clustering on graph with weighted edges

G : csr_matrix or csc_matrix
A sparse NxN matrix where each nonzero entry G[i,j] is the distance between nodes i and j.
seeds : {int, array}
If seeds is an integer, then its value determines the number of clusters. Otherwise, seeds is an array of unique integers between 0 and N-1 that will be used as the initial seeds for clustering.
maxiter : int
The maximum number of iterations to perform.

If G has complex values, abs(G) is used instead.

pyamg.graph.connected_components(G)[source]

Compute the connected components of a graph

The connected components of a graph G, which is represented by a symmetric sparse matrix, are labeled with the integers 0,1,..(K-1) where K is the number of components.

G : symmetric matrix, preferably in sparse CSR or CSC format
The nonzeros of G represent the edges of an undirected graph.
components : ndarray
An array of component labels for each vertex of the graph.

If the nonzero structure of G is not symmetric, then the result is undefined.

>>> from pyamg.graph import connected_components
>>> print connected_components( [[0,1,0],[1,0,1],[0,1,0]] )
[0 0 0]
>>> print connected_components( [[0,1,0],[1,0,0],[0,0,0]] )
[0 0 1]
>>> print connected_components( [[0,0,0],[0,0,0],[0,0,0]] )
[0 1 2]
>>> print connected_components( [[0,1,0,0],[1,0,0,0],[0,0,0,1],[0,0,1,0]] )
[0 0 1 1]

The multilevel Module

Generic AMG solver

class pyamg.multilevel.multilevel_solver(levels, coarse_solver='pinv2')[source]

Stores multigrid hierarchy and implements the multigrid cycle

The class constructs the cycling process and points to the methods for coarse grid solves. A multilevel_solver object is typically returned from a particular AMG method (see ruge_stuben_solver or smoothed_aggregation_solver for example). A call to multilevel_solver.solve() is a typical access point. The class also defines methods for constructing operator, cycle, and grid complexities.

levels : level array
Array of level objects that contain A, R, and P.
coarse_solver : string
String passed to coarse_grid_solver indicating the solve type
aspreconditioner()
Create a preconditioner using this multigrid cycle
cycle_complexity()
A measure of the cost of a single multigrid cycle.
grid_complexity()
A measure of the rate of coarsening.
operator_complexity()
A measure of the size of the multigrid hierarchy.
solve()
Iteratively solves a linear system for the right hand side.
aspreconditioner(cycle='V')[source]

Create a preconditioner using this multigrid cycle

cycle : {‘V’,’W’,’F’,’AMLI’}
Type of multigrid cycle to perform in each iteration.
precond : LinearOperator
Preconditioner suitable for the iterative solvers in defined in the scipy.sparse.linalg module (e.g. cg, gmres) and any other solver that uses the LinearOperator interface. Refer to the LinearOperator documentation in scipy.sparse.linalg

multilevel_solver.solve, scipy.sparse.linalg.LinearOperator

>>> from pyamg.aggregation import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> from scipy import rand
>>> A = poisson((100, 100), format='csr')          # matrix
>>> b = rand(A.shape[0])                           # random RHS
>>> ml = smoothed_aggregation_solver(A)            # AMG solver
>>> M = ml.aspreconditioner(cycle='V')             # preconditioner
>>> x, info = cg(A, b, tol=1e-8, maxiter=30, M=M)  # solve with CG
cycle_complexity(cycle='V')[source]
Cycle complexity of this multigrid hierarchy for V(1,1), W(1,1),
AMLI(1,1) and F(1,1) cycles when simple relaxation methods are used.

Cycle complexity is an approximate measure of the number of floating point operations (FLOPs) required to perform a single multigrid cycle relative to the cost a single smoothing operation.

cycle : {‘V’,’W’,’F’,’AMLI’}
Type of multigrid cycle to perform in each iteration.
cc : float
Defined as F_sum / F_0, where F_sum is the total number of nonzeros in the matrix on all levels encountered during a cycle and F_0 is the number of nonzeros in the matrix on the finest level.

This is only a rough estimate of the true cycle complexity. The estimate assumes that the cost of pre and post-smoothing are (each) equal to the number of nonzeros in the matrix on that level. This assumption holds for smoothers like Jacobi and Gauss-Seidel. However, the true cycle complexity of cycle using more expensive methods, like block Gauss-Seidel will be underestimated.

Additionally, if the cycle used in practice isn’t a (1,1)-cycle, then this cost estimate will be off.

grid_complexity()[source]

Grid complexity of this multigrid hierarchy

Defined as:
Number of unknowns on all levels / Number of unknowns on the finest level
class level[source]

Stores one level of the multigrid hierarchy

All level objects will have an ‘A’ attribute referencing the matrix of that level. All levels, except for the coarsest level, will also have ‘P’ and ‘R’ attributes referencing the prolongation and restriction operators that act between each level and the next coarser level.

A : csr_matrix
Problem matrix for Ax=b
R : csr_matrix
Restriction matrix between levels (often R = P.T)
P : csr_matrix
Prolongation or Interpolation matrix.

The functionality of this class is a struct

multilevel_solver.operator_complexity()[source]

Operator complexity of this multigrid hierarchy

Defined as:
Number of nonzeros in the matrix on all levels / Number of nonzeros in the matrix on the finest level
multilevel_solver.psolve(b)[source]
multilevel_solver.solve(b, x0=None, tol=1e-05, maxiter=100, cycle='V', accel=None, callback=None, residuals=None, return_residuals=False)[source]

Main solution call to execute multigrid cycling.

b : array
Right hand side.
x0 : array
Initial guess.
tol : float
Stopping criteria: relative residual r[k]/r[0] tolerance.
maxiter : int
Stopping criteria: maximum number of allowable iterations.
cycle : {‘V’,’W’,’F’,’AMLI’}
Type of multigrid cycle to perform in each iteration.
accel : {string, function}
Defines acceleration method. Can be a string such as ‘cg’ or ‘gmres’ which is the name of an iterative solver in pyamg.krylov (preferred) or scipy.sparse.linalg.isolve. If accel is not a string, it will be treated like a function with the same interface provided by the iterative solvers in SciPy.
callback : function
User-defined function called after each iteration. It is called as callback(xk) where xk is the k-th iterate vector.
residuals : list
List to contain residual norms at each iteration.
x : array
Approximate solution to Ax=b

aspreconditioner

>>> from numpy import ones
>>> from pyamg import ruge_stuben_solver
>>> from pyamg.gallery import poisson
>>> A = poisson((100, 100), format='csr')
>>> b = A * ones(A.shape[0])
>>> ml = ruge_stuben_solver(A, max_coarse=10)
>>> residuals = []
>>> x = ml.solve(b, tol=1e-12, residuals=residuals) # standalone solver
pyamg.multilevel.coarse_grid_solver(solver)[source]

Return a coarse grid solver suitable for multilevel_solver

solver : {string, callable, tuple}

The solver method is either (1) a string such as ‘splu’ or ‘pinv’ of a callable object which receives only parameters (A, b) and returns an (approximate or exact) solution to the linear system Ax = b, or (2) a callable object that takes parameters (A,b) and returns an (approximate or exact) solution to Ax = b, or (3) a tuple of the form (string|callable, args), where args is a dictionary of arguments to be passed to the function denoted by string or callable.

The set of valid string arguments is:
  • Sparse direct methods:
    • splu : sparse LU solver
  • Sparse iterative methods:
    • the name of any method in scipy.sparse.linalg.isolve or pyamg.krylov (e.g. ‘cg’). Methods in pyamg.krylov take precedence.
    • relaxation method, such as ‘gauss_seidel’ or ‘jacobi’, present in pyamg.relaxation
  • Dense methods:
    • pinv : pseudoinverse (QR)
    • pinv2 : pseudoinverse (SVD)
    • lu : LU factorization
    • cholesky : Cholesky factorization
ptr : generic_solver
A class for use as a standalone or coarse grids solver
>>> from numpy import ones
>>> from scipy.sparse import spdiags
>>> from pyamg.gallery import poisson
>>> from pyamg import coarse_grid_solver
>>> A = poisson((10, 10), format='csr')
>>> b = A * ones(A.shape[0])
>>> cgs = coarse_grid_solver('lu')
>>> x = cgs(A, b)

The setup Module

pyamg.setup.configuration(parent_package='', top_path=None)[source]

The strength Module

Strength of Connection functions

Requirements for the strength matrix C are:
  1. Nonzero diagonal whenever A has a nonzero diagonal
  2. Non-negative entries (float or bool) in [0,1]
  3. Large entries denoting stronger connections
  4. C denotes nodal connections, i.e., if A is an nxn BSR matrix with row block size of m, then C is (n/m) x (n/m)
pyamg.strength.classical_strength_of_connection(A, theta=0.0)[source]

Return a strength of connection matrix using the classical AMG measure An off-diagonal entry A[i,j] is a strong connection iff:

| A[i,j] | >= theta * max(| A[i,k] |), where k != i
A : csr_matrix or bsr_matrix
Square, sparse matrix in CSR or BSR format
theta : float
Threshold parameter in [0,1].
S : csr_matrix
Matrix graph defining strong connections. S[i,j]=1 if vertex i is strongly influenced by vertex j.

symmetric_strength_of_connection : symmetric measure used in SA evolution_strength_of_connection : relaxation based strength measure

  • A symmetric A does not necessarily yield a symmetric strength matrix S
  • Calls C++ function classical_strength_of_connection
  • The version as implemented is designed form M-matrices. Trottenberg et al. use max A[i,k] over all negative entries, which is the same. A positive edge weight never indicates a strong connection.
[1]Briggs, W. L., Henson, V. E., McCormick, S. F., “A multigrid tutorial”, Second edition. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000. xii+193 pp. ISBN: 0-89871-462-1
[2]Trottenberg, U., Oosterlee, C. W., Schuller, A., “Multigrid”, Academic Press, Inc., San Diego, CA, 2001. xvi+631 pp. ISBN: 0-12-701070-X
>>> import numpy
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.strength import classical_strength_of_connection
>>> n=3
>>> stencil = numpy.array([[-1.0,-1.0,-1.0],
...                        [-1.0, 8.0,-1.0],
...                        [-1.0,-1.0,-1.0]])
>>> A = stencil_grid(stencil, (n,n), format='csr')
>>> S = classical_strength_of_connection(A, 0.0)
pyamg.strength.symmetric_strength_of_connection(A, theta=0)[source]

Compute strength of connection matrix using the standard symmetric measure

An off-diagonal connection A[i,j] is strong iff:

abs(A[i,j]) >= theta * sqrt( abs(A[i,i]) * abs(A[j,j]) )
A : csr_matrix
Matrix graph defined in sparse format. Entry A[i,j] describes the strength of edge [i,j]
theta : float
Threshold parameter (positive).
S : csr_matrix
Matrix graph defining strong connections. S[i,j]=1 if vertex i is strongly influenced by vertex j.

symmetric_strength_of_connection : symmetric measure used in SA evolution_strength_of_connection : relaxation based strength measure

  • For vector problems, standard strength measures may produce undesirable aggregates. A “block approach” from Vanek et al. is used to replace vertex comparisons with block-type comparisons. A connection between nodes i and j in the block case is strong if:

    ||AB[i,j]|| >= theta * sqrt( ||AB[i,i]||*||AB[j,j]|| ) where AB[k,l]
    

    is the matrix block (degrees of freedom) associated with nodes k and l and ||.|| is a matrix norm, such a Frobenius.

[1]Vanek, P. and Mandel, J. and Brezina, M., “Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems”, Computing, vol. 56, no. 3, pp. 179–196, 1996. http://citeseer.ist.psu.edu/vanek96algebraic.html
>>> import numpy
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.strength import symmetric_strength_of_connection
>>> n=3
>>> stencil = numpy.array([[-1.0,-1.0,-1.0],
...                        [-1.0, 8.0,-1.0],
...                        [-1.0,-1.0,-1.0]])
>>> A = stencil_grid(stencil, (n,n), format='csr')
>>> S = symmetric_strength_of_connection(A, 0.0)
pyamg.strength.evolution_strength_of_connection(A, B='ones', epsilon=4.0, k=2, proj_type='l2', block_flag=False, symmetrize_measure=True)[source]

Construct strength of connection matrix using an Evolution-based measure

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix
B : {string, array}
If B=’ones’, then the near nullspace vector used is all ones. If B is an (NxK) array, then B is taken to be the near nullspace vectors.
epsilon : scalar
Drop tolerance
k : integer
ODE num time steps, step size is assumed to be 1/rho(DinvA)
proj_type : {‘l2’,’D_A’}
Define norm for constrained min prob, i.e. define projection
block_flag : {boolean}
If True, use a block D inverse as preconditioner for A during weighted-Jacobi
Atilde : {csr_matrix}
Sparse matrix of strength values
[1]Olson, L. N., Schroder, J., Tuminaro, R. S., “A New Perspective on Strength Measures in Algebraic Multigrid”, submitted, June, 2008.
>>> import numpy
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.strength import evolution_strength_of_connection
>>> n=3
>>> stencil = numpy.array([[-1.0,-1.0,-1.0],
...                        [-1.0, 8.0,-1.0],
...                        [-1.0,-1.0,-1.0]])
>>> A = stencil_grid(stencil, (n,n), format='csr')
>>> S = evolution_strength_of_connection(A, numpy.ones((A.shape[0],1)))
pyamg.strength.distance_strength_of_connection(A, V, theta=2.0, relative_drop=True)[source]

Distance based strength-of-connection

A : csr_matrix or bsr_matrix
Square, sparse matrix in CSR or BSR format
V : array
Coordinates of the vertices of the graph of A
relative_drop : bool
If false, then a connection must be within a distance of theta from a point to be strongly connected. If true, then the closest connection is always strong, and other points must be within theta times the smallest distance to be strong
C : csr_matrix
C(i,j) = distance(point_i, point_j) Strength of connection matrix where strength values are distances, i.e. the smaller the value, the stronger the connection. Sparsity pattern of C is copied from A.
  • theta is a drop tolerance that is applied row-wise
  • If a BSR matrix given, then the return matrix is still CSR. The strength is given between super nodes based on the BSR block size.
>>> from pyamg.gallery import load_example
>>> from pyamg.strength import distance_strength_of_connection
>>> data = load_example('airfoil')
>>> A = data['A'].tocsr()
>>> S = distance_strength_of_connection(data['A'], data['vertices'])
pyamg.strength.algebraic_distance(A, alpha=0.5, R=5, k=20, theta=0.1, p=2)[source]

Construct an AMG strength of connection matrix using an algebraic distance measure.

A : {csr_matrix}
Sparse NxN matrix
alpha : scalar
Weight for Jacobi
R : integer
Number of random vectors
k : integer
Number of relaxation passes
theta : scalar
Drop values larger than theta
p : scalar or inf
p-norm of the measure
C : {csr_matrix}
Sparse matrix of strength values
[1]“Advanced Coarsening Schemes for Graph Partitioning” by Ilya Safro, Peter Sanders, and Christian Schulz

No unit testing yet.

Does not handle BSR matrices yet.

pyamg.strength.ode_strength_of_connection(*args, **kwds)[source]

ode_strength_of_connection is deprecated!

Use evolution_strength_of_connection instead