aggregation Package

aggregation Package

Aggregation-based AMG

pyamg.aggregation.adaptive_sa_solver(A, initial_candidates=None, symmetry='hermitian', pdef=True, num_candidates=1, candidate_iters=5, improvement_iters=0, epsilon=0.1, max_levels=10, max_coarse=100, aggregate='standard', prepostsmoother=('gauss_seidel', {'sweep': 'symmetric'}), smooth=('jacobi', {}), strength='symmetric', coarse_solver='pinv2', eliminate_local=(False, {'Ca': 1.0}), keep=False, **kwargs)

Create a multilevel solver using Adaptive Smoothed Aggregation (aSA)

A : {csr_matrix, bsr_matrix}
Square matrix in CSR or BSR format
initial_candidates : {None, n x m dense matrix}
If a matrix, then this forms the basis for the first m candidates. Also in this case, the initial setup stage is skipped, because this provides the first candidate(s). If None, then a random initial guess and relaxation are used to inform the initial candidate.
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian Note that for the strictly real case, these two options are the same Note that this flag does not denote definiteness of the operator
pdef : {bool}
True or False, whether A is known to be positive definite.
num_candidates : {integer} : default 1
Number of near-nullspace candidates to generate
candidate_iters : {integer} : default 5
Number of smoothing passes/multigrid cycles used at each level of the adaptive setup phase
improvement_iters : {integer} : default 0
Number of times each candidate is improved
epsilon : {float} : default 0.1
Target convergence factor
max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
prepostsmoother : {string or dict}
Pre- and post-smoother used in the adaptive method
strength : [‘symmetric’, ‘classical’, ‘evolution’,
(‘predefined’, {‘C’: csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. See smoothed_aggregation_solver(...) documentation.

aggregate : [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’: csr_matrix})]

Method used to aggregate nodes. See smoothed_aggregation_solver(...) documentation.

smooth : [‘jacobi’, ‘richardson’, ‘energy’, None]
Method used used to smooth the tentative prolongator. See smoothed_aggregation_solver(...) documentation
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

eliminate_local : {tuple}
Length 2 tuple. If the first entry is True, then eliminate candidates where they aren’t needed locally, using the second entry of the tuple to contain arguments to local elimination routine. Given the rigid sparse data structures, this doesn’t help much, if at all, with complexity. Its more of a diagnostic utility.
keep: {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), and aggregation (AggOp) are kept.
multilevel_solver : multilevel_solver
Smoothed aggregation solver with adaptively generated candidates
  • Floating point value representing the “work” required to generate the solver. This value is the total cost of just relaxation, relative to the fine grid. The relaxation method used is assumed to symmetric Gauss-Seidel.
  • Unlike the standard Smoothed Aggregation (SA) method, adaptive SA does not require knowledge of near-nullspace candidate vectors. Instead, an adaptive procedure computes one or more candidates ‘from scratch’. This approach is useful when no candidates are known or the candidates have been invalidated due to changes to matrix A.
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.aggregation import adaptive_sa_solver
>>> import numpy
>>> A=stencil_grid([[-1,-1,-1],[-1,8.0,-1],[-1,-1,-1]],                       (31,31),format='csr')
>>> [asa,work] = adaptive_sa_solver(A,num_candidates=1)
>>> residuals=[]
>>> x=asa.solve(b=numpy.ones((A.shape[0],)), x0=numpy.ones((A.shape[0],)),                    residuals=residuals)
[1]Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge “Adaptive Smoothed Aggregation ($lpha$SA) Multigrid” SIAM Review Volume 47, Issue 2 (2005) http://www.cs.umn.edu/~maclach/research/aSA2.pdf
pyamg.aggregation.energy_prolongation_smoother(A, T, Atilde, B, Bf, Cpt_params, krylov='cg', maxiter=4, tol=1e-08, degree=1, weighting='local')

Minimize the energy of the coarse basis functions (columns of T). Both root-node and non-root-node style prolongation smoothing is available, see Cpt_params description below.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix
T : {bsr_matrix}
Tentative prolongator, a NxM sparse matrix (M < N)
Atilde : {csr_matrix}
Strength of connection matrix
B : {array}
Near-nullspace modes for coarse grid. Has shape (M,k) where k is the number of coarse candidate vectors.
Bf : {array}
Near-nullspace modes for fine grid. Has shape (N,k) where k is the number of coarse candidate vectors.
Cpt_params : {tuple}
Tuple of the form (bool, dict). If the Cpt_params[0] = False, then the standard SA prolongation smoothing is carried out. If True, then root-node style prolongation smoothing is carried out. The dict must be a dictionary of parameters containing, (1) P_I: P_I.T is the injection matrix for the Cpts, (2) I_F: an identity matrix for only the F-points (i.e. I, but with zero rows and columns for C-points) and I_C: the C-point analogue to I_F. See Notes below for more information.
krylov : {string}

‘cg’ : for SPD systems. Solve A T = 0 in a constraint space with CG ‘cgnr’ : for nonsymmetric and/or indefinite systems.

Solve A T = 0 in a constraint space with CGNR
‘gmres’ : for nonsymmetric and/or indefinite systems.
Solve A T = 0 in a constraint space with GMRES
maxiter : integer
Number of energy minimization steps to apply to the prolongator
tol : {scalar}
Minimization tolerance
degree : {int}
Generate sparsity pattern for P based on (Atilde^degree T)
weighting : {string}
‘block’, ‘diagonal’ or ‘local’ construction of the
diagonal preconditioning
‘local’: Uses a local row-wise weight based on the Gershgorin estimate.
Avoids any potential under-damping due to inaccurate spectral radius estimates.

‘block’: If A is a BSR matrix, use a block diagonal inverse of A ‘diagonal’: Use inverse of the diagonal of A

T : {bsr_matrix}
Smoothed prolongator

Only ‘diagonal’ weighting is supported for the CGNR method, because we are working with A^* A and not A.

When Cpt_params[0] == True, root-node style prolongation smoothing is used to minimize the energy of columns of T. Essentially, an identity block is maintained in T, corresponding to injection from the coarse-grid to the fine-grid root-nodes. See [2] for more details, and see util.utils.get_Cpt_params for the helper function to generate Cpt_params.

If Cpt_params[0] == False, the energy of columns of T are still minimized, but without maintaining the identity block.

>>> from pyamg.aggregation import energy_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> print T.todense()
[[ 1.  0.]
 [ 1.  0.]
 [ 1.  0.]
 [ 0.  1.]
 [ 0.  1.]
 [ 0.  1.]]
>>> A = poisson((6,),format='csr')
>>> B = numpy.ones((2,1),dtype=float)
>>> P = energy_prolongation_smoother(A,T,A,B, None, (False,{}))
>>> print P.todense()
[[ 1.          0.        ]
 [ 1.          0.        ]
 [ 0.66666667  0.33333333]
 [ 0.33333333  0.66666667]
 [ 0.          1.        ]
 [ 0.          1.        ]]
[1]Jan Mandel, Marian Brezina, and Petr Vanek “Energy Optimization of Algebraic Multigrid Bases” Computing 62, 205-228, 1999 http://dx.doi.org/10.1007/s006070050022
[2]Olson, L. and Schroder, J. and Tuminaro, R., “A general interpolation strategy for algebraic multigrid using energy minimization”, SIAM Journal on Scientific Computing (SISC), vol. 33, pp. 966–991, 2011.
pyamg.aggregation.fit_candidates(AggOp, B, tol=1e-10)

Fit near-nullspace candidates to form the tentative prolongator

AggOp : csr_matrix
Describes the sparsity pattern of the tentative prolongator. Has dimension (#blocks, #aggregates)
B : array
The near-nullspace candidates stored in column-wise fashion. Has dimension (#blocks * blocksize, #candidates)
tol : scalar
Threshold for eliminating local basis functions. If after orthogonalization a local basis function Q[:, j] is small, i.e. ||Q[:, j]|| < tol, then Q[:, j] is set to zero.
(Q, R) : (bsr_matrix, array)
The tentative prolongator Q is a sparse block matrix with dimensions (#blocks * blocksize, #aggregates * #candidates) formed by dense blocks of size (blocksize, #candidates). The coarse level candidates are stored in R which has dimensions (#aggregates * #candidates, #candidates).

amg_core.fit_candidates

Assuming that each row of AggOp contains exactly one non-zero entry, i.e. all unknowns belong to an aggregate, then Q and R satisfy the relationship B = Q*R. In other words, the near-nullspace candidates are represented exactly by the tentative prolongator.

If AggOp contains rows with no non-zero entries, then the range of the tentative prolongator will not include those degrees of freedom. This situation is illustrated in the examples below.

[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
>>> from scipy.sparse import csr_matrix
>>> from pyamg.aggregation.tentative import fit_candidates
>>> # four nodes divided into two aggregates
... AggOp = csr_matrix( [[1, 0],
...                      [1, 0],
...                      [0, 1],
...                      [0, 1]] )
>>> # B contains one candidate, the constant vector
... B = [[1],
...      [1],
...      [1],
...      [1]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678,  0.        ],
        [ 0.70710678,  0.        ],
        [ 0.        ,  0.70710678],
        [ 0.        ,  0.70710678]])
>>> R
array([[ 1.41421356],
       [ 1.41421356]])
>>> # Two candidates, the constant vector and a linear function
... B = [[1, 0],
...      [1, 1],
...      [1, 2],
...      [1, 3]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678, -0.70710678,  0.        ,  0.        ],
        [ 0.70710678,  0.70710678,  0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.70710678, -0.70710678],
        [ 0.        ,  0.        ,  0.70710678,  0.70710678]])
>>> R
array([[ 1.41421356,  0.70710678],
       [ 0.        ,  0.70710678],
       [ 1.41421356,  3.53553391],
       [ 0.        ,  0.70710678]])
>>> # aggregation excludes the third node
... AggOp = csr_matrix( [[1, 0],
...                      [1, 0],
...                      [0, 0],
...                      [0, 1]] )
>>> B = [[1],
...      [1],
...      [1],
...      [1]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678,  0.        ],
        [ 0.70710678,  0.        ],
        [ 0.        ,  0.        ],
        [ 0.        ,  1.        ]])
