pyamg.krylov#

Krylov solvers.

This module contains several Krylov subspace methods, in addition to two simple iterations, to solve linear systems iteratively. These methods often use multigrid as a preconditioner to accelerate convergence to the solution. See [1] and [2].

Functions#

gmres

fgmres

cgne

cgnr

cg

bicgstab

steepest descent, (simple iteration)

minimal residual (MR), (simple iteration)

References#

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 231-234, 2003 http://www-users.cs.umn.edu/~saad/books.html

[2]

Richard Barrett et al. “Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition”, SIAM http://www.netlib.org/linalg/html_templates/Templates.html http://www.netlib.org/templates/

pyamg.krylov.bicgstab(A, b, x0=None, tol=1e-05, criteria='rr', maxiter=None, M=None, callback=None, residuals=None)[source]#

Biconjugate Gradient Algorithm with Stabilization.

Solves the linear system Ax = b. Left preconditioning is supported.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria
criteriastring: Stopping criteria, let r=r_k, x=x_k ‘rr’: ||r|| < tol ||b|| ‘rr+’: ||r|| < tol (||b|| + ||A||_F ||x||) if ||b||=0, then set ||b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverted preconditioner, i.e. solve M A x = M b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: residual history in the 2-norm, including the initial residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 231-234, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import bicgstab
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = bicgstab(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
4.68163

pyamg.krylov.cg(A, b, x0=None, tol=1e-05, criteria='rr', maxiter=None, M=None, callback=None, residuals=None)[source]#

Conjugate Gradient algorithm.

Solves the linear system Ax = b. Left preconditioning is supported.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria
criteriastring: Stopping criteria, let r=r_k, x=x_k ‘rr’: ||r|| < tol ||b|| ‘rr+’: ||r|| < tol (||b|| + ||A||_F ||x||) ‘MrMr’: ||M r|| < tol ||M b|| ‘rMr’: <r, Mr>^1/2 < tol if ||b||=0, then set ||b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverse preconditioner, i.e. solve M A x = M b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: residual history in the 2-norm, including the initial residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 262-67, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import cg
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = cg(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
10.9371

pyamg.krylov.cgne(A, b, x0=None, tol=1e-05, criteria='rr', maxiter=None, M=None, callback=None, residuals=None)[source]#

Conjugate Gradient, Normal Error algorithm.

Applies CG to the normal equations, A A.H x = b. Left preconditioning is supported. Note that unless A is well-conditioned, the use of CGNE is inadvisable

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria
criteriastring: Stopping criteria, let r=r_k, x=x_k ‘rr’: ||r|| < tol ||b|| ‘rr+’: ||r|| < tol (||b|| + ||A||_F ||x||) ‘MrMr’: ||M r|| < tol ||M b|| ‘rMr’: <r, Mr>^1/2 < tol if ||b||=0, then set ||b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverted preconditioner, i.e. solve M A A.H x = M b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: residual history in the 2-norm, including the initial residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 276-7, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import cgne
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = cgne(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
46.1547

pyamg.krylov.cgnr(A, b, x0=None, tol=1e-05, criteria='rr', maxiter=None, M=None, callback=None, residuals=None)[source]#

Conjugate Gradient, Normal Residual algorithm.

Applies CG to the normal equations, A.H A x = A.H b. Left preconditioning is supported. Note that unless A is well-conditioned, the use of CGNR is inadvisable

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria
criteriastring: Stopping criteria, let r=r_k, x=x_k ‘rr’: ||r|| < tol ||b|| ‘rr+’: ||r|| < tol (||b|| + ||A||_F ||x||) ‘MrMr’: ||M r|| < tol ||M b|| ‘rMr’: <r, Mr>^1/2 < tol if ||b||=0, then set ||b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverted preconditioner, i.e. solve M A.H A x = M A.H b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: residual history in the 2-norm, including the initial residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 276-7, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import cgnr
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = cgnr(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
9.39102

pyamg.krylov.cr(A, b, x0=None, tol=1e-05, criteria='rr', maxiter=None, M=None, callback=None, residuals=None)[source]#

Conjugate Residual algorithm.

Solves the linear system Ax = b. Left preconditioning is supported. The matrix A must be Hermitian symmetric (but not necessarily definite).

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria
criteriastring: Stopping criteria, let r=r_k, x=x_k ‘rr’: ||r|| < tol ||b|| ‘rr+’: ||r|| < tol (||b|| + ||A||_F ||x||) ‘MrMr’: ||M r|| < tol ||M b|| if ||b||=0, then set ||b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverted preconditioner, i.e. solve M A x = M b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: residual history in the 2-norm, including the initial residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 262-67, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import cr
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = cr(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
6.54282

pyamg.krylov.fgmres(A, b, x0=None, tol=1e-05, restart=None, maxiter=None, M=None, callback=None, residuals=None, restrt=None)[source]#

Flexible Generalized Minimum Residual Method (fGMRES).

fGMRES iteratively refines the initial solution guess to the system Ax = b. fGMRES is flexible in the sense that the right preconditioner (M) can vary from iteration to iteration.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator

n x n, linear system to solve

barray, matrix

right hand side, shape is (n,) or (n,1)

x0array, matrix

initial guess, default is a vector of zeros

tolfloat

Tolerance for stopping criteria, let r=r_k ||r|| < tol ||b|| if ||b||=0, then set ||b||=1 for these tests.

restartNone, int

if int, restart is max number of inner iterations and maxiter is the max number of outer iterations
if None, do not restart GMRES, and max number of inner iterations is maxiter

maxiterNone, int

if restart is None, maxiter is the max number of inner iterations and GMRES does not restart
if restart is int, maxiter is the max number of outer iterations, and restart is the max number of inner iterations
defaults to min(n,40) if restart=None

Marray, matrix, sparse matrix, LinearOperator

n x n, inverted preconditioner, i.e. solve M A x = M b. M need not be stationary for fgmres

callbackfunction

User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector

residualslist

residual history in the 2-norm, including the initial residual

reorthboolean

If True, then a check is made whether to re-orthogonalize the Krylov space each GMRES iteration

restrtNone, int

Deprecated. See restart.

Returns:

xk, info

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

fGMRES allows for non-stationary preconditioners, as opposed to GMRES

For robustness, Householder reflections are used to orthonormalize the Krylov Space Givens Rotations are used to provide the residual norm each iteration Flexibility implies that the right preconditioner, M, can vary from iteration to iteration

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 151-172, pp. 272-275, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import fgmres
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = fgmres(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
6.54282

pyamg.krylov.gmres(A, b, x0=None, tol=1e-05, restart=None, maxiter=None, M=None, callback=None, residuals=None, orthog='householder', restrt=None, **kwargs)[source]#

Generalized Minimum Residual Method (GMRES).

GMRES iteratively refines the initial solution guess to the system Ax = b. Left preconditioned. Residuals are preconditioned residuals.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator

n x n, linear system to solve

barray, matrix

right hand side, shape is (n,) or (n,1)

x0array, matrix

initial guess, default is a vector of zeros

tolfloat

Tolerance for stopping criteria, let r=r_k ||M r|| < tol ||M b|| if ||b||=0, then set ||M b||=1 for these tests.

restartNone, int

if int, restart is max number of inner iterations and maxiter is the max number of outer iterations
if None, do not restart GMRES, and max number of inner iterations is maxiter

maxiterNone, int

if restart is None, maxiter is the max number of inner iterations and GMRES does not restart
if restart is int, maxiter is the max number of outer iterations, and restart is the max number of inner iterations
defaults to min(n,40) if restart=None

Marray, matrix, sparse matrix, LinearOperator

n x n, inverted preconditioner, i.e. solve M A x = M b.

callbackfunction

User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector

residualslist

preconditioned residual history in the 2-norm, including the initial residual

orthogstring

‘householder’ calls _gmres_householder which uses Householder reflections to find the orthogonal basis for the Krylov space. ‘mgs’ calls _gmres_mgs which uses modified Gram-Schmidt to find the orthogonal basis for the Krylov space

restrtNone, int

Deprecated. See restart.

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

The orthogonalization method, orthog=’householder’, is more robust than orthog=’mgs’, however for the majority of problems your problem will converge before ‘mgs’ loses orthogonality in your basis.

orthog=’householder’ has been more rigorously tested, and is therefore currently the default

The residual is the preconditioned residual.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 151-172, pp. 272-275, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import gmres
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = gmres(A,b, maxiter=2, tol=1e-8)
>>> print(f'{norm(b - A*x):.6}')
6.54282

pyamg.krylov.gmres_householder(A, b, x0=None, tol=1e-05, restart=None, maxiter=None, M=None, callback=None, residuals=None, restrt=None)[source]#

Generalized Minimum Residual Method (GMRES) based on Housholder.

GMRES iteratively refines the initial solution guess to the system Ax = b. Householder reflections are used for orthogonalization. Left preconditioning, leading to preconditioned residuals.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator

n x n, linear system to solve

barray, matrix

right hand side, shape is (n,) or (n,1)

x0array, matrix

initial guess, default is a vector of zeros

tolfloat

Tolerance for stopping criteria, let r=r_k ||M r|| < tol ||M b|| if ||b||=0, then set ||M b||=1 for these tests.

restartNone, int

if int, restart is max number of inner iterations and maxiter is the max number of outer iterations
if None, do not restart GMRES, and max number of inner iterations is maxiter

maxiterNone, int

if restart is None, maxiter is the max number of inner iterations and GMRES does not restart
if restart is int, maxiter is the max number of outer iterations, and restart is the max number of inner iterations
defaults to min(n,40) if restart=None

Marray, matrix, sparse matrix, LinearOperator

n x n, inverted preconditioner, i.e. solve M A x = M b.

callbackfunction

User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector

residualslist

preconditioned residual history in the 2-norm, including the initial preconditioned residual

reorthboolean

If True, then a check is made whether to re-orthogonalize the Krylov space each GMRES iteration

restrtNone, int

Deprecated. See restart.

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

For robustness, modified Gram-Schmidt is used to orthogonalize the Krylov Space Givens Rotations are used to provide the residual norm each iteration

The residual is the preconditioned residual.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 151-172, pp. 272-275, 2003 http://www-users.cs.umn.edu/~saad/books.html

Examples

>>> from pyamg.krylov import gmres
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10, 10))
>>> b = np.ones((A.shape[0],))
>>> (x, flag) = gmres(A, b, maxiter=2, tol=1e-8, orthog='householder')
>>> print(f'{norm(b - A*x):.6}')
6.54282

pyamg.krylov.gmres_mgs(A, b, x0=None, tol=1e-05, restart=None, maxiter=None, M=None, callback=None, residuals=None, reorth=False, restrt=None)[source]#

Generalized Minimum Residual Method (GMRES) based on MGS.

GMRES iteratively refines the initial solution guess to the system Ax = b. Modified Gram-Schmidt version. Left preconditioning, leading to preconditioned residuals.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator

n x n, linear system to solve

barray, matrix

right hand side, shape is (n,) or (n,1)

x0array, matrix

initial guess, default is a vector of zeros

tolfloat

Tolerance for stopping criteria, let r=r_k ||M r|| < tol ||M b|| if ||b||=0, then set ||M b||=1 for these tests.

restartNone, int

if int, restart is max number of inner iterations and maxiter is the max number of outer iterations
if None, do not restart GMRES, and max number of inner iterations is maxiter

maxiterNone, int

if restart is None, maxiter is the max number of inner iterations and GMRES does not restart
if restart is int, maxiter is the max number of outer iterations, and restart is the max number of inner iterations
defaults to min(n,40) if restart=None

Marray, matrix, sparse matrix, LinearOperator

n x n, inverted preconditioner, i.e. solve M A x = M b.

callbackfunction

User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector

residualslist

preconditioned residual history in the 2-norm, including the initial preconditioned residual

reorthboolean

If True, then a check is made whether to re-orthogonalize the Krylov space each GMRES iteration

restrtNone, int

Deprecated. See restart.

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input

Notes

The LinearOperator class is in scipy.sparse.linalg. Use this class if you prefer to define A or M as a mat-vec routine as opposed to explicitly constructing the matrix.

For robustness, modified Gram-Schmidt is used to orthogonalize the Krylov Space Givens Rotations are used to provide the residual norm each iteration

The residual is the preconditioned residual.

References

[1]

Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 151-172, pp. 272-275, 2003 http://www-users.cs.umn.edu/~saad/books.html

[2]

1. Kelley, http://www4.ncsu.edu/~ctk/matlab_roots.html

Examples

>>> from pyamg.krylov import gmres
>>> from pyamg.util.linalg import norm
>>> import numpy as np
>>> from pyamg.gallery import poisson
>>> A = poisson((10,10))
>>> b = np.ones((A.shape[0],))
>>> (x,flag) = gmres(A,b, maxiter=2, tol=1e-8, orthog='mgs')
>>> print(f'{norm(b - A*x):.6}')
6.54282

pyamg.krylov.minimal_residual(A, b, x0=None, tol=1e-05, maxiter=None, M=None, callback=None, residuals=None)[source]#

Minimal residual (MR) algorithm. 1D projection method.

Solves the linear system Ax = b. Left preconditioning is supported.

Parameters:

Aarray, matrix, sparse matrix, LinearOperator: n x n, linear system to solve
barray, matrix: right hand side, shape is (n,) or (n,1)
x0array, matrix: initial guess, default is a vector of zeros
tolfloat: Tolerance for stopping criteria, let r=r_k ||M r|| < tol ||M b|| if ||b||=0, then set ||M b||=1 for these tests.
maxiterint: maximum number of iterations allowed
Marray, matrix, sparse matrix, LinearOperator: n x n, inverted preconditioner, i.e. solve M A x = M b.
callbackfunction: User-supplied function is called after each iteration as callback(xk), where xk is the current solution vector
residualslist: preconditioned residual history in the 2-norm, including the initial preconditioned residual

Returns:

(xk, info)

xkan updated guess after k iterations to the solution of Ax = b

infohalting status

0	successful exit
>0	convergence to tolerance not achieved, return iteration count instead.
<0	numerical breakdown, or illegal input