>>> R
array([[ 1.41421356],
       [ 1.        ]])
pyamg.aggregation.jacobi_prolongation_smoother(S, T, C, B, omega=1.3333333333333333, degree=1, filter=False, weighting='diagonal')

Jacobi prolongation smoother

S : {csr_matrix, bsr_matrix}
Sparse NxN matrix used for smoothing. Typically, A.
T : {csr_matrix, bsr_matrix}
Tentative prolongator
C : {csr_matrix, bsr_matrix}
Strength-of-connection matrix
B : {array}
Near nullspace modes for the coarse grid such that T*B exactly reproduces the fine grid near nullspace modes
omega : {scalar}
Damping parameter
filter : {boolean}
If true, filter S before smoothing T. This option can greatly control complexity.
weighting : {string}

‘block’, ‘diagonal’ or ‘local’ weighting for constructing the Jacobi D ‘local’: Uses a local row-wise weight based on the Gershgorin estimate.

Avoids any potential under-damping due to inaccurate spectral radius estimates.

‘block’: If A is a BSR matrix, use a block diagonal inverse of A ‘diagonal’: Classic Jacobi D = diagonal(A)

P : {csr_matrix, bsr_matrix}
Smoothed (final) prolongator defined by P = (I - omega/rho(K) K) * T where K = diag(S)^-1 * S and rho(K) is an approximation to the spectral radius of K.

If weighting is not ‘local’, then results using Jacobi prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.

>>> from pyamg.aggregation import jacobi_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1.,  0.],
        [ 1.,  0.],
        [ 1.,  0.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.]])
>>> A = poisson((6,),format='csr')
>>> P = jacobi_prolongation_smoother(A,T,A,numpy.ones((2,1)))
>>> P.todense()
matrix([[ 0.64930164,  0.        ],
        [ 1.        ,  0.        ],
        [ 0.64930164,  0.35069836],
        [ 0.35069836,  0.64930164],
        [ 0.        ,  1.        ],
        [ 0.        ,  0.64930164]])
pyamg.aggregation.lloyd_aggregation(C, ratio=0.03, distance='unit', maxiter=10)

Aggregated nodes using Lloyd Clustering

C : csr_matrix
strength of connection matrix
ratio : scalar
Fraction of the nodes which will be seeds.
distance : [‘unit’,’abs’,’inv’,None]

Distance assigned to each edge of the graph G used in Lloyd clustering

For each nonzero value C[i,j]:

‘unit’ G[i,j] = 1
‘abs’ G[i,j] = abs(C[i,j])
‘inv’ G[i,j] = 1.0/abs(C[i,j])
‘same’ G[i,j] = C[i,j]
‘sub’ G[i,j] = C[i,j] - min(C)
maxiter : int
Maximum number of iterations to perform
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
seeds : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i

amg_core.standard_aggregation

>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import lloyd_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> lloyd_aggregation(A)[0].todense() # one aggregate
matrix([[1],
        [1],
        [1],
        [1]], dtype=int8)
>>> # more seeding for two aggregates
>>> Agg = lloyd_aggregation(A,ratio=0.5)[0].todense()
pyamg.aggregation.naive_aggregation(C)

Compute the sparsity pattern of the tentative prolongator

C : csr_matrix
strength of connection matrix
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
Cpts : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import naive_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [1, 0],
        [0, 1],
        [0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense()                      # first vertex is isolated
matrix([[1, 0, 0],
        [0, 1, 1],
        [0, 1, 1]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [0, 1],
        [0, 1]], dtype=int8)

amg_core.naive_aggregation

Differs from standard aggregation. Each dof is considered. If it has been aggregated, skip over. Otherwise, put dof and any unaggregated neighbors in an aggregate. Results in possibly much higher complexities than standard aggregation.

pyamg.aggregation.richardson_prolongation_smoother(S, T, omega=1.3333333333333333, degree=1)

Richardson prolongation smoother

S : {csr_matrix, bsr_matrix}
Sparse NxN matrix used for smoothing. Typically, A or the “filtered matrix” obtained from A by lumping weak connections onto the diagonal of A.
T : {csr_matrix, bsr_matrix}
Tentative prolongator
omega : {scalar}
Damping parameter
P : {csr_matrix, bsr_matrix}
Smoothed (final) prolongator defined by P = (I - omega/rho(S) S) * T where rho(S) is an approximation to the spectral radius of S.

Results using Richardson prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.

>>> from pyamg.aggregation import richardson_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1.,  0.],
        [ 1.,  0.],
        [ 1.,  0.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.]])
>>> A = poisson((6,),format='csr')
>>> P = richardson_prolongation_smoother(A,T)
>>> P.todense()
matrix([[ 0.64930164,  0.        ],
        [ 1.        ,  0.        ],
        [ 0.64930164,  0.35069836],
        [ 0.35069836,  0.64930164],
        [ 0.        ,  1.        ],
        [ 0.        ,  0.64930164]])
pyamg.aggregation.rootnode_solver(A, B=None, BH=None, symmetry='hermitian', strength='symmetric', aggregate='standard', smooth='energy', 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=500, diagonal_dominance=False, keep=False, **kwargs)

Create a multilevel solver using root-node based Smoothed Aggregation (SA). See the notes below, for the major differences with the classical-style smoothed aggregation solver in aggregation.smoothed_aggregation_solver.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix in CSR or BSR format
B : {None, array_like}
Right near-nullspace candidates stored in the columns of an NxK array. K must be >= the blocksize of A (see reference [2]). The default value B=None is equivalent to choosing the constant over each block-variable, B=np.kron(np.ones((A.shape[0]/blocksize(A), 1)), np.eye(blocksize(A)))
BH : {None, array_like}
Left near-nullspace candidates stored in the columns of an NxK array. BH is only used if symmetry is ‘nonsymmetric’. K must be >= the blocksize of A (see reference [2]). The default value B=None is equivalent to choosing the constant over each block-variable, B=np.kron(np.ones((A.shape[0]/blocksize(A), 1)), np.eye(blocksize(A)))
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian ‘nonsymmetric’ i.e. nonsymmetric in a hermitian sense Note that for the strictly real case, symmetric and hermitian are the same Note that this flag does not denote definiteness of the operator.
strength : {list} : default
[‘symmetric’, ‘classical’, ‘evolution’, (‘predefined’,
{‘C’ : csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. Method-specific parameters may be passed in using a tuple, e.g. strength=(‘symmetric’,{‘theta’ : 0.25 }). If strength=None, all nonzero entries of the matrix are considered strong. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined strength matrix on each level.

aggregate : {list} : default [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’ : csr_matrix})]

Method used to aggregate nodes. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined aggregation on each level.

smooth : {list} : default [‘energy’, None]
Method used to smooth the tentative prolongator. Method-specific parameters may be passed in using a tuple, e.g. smooth= (‘energy’,{‘krylov’ : ‘gmres’}). Only ‘energy’ and None are valid prolongation smoothing options. See notes below for varying this parameter on a per level basis.
presmoother : {tuple, string, list} : default (‘block_gauss_seidel’,
{‘sweep’:’symmetric’})

Defines the presmoother for the multilevel cycling. The default block Gauss-Seidel option defaults to point-wise Gauss-Seidel, if the matrix is CSR or is a BSR matrix with blocksize of 1. See notes below for varying this parameter on a per level basis.

postsmoother : {tuple, string, list}
Same as presmoother, except defines the postsmoother.
improve_candidates : {tuple, string, list} : default
[(‘block_gauss_seidel’,
{‘sweep’: ‘symmetric’, ‘iterations’: 4}), None]

The ith entry defines the method used to improve the candidates B on level i. If the list is shorter than max_levels, then the last entry will define the method for all levels lower. If tuple or string, then this single relaxation descriptor defines improve_candidates on all levels. The list elements are relaxation descriptors of the form used for presmoother and postsmoother. A value of None implies no action on B.

max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
diagonal_dominance : {bool, tuple} : default False
If True (or the first tuple entry is True), then avoid coarsening diagonally dominant rows. The second tuple entry requires a dictionary, where the key value ‘theta’ is used to tune the diagonal dominance threshold.
keep : {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), aggregation (AggOp), and arrays storing the C-points (Cpts) and F-points (Fpts) are kept at each level.
cycle_type : [‘V’,’W’,’F’]
Structrure of multigrid cycle
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

ml : multilevel_solver
Multigrid hierarchy of matrices and prolongation operators

multilevel_solver, aggregation.smoothed_aggregation_solver, classical.ruge_stuben_solver

  • Root-node style SA differs from classical SA primarily by preserving and identity block in the interpolation operator, P. Each aggregate has a “root-node” or “center-node” associated with it, and this root-node is injected from the coarse grid to the fine grid. The injection corresponds to the identity block.
  • Only smooth={‘energy’, None} is supported for prolongation smoothing. See reference [2] below for more details on why the ‘energy’ prolongation smoother is the natural counterpart to root-node style SA.
  • The additional parameters are passed through as arguments to multilevel_solver. Refer to pyamg.multilevel_solver for additional documentation.

  • At each level, four steps are executed in order to define the coarser level operator.

    1. Matrix A is given and used to derive a strength matrix, C.
    2. Based on the strength matrix, indices are grouped or aggregated.
    3. The aggregates define coarse nodes and a tentative prolongation operator T is defined by injection
    4. The tentative prolongation operator is smoothed by a relaxation scheme to improve the quality and extent of interpolation from the aggregates to fine nodes.
  • The parameters smooth, strength, aggregate, presmoother, postsmoother can be varied on a per level basis. For different methods on different levels, use a list as input so that the i-th entry defines the method at the i-th level. If there are more levels in the hierarchy than list entries, the last entry will define the method for all levels lower.

    Examples are: smooth=[(‘jacobi’, {‘omega’:1.0}), None, ‘jacobi’] presmoother=[(‘block_gauss_seidel’, {‘sweep’:symmetric}), ‘sor’] aggregate=[‘standard’, ‘naive’] strength=[(‘symmetric’, {‘theta’:0.25}),

    (‘symmetric’, {‘theta’:0.08})]

  • Predefined strength of connection and aggregation schemes can be specified. These options are best used together, but aggregation can be predefined while strength of connection is not.

    For predefined strength of connection, use a list consisting of tuples of the form (‘predefined’, {‘C’ : C0}), where C0 is a csr_matrix and each degree-of-freedom in C0 represents a supernode. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘C’ : C0}), (‘predefined’, {‘C’ : C1}) ].

    Similarly for predefined aggregation, use a list of tuples. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘AggOp’ : Agg0}), (‘predefined’, {‘AggOp’ : Agg1}) ], where the dimensions of A, Agg0 and Agg1 are compatible, i.e. Agg0.shape[1] == A.shape[0] and Agg1.shape[1] == Agg0.shape[0]. Each AggOp is a csr_matrix.

    Because this is a root-nodes solver, if a member of the predefined aggregation list is predefined, it must be of the form (‘predefined’, {‘AggOp’ : Agg, ‘Cnodes’ : Cnodes}).

>>> from pyamg import rootnode_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> import numpy
>>> A = poisson((100, 100), format='csr')           # matrix
>>> b = numpy.ones((A.shape[0]))                   # RHS
>>> ml = rootnode_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
[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
[2]Olson, L. and Schroder, J. and Tuminaro, R., “A general interpolation strategy for algebraic multigrid using energy minimization”, SIAM Journal on Scientific Computing (SISC), vol. 33, pp. 966–991, 2011.
pyamg.aggregation.smoothed_aggregation_solver(A, B=None, BH=None, symmetry='hermitian', strength='symmetric', aggregate='standard', smooth=('jacobi', {'omega': 1.3333333333333333}), 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=500, diagonal_dominance=False, keep=False, **kwargs)

Create a multilevel solver using classical-style Smoothed Aggregation (SA)

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix in CSR or BSR format
B : {None, array_like}
Right near-nullspace candidates stored in the columns of an NxK array. The default value B=None is equivalent to B=ones((N,1))
BH : {None, array_like}
Left near-nullspace candidates stored in the columns of an NxK array. BH is only used if symmetry is ‘nonsymmetric’. The default value B=None is equivalent to BH=B.copy()
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian ‘nonsymmetric’ i.e. nonsymmetric in a hermitian sense Note, in the strictly real case, symmetric and hermitian are the same Note, this flag does not denote definiteness of the operator.
strength : {list} : default [‘symmetric’, ‘classical’, ‘evolution’,
(‘predefined’, {‘C’ : csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. Method-specific parameters may be passed in using a tuple, e.g. strength=(‘symmetric’,{‘theta’ : 0.25 }). If strength=None, all nonzero entries of the matrix are considered strong. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined strength matrix on each level.

aggregate : {list} : default [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’ : csr_matrix})]

Method used to aggregate nodes. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined aggregation on each level.

smooth : {list} : default [‘jacobi’, ‘richardson’, ‘energy’, None]
Method used to smooth the tentative prolongator. Method-specific parameters may be passed in using a tuple, e.g. smooth= (‘jacobi’,{‘filter’ : True }). See notes below for varying this parameter on a per level basis.
presmoother : {tuple, string, list} : default (‘block_gauss_seidel’,
{‘sweep’:’symmetric’})

Defines the presmoother for the multilevel cycling. The default block Gauss-Seidel option defaults to point-wise Gauss-Seidel, if the matrix is CSR or is a BSR matrix with blocksize of 1. See notes below for varying this parameter on a per level basis.

postsmoother : {tuple, string, list}
Same as presmoother, except defines the postsmoother.
improve_candidates : {tuple, string, list} : default
[(‘block_gauss_seidel’,
{‘sweep’: ‘symmetric’, ‘iterations’: 4}), None]

The ith entry defines the method used to improve the candidates B on level i. If the list is shorter than max_levels, then the last entry will define the method for all levels lower. If tuple or string, then this single relaxation descriptor defines improve_candidates on all levels. The list elements are relaxation descriptors of the form used for presmoother and postsmoother. A value of None implies no action on B.

max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
diagonal_dominance : {bool, tuple} : default False
If True (or the first tuple entry is True), then avoid coarsening diagonally dominant rows. The second tuple entry requires a dictionary, where the key value ‘theta’ is used to tune the diagonal dominance threshold.
keep : {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), and aggregation (AggOp) are kept.
cycle_type : [‘V’,’W’,’F’]
Structrure of multigrid cycle
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

ml : multilevel_solver
Multigrid hierarchy of matrices and prolongation operators

multilevel_solver, classical.ruge_stuben_solver, aggregation.smoothed_aggregation_solver

  • This method implements classical-style SA, not root-node style SA (see aggregation.rootnode_solver).

  • The additional parameters are passed through as arguments to multilevel_solver. Refer to pyamg.multilevel_solver for additional documentation.

  • At each level, four steps are executed in order to define the coarser level operator.

    1. Matrix A is given and used to derive a strength matrix, C.
    2. Based on the strength matrix, indices are grouped or aggregated.
    3. The aggregates define coarse nodes and a tentative prolongation operator T is defined by injection
    4. The tentative prolongation operator is smoothed by a relaxation scheme to improve the quality and extent of interpolation from the aggregates to fine nodes.
  • The parameters smooth, strength, aggregate, presmoother, postsmoother can be varied on a per level basis. For different methods on different levels, use a list as input so that the i-th entry defines the method at the i-th level. If there are more levels in the hierarchy than list entries, the last entry will define the method for all levels lower.

    Examples are: smooth=[(‘jacobi’, {‘omega’:1.0}), None, ‘jacobi’] presmoother=[(‘block_gauss_seidel’, {‘sweep’:symmetric}), ‘sor’] aggregate=[‘standard’, ‘naive’] strength=[(‘symmetric’, {‘theta’:0.25}), (‘symmetric’,

    {‘theta’:0.08})]

  • Predefined strength of connection and aggregation schemes can be specified. These options are best used together, but aggregation can be predefined while strength of connection is not.

    For predefined strength of connection, use a list consisting of tuples of the form (‘predefined’, {‘C’ : C0}), where C0 is a csr_matrix and each degree-of-freedom in C0 represents a supernode. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘C’ : C0}), (‘predefined’, {‘C’ : C1}) ].

    Similarly for predefined aggregation, use a list of tuples. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘AggOp’ : Agg0}), (‘predefined’, {‘AggOp’ : Agg1}) ], where the dimensions of A, Agg0 and Agg1 are compatible, i.e. Agg0.shape[1] == A.shape[0] and Agg1.shape[1] == Agg0.shape[0]. Each AggOp is a csr_matrix.

>>> from pyamg import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> import numpy
>>> A = poisson((100,100), format='csr')           # matrix
>>> b = numpy.ones((A.shape[0]))                   # 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
[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
pyamg.aggregation.standard_aggregation(C)

Compute the sparsity pattern of the tentative prolongator

C : csr_matrix
strength of connection matrix
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
Cpts : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import standard_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> standard_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [1, 0],
        [0, 1],
        [0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense()                      # first vertex is isolated
matrix([[1, 0, 0],
        [0, 1, 1],
        [0, 1, 1]])
>>> standard_aggregation(A)[0].todense() # one aggregate
matrix([[0],
        [1],
        [1]], dtype=int8)

amg_core.standard_aggregation

adaptive Module

Adaptive Smoothed Aggregation

pyamg.aggregation.adaptive.adaptive_sa_solver(A, initial_candidates=None, symmetry='hermitian', pdef=True, num_candidates=1, candidate_iters=5, improvement_iters=0, epsilon=0.1, max_levels=10, max_coarse=100, aggregate='standard', prepostsmoother=('gauss_seidel', {'sweep': 'symmetric'}), smooth=('jacobi', {}), strength='symmetric', coarse_solver='pinv2', eliminate_local=(False, {'Ca': 1.0}), keep=False, **kwargs)[source]

Create a multilevel solver using Adaptive Smoothed Aggregation (aSA)

A : {csr_matrix, bsr_matrix}
Square matrix in CSR or BSR format
initial_candidates : {None, n x m dense matrix}
If a matrix, then this forms the basis for the first m candidates. Also in this case, the initial setup stage is skipped, because this provides the first candidate(s). If None, then a random initial guess and relaxation are used to inform the initial candidate.
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian Note that for the strictly real case, these two options are the same Note that this flag does not denote definiteness of the operator
pdef : {bool}
True or False, whether A is known to be positive definite.
num_candidates : {integer} : default 1
Number of near-nullspace candidates to generate
candidate_iters : {integer} : default 5
Number of smoothing passes/multigrid cycles used at each level of the adaptive setup phase
improvement_iters : {integer} : default 0
Number of times each candidate is improved
epsilon : {float} : default 0.1
Target convergence factor
max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
prepostsmoother : {string or dict}
Pre- and post-smoother used in the adaptive method
strength : [‘symmetric’, ‘classical’, ‘evolution’,
(‘predefined’, {‘C’: csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. See smoothed_aggregation_solver(...) documentation.

aggregate : [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’: csr_matrix})]

Method used to aggregate nodes. See smoothed_aggregation_solver(...) documentation.

smooth : [‘jacobi’, ‘richardson’, ‘energy’, None]
Method used used to smooth the tentative prolongator. See smoothed_aggregation_solver(...) documentation
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

eliminate_local : {tuple}
Length 2 tuple. If the first entry is True, then eliminate candidates where they aren’t needed locally, using the second entry of the tuple to contain arguments to local elimination routine. Given the rigid sparse data structures, this doesn’t help much, if at all, with complexity. Its more of a diagnostic utility.
keep: {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), and aggregation (AggOp) are kept.
multilevel_solver : multilevel_solver
Smoothed aggregation solver with adaptively generated candidates
  • Floating point value representing the “work” required to generate the solver. This value is the total cost of just relaxation, relative to the fine grid. The relaxation method used is assumed to symmetric Gauss-Seidel.
  • Unlike the standard Smoothed Aggregation (SA) method, adaptive SA does not require knowledge of near-nullspace candidate vectors. Instead, an adaptive procedure computes one or more candidates ‘from scratch’. This approach is useful when no candidates are known or the candidates have been invalidated due to changes to matrix A.
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.aggregation import adaptive_sa_solver
>>> import numpy
>>> A=stencil_grid([[-1,-1,-1],[-1,8.0,-1],[-1,-1,-1]],                       (31,31),format='csr')
>>> [asa,work] = adaptive_sa_solver(A,num_candidates=1)
>>> residuals=[]
>>> x=asa.solve(b=numpy.ones((A.shape[0],)), x0=numpy.ones((A.shape[0],)),                    residuals=residuals)
[1]Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge “Adaptive Smoothed Aggregation ($lpha$SA) Multigrid” SIAM Review Volume 47, Issue 2 (2005) http://www.cs.umn.edu/~maclach/research/aSA2.pdf

aggregate Module

Aggregation methods

pyamg.aggregation.aggregate.standard_aggregation(C)[source]

Compute the sparsity pattern of the tentative prolongator

C : csr_matrix
strength of connection matrix
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
Cpts : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import standard_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> standard_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [1, 0],
        [0, 1],
        [0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense()                      # first vertex is isolated
matrix([[1, 0, 0],
        [0, 1, 1],
        [0, 1, 1]])
>>> standard_aggregation(A)[0].todense() # one aggregate
matrix([[0],
        [1],
        [1]], dtype=int8)

amg_core.standard_aggregation

pyamg.aggregation.aggregate.naive_aggregation(C)[source]

Compute the sparsity pattern of the tentative prolongator

C : csr_matrix
strength of connection matrix
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
Cpts : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import naive_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [1, 0],
        [0, 1],
        [0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense()                      # first vertex is isolated
matrix([[1, 0, 0],
        [0, 1, 1],
        [0, 1, 1]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
        [0, 1],
        [0, 1]], dtype=int8)

amg_core.naive_aggregation

Differs from standard aggregation. Each dof is considered. If it has been aggregated, skip over. Otherwise, put dof and any unaggregated neighbors in an aggregate. Results in possibly much higher complexities than standard aggregation.

pyamg.aggregation.aggregate.lloyd_aggregation(C, ratio=0.03, distance='unit', maxiter=10)[source]

Aggregated nodes using Lloyd Clustering

C : csr_matrix
strength of connection matrix
ratio : scalar
Fraction of the nodes which will be seeds.
distance : [‘unit’,’abs’,’inv’,None]

Distance assigned to each edge of the graph G used in Lloyd clustering

For each nonzero value C[i,j]:

‘unit’ G[i,j] = 1
‘abs’ G[i,j] = abs(C[i,j])
‘inv’ G[i,j] = 1.0/abs(C[i,j])
‘same’ G[i,j] = C[i,j]
‘sub’ G[i,j] = C[i,j] - min(C)
maxiter : int
Maximum number of iterations to perform
AggOp : csr_matrix
aggregation operator which determines the sparsity pattern of the tentative prolongator
seeds : array
array of Cpts, i.e., Cpts[i] = root node of aggregate i

amg_core.standard_aggregation

>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import lloyd_aggregation
>>> A = poisson((4,), format='csr')   # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1.,  0.,  0.],
        [-1.,  2., -1.,  0.],
        [ 0., -1.,  2., -1.],
        [ 0.,  0., -1.,  2.]])
>>> lloyd_aggregation(A)[0].todense() # one aggregate
matrix([[1],
        [1],
        [1],
        [1]], dtype=int8)
>>> # more seeding for two aggregates
>>> Agg = lloyd_aggregation(A,ratio=0.5)[0].todense()

aggregation Module

Support for aggregation-based AMG

pyamg.aggregation.aggregation.smoothed_aggregation_solver(A, B=None, BH=None, symmetry='hermitian', strength='symmetric', aggregate='standard', smooth=('jacobi', {'omega': 1.3333333333333333}), 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=500, diagonal_dominance=False, keep=False, **kwargs)[source]

Create a multilevel solver using classical-style Smoothed Aggregation (SA)

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix in CSR or BSR format
B : {None, array_like}
Right near-nullspace candidates stored in the columns of an NxK array. The default value B=None is equivalent to B=ones((N,1))
BH : {None, array_like}
Left near-nullspace candidates stored in the columns of an NxK array. BH is only used if symmetry is ‘nonsymmetric’. The default value B=None is equivalent to BH=B.copy()
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian ‘nonsymmetric’ i.e. nonsymmetric in a hermitian sense Note, in the strictly real case, symmetric and hermitian are the same Note, this flag does not denote definiteness of the operator.
strength : {list} : default [‘symmetric’, ‘classical’, ‘evolution’,
(‘predefined’, {‘C’ : csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. Method-specific parameters may be passed in using a tuple, e.g. strength=(‘symmetric’,{‘theta’ : 0.25 }). If strength=None, all nonzero entries of the matrix are considered strong. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined strength matrix on each level.

aggregate : {list} : default [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’ : csr_matrix})]

Method used to aggregate nodes. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined aggregation on each level.

smooth : {list} : default [‘jacobi’, ‘richardson’, ‘energy’, None]
Method used to smooth the tentative prolongator. Method-specific parameters may be passed in using a tuple, e.g. smooth= (‘jacobi’,{‘filter’ : True }). See notes below for varying this parameter on a per level basis.
presmoother : {tuple, string, list} : default (‘block_gauss_seidel’,
{‘sweep’:’symmetric’})

Defines the presmoother for the multilevel cycling. The default block Gauss-Seidel option defaults to point-wise Gauss-Seidel, if the matrix is CSR or is a BSR matrix with blocksize of 1. See notes below for varying this parameter on a per level basis.

postsmoother : {tuple, string, list}
Same as presmoother, except defines the postsmoother.
improve_candidates : {tuple, string, list} : default
[(‘block_gauss_seidel’,
{‘sweep’: ‘symmetric’, ‘iterations’: 4}), None]

The ith entry defines the method used to improve the candidates B on level i. If the list is shorter than max_levels, then the last entry will define the method for all levels lower. If tuple or string, then this single relaxation descriptor defines improve_candidates on all levels. The list elements are relaxation descriptors of the form used for presmoother and postsmoother. A value of None implies no action on B.

max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
diagonal_dominance : {bool, tuple} : default False
If True (or the first tuple entry is True), then avoid coarsening diagonally dominant rows. The second tuple entry requires a dictionary, where the key value ‘theta’ is used to tune the diagonal dominance threshold.
keep : {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), and aggregation (AggOp) are kept.
cycle_type : [‘V’,’W’,’F’]
Structrure of multigrid cycle
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

ml : multilevel_solver
Multigrid hierarchy of matrices and prolongation operators

multilevel_solver, classical.ruge_stuben_solver, aggregation.smoothed_aggregation_solver

  • This method implements classical-style SA, not root-node style SA (see aggregation.rootnode_solver).

  • The additional parameters are passed through as arguments to multilevel_solver. Refer to pyamg.multilevel_solver for additional documentation.

  • At each level, four steps are executed in order to define the coarser level operator.

    1. Matrix A is given and used to derive a strength matrix, C.
    2. Based on the strength matrix, indices are grouped or aggregated.
    3. The aggregates define coarse nodes and a tentative prolongation operator T is defined by injection
    4. The tentative prolongation operator is smoothed by a relaxation scheme to improve the quality and extent of interpolation from the aggregates to fine nodes.
  • The parameters smooth, strength, aggregate, presmoother, postsmoother can be varied on a per level basis. For different methods on different levels, use a list as input so that the i-th entry defines the method at the i-th level. If there are more levels in the hierarchy than list entries, the last entry will define the method for all levels lower.

    Examples are: smooth=[(‘jacobi’, {‘omega’:1.0}), None, ‘jacobi’] presmoother=[(‘block_gauss_seidel’, {‘sweep’:symmetric}), ‘sor’] aggregate=[‘standard’, ‘naive’] strength=[(‘symmetric’, {‘theta’:0.25}), (‘symmetric’,

    {‘theta’:0.08})]

  • Predefined strength of connection and aggregation schemes can be specified. These options are best used together, but aggregation can be predefined while strength of connection is not.

    For predefined strength of connection, use a list consisting of tuples of the form (‘predefined’, {‘C’ : C0}), where C0 is a csr_matrix and each degree-of-freedom in C0 represents a supernode. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘C’ : C0}), (‘predefined’, {‘C’ : C1}) ].

    Similarly for predefined aggregation, use a list of tuples. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘AggOp’ : Agg0}), (‘predefined’, {‘AggOp’ : Agg1}) ], where the dimensions of A, Agg0 and Agg1 are compatible, i.e. Agg0.shape[1] == A.shape[0] and Agg1.shape[1] == Agg0.shape[0]. Each AggOp is a csr_matrix.

>>> from pyamg import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> import numpy
>>> A = poisson((100,100), format='csr')           # matrix
>>> b = numpy.ones((A.shape[0]))                   # 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
[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

rootnode Module

Support for aggregation-based AMG

pyamg.aggregation.rootnode.rootnode_solver(A, B=None, BH=None, symmetry='hermitian', strength='symmetric', aggregate='standard', smooth='energy', 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=500, diagonal_dominance=False, keep=False, **kwargs)[source]

Create a multilevel solver using root-node based Smoothed Aggregation (SA). See the notes below, for the major differences with the classical-style smoothed aggregation solver in aggregation.smoothed_aggregation_solver.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix in CSR or BSR format
B : {None, array_like}
Right near-nullspace candidates stored in the columns of an NxK array. K must be >= the blocksize of A (see reference [2]). The default value B=None is equivalent to choosing the constant over each block-variable, B=np.kron(np.ones((A.shape[0]/blocksize(A), 1)), np.eye(blocksize(A)))
BH : {None, array_like}
Left near-nullspace candidates stored in the columns of an NxK array. BH is only used if symmetry is ‘nonsymmetric’. K must be >= the blocksize of A (see reference [2]). The default value B=None is equivalent to choosing the constant over each block-variable, B=np.kron(np.ones((A.shape[0]/blocksize(A), 1)), np.eye(blocksize(A)))
symmetry : {string}
‘symmetric’ refers to both real and complex symmetric ‘hermitian’ refers to both complex Hermitian and real Hermitian ‘nonsymmetric’ i.e. nonsymmetric in a hermitian sense Note that for the strictly real case, symmetric and hermitian are the same Note that this flag does not denote definiteness of the operator.
strength : {list} : default
[‘symmetric’, ‘classical’, ‘evolution’, (‘predefined’,
{‘C’ : csr_matrix}), None]

Method used to determine the strength of connection between unknowns of the linear system. Method-specific parameters may be passed in using a tuple, e.g. strength=(‘symmetric’,{‘theta’ : 0.25 }). If strength=None, all nonzero entries of the matrix are considered strong. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined strength matrix on each level.

aggregate : {list} : default [‘standard’, ‘lloyd’, ‘naive’,
(‘predefined’, {‘AggOp’ : csr_matrix})]

Method used to aggregate nodes. See notes below for varying this parameter on a per level basis. Also, see notes below for using a predefined aggregation on each level.

smooth : {list} : default [‘energy’, None]
Method used to smooth the tentative prolongator. Method-specific parameters may be passed in using a tuple, e.g. smooth= (‘energy’,{‘krylov’ : ‘gmres’}). Only ‘energy’ and None are valid prolongation smoothing options. See notes below for varying this parameter on a per level basis.
presmoother : {tuple, string, list} : default (‘block_gauss_seidel’,
{‘sweep’:’symmetric’})

Defines the presmoother for the multilevel cycling. The default block Gauss-Seidel option defaults to point-wise Gauss-Seidel, if the matrix is CSR or is a BSR matrix with blocksize of 1. See notes below for varying this parameter on a per level basis.

postsmoother : {tuple, string, list}
Same as presmoother, except defines the postsmoother.
improve_candidates : {tuple, string, list} : default
[(‘block_gauss_seidel’,
{‘sweep’: ‘symmetric’, ‘iterations’: 4}), None]

The ith entry defines the method used to improve the candidates B on level i. If the list is shorter than max_levels, then the last entry will define the method for all levels lower. If tuple or string, then this single relaxation descriptor defines improve_candidates on all levels. The list elements are relaxation descriptors of the form used for presmoother and postsmoother. A value of None implies no action on B.

max_levels : {integer} : default 10
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 500
Maximum number of variables permitted on the coarse grid.
diagonal_dominance : {bool, tuple} : default False
If True (or the first tuple entry is True), then avoid coarsening diagonally dominant rows. The second tuple entry requires a dictionary, where the key value ‘theta’ is used to tune the diagonal dominance threshold.
keep : {bool} : default False
Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C), tentative prolongation (T), aggregation (AggOp), and arrays storing the C-points (Cpts) and F-points (Fpts) are kept at each level.
cycle_type : [‘V’,’W’,’F’]
Structrure of multigrid cycle
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
Solver used at the coarsest level of the MG hierarchy.
Optionally, may be a tuple (fn, args), where fn is a string such as

[‘splu’, ‘lu’, ...] or a callable function, and args is a dictionary of arguments to be passed to fn.

ml : multilevel_solver
Multigrid hierarchy of matrices and prolongation operators

multilevel_solver, aggregation.smoothed_aggregation_solver, classical.ruge_stuben_solver

  • Root-node style SA differs from classical SA primarily by preserving and identity block in the interpolation operator, P. Each aggregate has a “root-node” or “center-node” associated with it, and this root-node is injected from the coarse grid to the fine grid. The injection corresponds to the identity block.
  • Only smooth={‘energy’, None} is supported for prolongation smoothing. See reference [2] below for more details on why the ‘energy’ prolongation smoother is the natural counterpart to root-node style SA.
  • The additional parameters are passed through as arguments to multilevel_solver. Refer to pyamg.multilevel_solver for additional documentation.

  • At each level, four steps are executed in order to define the coarser level operator.

    1. Matrix A is given and used to derive a strength matrix, C.
    2. Based on the strength matrix, indices are grouped or aggregated.
    3. The aggregates define coarse nodes and a tentative prolongation operator T is defined by injection
    4. The tentative prolongation operator is smoothed by a relaxation scheme to improve the quality and extent of interpolation from the aggregates to fine nodes.
  • The parameters smooth, strength, aggregate, presmoother, postsmoother can be varied on a per level basis. For different methods on different levels, use a list as input so that the i-th entry defines the method at the i-th level. If there are more levels in the hierarchy than list entries, the last entry will define the method for all levels lower.

    Examples are: smooth=[(‘jacobi’, {‘omega’:1.0}), None, ‘jacobi’] presmoother=[(‘block_gauss_seidel’, {‘sweep’:symmetric}), ‘sor’] aggregate=[‘standard’, ‘naive’] strength=[(‘symmetric’, {‘theta’:0.25}),

    (‘symmetric’, {‘theta’:0.08})]

  • Predefined strength of connection and aggregation schemes can be specified. These options are best used together, but aggregation can be predefined while strength of connection is not.

    For predefined strength of connection, use a list consisting of tuples of the form (‘predefined’, {‘C’ : C0}), where C0 is a csr_matrix and each degree-of-freedom in C0 represents a supernode. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘C’ : C0}), (‘predefined’, {‘C’ : C1}) ].

    Similarly for predefined aggregation, use a list of tuples. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘AggOp’ : Agg0}), (‘predefined’, {‘AggOp’ : Agg1}) ], where the dimensions of A, Agg0 and Agg1 are compatible, i.e. Agg0.shape[1] == A.shape[0] and Agg1.shape[1] == Agg0.shape[0]. Each AggOp is a csr_matrix.

    Because this is a root-nodes solver, if a member of the predefined aggregation list is predefined, it must be of the form (‘predefined’, {‘AggOp’ : Agg, ‘Cnodes’ : Cnodes}).

>>> from pyamg import rootnode_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> import numpy
>>> A = poisson((100, 100), format='csr')           # matrix
>>> b = numpy.ones((A.shape[0]))                   # RHS
>>> ml = rootnode_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
[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
[2]Olson, L. and Schroder, J. and Tuminaro, R., “A general interpolation strategy for algebraic multigrid using energy minimization”, SIAM Journal on Scientific Computing (SISC), vol. 33, pp. 966–991, 2011.

setup Module

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

smooth Module

Methods to smooth tentative prolongation operators

pyamg.aggregation.smooth.jacobi_prolongation_smoother(S, T, C, B, omega=1.3333333333333333, degree=1, filter=False, weighting='diagonal')[source]

Jacobi prolongation smoother

S : {csr_matrix, bsr_matrix}
Sparse NxN matrix used for smoothing. Typically, A.
T : {csr_matrix, bsr_matrix}
Tentative prolongator
C : {csr_matrix, bsr_matrix}
Strength-of-connection matrix
B : {array}
Near nullspace modes for the coarse grid such that T*B exactly reproduces the fine grid near nullspace modes
omega : {scalar}
Damping parameter
filter : {boolean}
If true, filter S before smoothing T. This option can greatly control complexity.
weighting : {string}

‘block’, ‘diagonal’ or ‘local’ weighting for constructing the Jacobi D ‘local’: Uses a local row-wise weight based on the Gershgorin estimate.

Avoids any potential under-damping due to inaccurate spectral radius estimates.

‘block’: If A is a BSR matrix, use a block diagonal inverse of A ‘diagonal’: Classic Jacobi D = diagonal(A)

P : {csr_matrix, bsr_matrix}
Smoothed (final) prolongator defined by P = (I - omega/rho(K) K) * T where K = diag(S)^-1 * S and rho(K) is an approximation to the spectral radius of K.

If weighting is not ‘local’, then results using Jacobi prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.

>>> from pyamg.aggregation import jacobi_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1.,  0.],
        [ 1.,  0.],
        [ 1.,  0.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.]])
>>> A = poisson((6,),format='csr')
>>> P = jacobi_prolongation_smoother(A,T,A,numpy.ones((2,1)))
>>> P.todense()
matrix([[ 0.64930164,  0.        ],
        [ 1.        ,  0.        ],
        [ 0.64930164,  0.35069836],
        [ 0.35069836,  0.64930164],
        [ 0.        ,  1.        ],
        [ 0.        ,  0.64930164]])
pyamg.aggregation.smooth.richardson_prolongation_smoother(S, T, omega=1.3333333333333333, degree=1)[source]

Richardson prolongation smoother

S : {csr_matrix, bsr_matrix}
Sparse NxN matrix used for smoothing. Typically, A or the “filtered matrix” obtained from A by lumping weak connections onto the diagonal of A.
T : {csr_matrix, bsr_matrix}
Tentative prolongator
omega : {scalar}
Damping parameter
P : {csr_matrix, bsr_matrix}
Smoothed (final) prolongator defined by P = (I - omega/rho(S) S) * T where rho(S) is an approximation to the spectral radius of S.

Results using Richardson prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.

>>> from pyamg.aggregation import richardson_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1.,  0.],
        [ 1.,  0.],
        [ 1.,  0.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.]])
>>> A = poisson((6,),format='csr')
>>> P = richardson_prolongation_smoother(A,T)
>>> P.todense()
matrix([[ 0.64930164,  0.        ],
        [ 1.        ,  0.        ],
        [ 0.64930164,  0.35069836],
        [ 0.35069836,  0.64930164],
        [ 0.        ,  1.        ],
        [ 0.        ,  0.64930164]])
pyamg.aggregation.smooth.energy_prolongation_smoother(A, T, Atilde, B, Bf, Cpt_params, krylov='cg', maxiter=4, tol=1e-08, degree=1, weighting='local')[source]

Minimize the energy of the coarse basis functions (columns of T). Both root-node and non-root-node style prolongation smoothing is available, see Cpt_params description below.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix
T : {bsr_matrix}
Tentative prolongator, a NxM sparse matrix (M < N)
Atilde : {csr_matrix}
Strength of connection matrix
B : {array}
Near-nullspace modes for coarse grid. Has shape (M,k) where k is the number of coarse candidate vectors.
Bf : {array}
Near-nullspace modes for fine grid. Has shape (N,k) where k is the number of coarse candidate vectors.
Cpt_params : {tuple}
Tuple of the form (bool, dict). If the Cpt_params[0] = False, then the standard SA prolongation smoothing is carried out. If True, then root-node style prolongation smoothing is carried out. The dict must be a dictionary of parameters containing, (1) P_I: P_I.T is the injection matrix for the Cpts, (2) I_F: an identity matrix for only the F-points (i.e. I, but with zero rows and columns for C-points) and I_C: the C-point analogue to I_F. See Notes below for more information.
krylov : {string}

‘cg’ : for SPD systems. Solve A T = 0 in a constraint space with CG ‘cgnr’ : for nonsymmetric and/or indefinite systems.

Solve A T = 0 in a constraint space with CGNR
‘gmres’ : for nonsymmetric and/or indefinite systems.
Solve A T = 0 in a constraint space with GMRES
maxiter : integer
Number of energy minimization steps to apply to the prolongator
tol : {scalar}
Minimization tolerance
degree : {int}
Generate sparsity pattern for P based on (Atilde^degree T)
weighting : {string}
‘block’, ‘diagonal’ or ‘local’ construction of the
diagonal preconditioning
‘local’: Uses a local row-wise weight based on the Gershgorin estimate.
Avoids any potential under-damping due to inaccurate spectral radius estimates.

‘block’: If A is a BSR matrix, use a block diagonal inverse of A ‘diagonal’: Use inverse of the diagonal of A

T : {bsr_matrix}
Smoothed prolongator

Only ‘diagonal’ weighting is supported for the CGNR method, because we are working with A^* A and not A.

When Cpt_params[0] == True, root-node style prolongation smoothing is used to minimize the energy of columns of T. Essentially, an identity block is maintained in T, corresponding to injection from the coarse-grid to the fine-grid root-nodes. See [2] for more details, and see util.utils.get_Cpt_params for the helper function to generate Cpt_params.

If Cpt_params[0] == False, the energy of columns of T are still minimized, but without maintaining the identity block.

>>> from pyamg.aggregation import energy_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> print T.todense()
[[ 1.  0.]
 [ 1.  0.]
 [ 1.  0.]
 [ 0.  1.]
 [ 0.  1.]
 [ 0.  1.]]
>>> A = poisson((6,),format='csr')
>>> B = numpy.ones((2,1),dtype=float)
>>> P = energy_prolongation_smoother(A,T,A,B, None, (False,{}))
>>> print P.todense()
[[ 1.          0.        ]
 [ 1.          0.        ]
 [ 0.66666667  0.33333333]
 [ 0.33333333  0.66666667]
 [ 0.          1.        ]
 [ 0.          1.        ]]
[1]Jan Mandel, Marian Brezina, and Petr Vanek “Energy Optimization of Algebraic Multigrid Bases” Computing 62, 205-228, 1999 http://dx.doi.org/10.1007/s006070050022
[2]Olson, L. and Schroder, J. and Tuminaro, R., “A general interpolation strategy for algebraic multigrid using energy minimization”, SIAM Journal on Scientific Computing (SISC), vol. 33, pp. 966–991, 2011.

tentative Module

Tentative prolongator

pyamg.aggregation.tentative.fit_candidates(AggOp, B, tol=1e-10)[source]

Fit near-nullspace candidates to form the tentative prolongator

AggOp : csr_matrix
Describes the sparsity pattern of the tentative prolongator. Has dimension (#blocks, #aggregates)
B : array
The near-nullspace candidates stored in column-wise fashion. Has dimension (#blocks * blocksize, #candidates)
tol : scalar
Threshold for eliminating local basis functions. If after orthogonalization a local basis function Q[:, j] is small, i.e. ||Q[:, j]|| < tol, then Q[:, j] is set to zero.
(Q, R) : (bsr_matrix, array)
The tentative prolongator Q is a sparse block matrix with dimensions (#blocks * blocksize, #aggregates * #candidates) formed by dense blocks of size (blocksize, #candidates). The coarse level candidates are stored in R which has dimensions (#aggregates * #candidates, #candidates).

amg_core.fit_candidates

Assuming that each row of AggOp contains exactly one non-zero entry, i.e. all unknowns belong to an aggregate, then Q and R satisfy the relationship B = Q*R. In other words, the near-nullspace candidates are represented exactly by the tentative prolongator.

If AggOp contains rows with no non-zero entries, then the range of the tentative prolongator will not include those degrees of freedom. This situation is illustrated in the examples below.

[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
>>> from scipy.sparse import csr_matrix
>>> from pyamg.aggregation.tentative import fit_candidates
>>> # four nodes divided into two aggregates
... AggOp = csr_matrix( [[1, 0],
...                      [1, 0],
...                      [0, 1],
...                      [0, 1]] )
>>> # B contains one candidate, the constant vector
... B = [[1],
...      [1],
...      [1],
...      [1]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678,  0.        ],
        [ 0.70710678,  0.        ],
        [ 0.        ,  0.70710678],
        [ 0.        ,  0.70710678]])
>>> R
array([[ 1.41421356],
       [ 1.41421356]])
>>> # Two candidates, the constant vector and a linear function
... B = [[1, 0],
...      [1, 1],
...      [1, 2],
...      [1, 3]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678, -0.70710678,  0.        ,  0.        ],
        [ 0.70710678,  0.70710678,  0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.70710678, -0.70710678],
        [ 0.        ,  0.        ,  0.70710678,  0.70710678]])
>>> R
array([[ 1.41421356,  0.70710678],
       [ 0.        ,  0.70710678],
       [ 1.41421356,  3.53553391],
       [ 0.        ,  0.70710678]])
>>> # aggregation excludes the third node
... AggOp = csr_matrix( [[1, 0],
...                      [1, 0],
...                      [0, 0],
...                      [0, 1]] )
>>> B = [[1],
...      [1],
...      [1],
...      [1]]
>>> Q, R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678,  0.        ],
        [ 0.70710678,  0.        ],
        [ 0.        ,  0.        ],
        [ 0.        ,  1.        ]])
>>> R
array([[ 1.41421356],
       [ 1.        ]])