Sparse linear algebra (`scipy.sparse.linalg`)¶

Warning

This documentation is work-in-progress and unorganized.

Sparse Linear Algebra¶

The submodules of sparse.linalg:

eigen: sparse eigenvalue problem solvers
isolve: iterative methods for solving linear systems
dsolve: direct factorization methods for solving linear systems

Examples¶

class scipy.sparse.linalg.LinearOperator(shape, matvec, rmatvec=None, matmat=None, dtype=None)¶

Common interface for performing matrix vector products

Many iterative methods (e.g. cg, gmres) do not need to know the individual entries of a matrix to solve a linear system A*x=b. Such solvers only require the computation of matrix vector products, A*v where v is a dense vector. This class serves as an abstract interface between iterative solvers and matrix-like objects.

Parameters:

shape : tuple

Matrix dimensions (M,N)

matvec : callable f(v)

Returns returns A * v.

See also

aslinearoperator: Construct LinearOperators

Notes

The user-defined matvec() function must properly handle the case where v has shape (N,) as well as the (N,1) case. The shape of the return type is handled internally by LinearOperator.

Examples

>>> from scipy.sparse.linalg import LinearOperator
>>> from scipy import *
>>> def mv(v):
...     return array([ 2*v[0], 3*v[1]])
...
>>> A = LinearOperator( (2,2), matvec=mv )
>>> A
<2x2 LinearOperator with unspecified dtype>
>>> A.matvec( ones(2) )
array([ 2.,  3.])
>>> A * ones(2)
array([ 2.,  3.])

Methods

`matmat`(X)	Matrix-matrix multiplication
`matvec`(x)	Matrix-vector multiplication

matmat(X)¶

Matrix-matrix multiplication

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X : {matrix, ndarray}

An array with shape (N,K).

Returns:

Y : {matrix, ndarray}

A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Notes

This matmat wraps any user-specified matmat routine to ensure that y has the correct type.

matvec(x)¶

Matrix-vector multiplication

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or rank-1 array.

Parameters:

x : {matrix, ndarray}

An array with shape (N,) or (N,1).

Returns:

y : {matrix, ndarray}

A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Notes

This matvec wraps the user-specified matvec routine to ensure that y has the correct shape and type.

scipy.sparse.linalg.Tester¶: alias of NoseTester

scipy.sparse.linalg.aslinearoperator(A)¶

Return A as a LinearOperator.

‘A’ may be any of the following types:

ndarray
matrix
sparse matrix (e.g. csr_matrix, lil_matrix, etc.)
LinearOperator
An object with .shape and .matvec attributes

See the LinearOperator documentation for additonal information.

Examples

>>> from scipy import matrix
>>> M = matrix( [[1,2,3],[4,5,6]], dtype='int32' )
>>> aslinearoperator( M )
<2x3 LinearOperator with dtype=int32>

scipy.sparse.linalg.bicg(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M=None, callback=None)¶

Use BIConjugate Gradient iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

scipy.sparse.linalg.bicgstab(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M=None, callback=None)¶

Use BIConjugate Gradient STABilized iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

scipy.sparse.linalg.cg(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M=None, callback=None)¶

Use Conjugate Gradient iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

scipy.sparse.linalg.cgs(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M=None, callback=None)¶

Use Conjugate Gradient Squared iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

scipy.sparse.linalg.eigen(A, k=6, M=None, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True)¶

Find k eigenvalues and eigenvectors of the square matrix A.

Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for w[i] eigenvalues with corresponding eigenvectors x[i].

Parameters:

A : matrix, array, or object with matvec(x) method

An N x N matrix, array, or an object with matvec(x) method to perform the matrix vector product A * x. The sparse matrix formats in scipy.sparse are appropriate for A.

k : integer

The number of eigenvalues and eigenvectors desired

Returns:

w : array

Array of k eigenvalues

v : array

An array of k eigenvectors The v[i] is the eigenvector corresponding to the eigenvector w[i]

See also

eigen_symmetric: eigenvalues and eigenvectors for symmetric matrix A

scipy.sparse.linalg.eigen_symmetric(A, k=6, M=None, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True)¶

Find k eigenvalues and eigenvectors of the real symmetric square matrix A.

Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for w[i] eigenvalues with corresponding eigenvectors x[i].

Parameters:

A : matrix or array with real entries or object with matvec(x) method

An N x N real symmetric matrix or array or an object with matvec(x) method to perform the matrix vector product A * x. The sparse matrix formats in scipy.sparse are appropriate for A.

k : integer

The number of eigenvalues and eigenvectors desired

Returns:

w : array

Array of k eigenvalues

v : array

An array of k eigenvectors The v[i] is the eigenvector corresponding to the eigenvector w[i]

See also

eigen: eigenvalues and eigenvectors for a general (nonsymmetric) matrix A

scipy.sparse.linalg.factorized(A)¶

Return a fuction for solving a sparse linear system, with A pre-factorized.

Example:: solve = factorized( A ) # Makes LU decomposition. x1 = solve( rhs1 ) # Uses the LU factors. x2 = solve( rhs2 ) # Uses again the LU factors.

scipy.sparse.linalg.gmres(A, b, x0=None, tol=1.0000000000000001e-05, restrt=20, maxiter=None, xtype=None, M=None, callback=None)¶

Use Generalized Minimal RESidual iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

See also

LinearOperator

scipy.sparse.linalg.lgmres(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=1000, M=None, callback=None, inner_m=30, outer_k=3, outer_v=None, store_outer_Av=True)¶

Solve a matrix equation using the LGMRES algorithm.

The LGMRES algorithm [BJM] [BPh] is designed to avoid some problems in the convergence in restarted GMRES, and often converges in fewer iterations.

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

x0 : {array, matrix}

Starting guess for the solution.

tol : float

Tolerance to achieve. The algorithm terminates when either the relative or the absolute residual is below tol.

maxiter : integer

Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.

M : {sparse matrix, dense matrix, LinearOperator}

Preconditioner for A. The preconditioner should approximate the inverse of A. Effective preconditioning dramatically improves the rate of convergence, which implies that fewer iterations are needed to reach a given error tolerance.

callback : function

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

Returns:

x : array or matrix

The converged solution.

info : integer

Provides convergence information:

0 : successful exit >0 : convergence to tolerance not achieved, number of iterations <0 : illegal input or breakdown

Notes

The LGMRES algorithm [BJM] [BPh] is designed to avoid the slowing of convergence in restarted GMRES, due to alternating residual vectors. Typically, it often outperforms GMRES(m) of comparable memory requirements by some measure, or at least is not much worse.

Another advantage in this algorithm is that you can supply it with ‘guess’ vectors in the outer_v argument that augment the Krylov subspace. If the solution lies close to the span of these vectors, the algorithm converges faster. This can be useful if several very similar matrices need to be inverted one after another, such as in Newton-Krylov iteration where the Jacobian matrix often changes little in the nonlinear steps.

References

[BJM]

(1, 2, 3) A.H. Baker and E.R. Jessup and T. Manteuffel, SIAM J. Matrix Anal. Appl. 26, 962 (2005).

[BPh]

(1, 2, 3) A.H. Baker, PhD thesis, University of Colorado (2003). http://amath.colorado.edu/activities/thesis/allisonb/Thesis.ps

scipy.sparse.linalg.lobpcg(A, X, B=None, M=None, Y=None, tol=None, maxiter=20, largest=True, verbosityLevel=0, retLambdaHistory=False, retResidualNormsHistory=False)¶

Solve symmetric partial eigenproblems with optional preconditioning

This function implements the Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The symmetric linear operator of the problem, usually a sparse matrix. Often called the “stiffness matrix”.

X : array_like

Initial approximation to the k eigenvectors. If A has shape=(n,n) then X should have shape shape=(n,k).

Returns:

w : array

Array of k eigenvalues

v : array

An array of k eigenvectors. V has the same shape as X.

Notes

If both retLambdaHistory and retResidualNormsHistory are True, the return tuple has the following format (lambda, V, lambda history, residual norms history)

scipy.sparse.linalg.minres(A, b, x0=None, shift=0.0, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M=None, callback=None, show=False, check=False)¶

Use MINimum RESidual iteration to solve Ax=b

MINRES minimizes norm(A*x - b) for the symmetric matrix A. Unlike the Conjugate Gradient method, A can be indefinite or singular.

If shift != 0 then the method solves (A - shift*I)x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

Notes

THIS FUNCTION IS EXPERIMENTAL AND SUBJECT TO CHANGE!

References

Solution of sparse indefinite systems of linear equations,: C. C. Paige and M. A. Saunders (1975), SIAM J. Numer. Anal. 12(4), pp. 617-629. http://www.stanford.edu/group/SOL/software/minres.html
This file is a translation of the following MATLAB implementation:: http://www.stanford.edu/group/SOL/software/minres/matlab/

scipy.sparse.linalg.qmr(A, b, x0=None, tol=1.0000000000000001e-05, maxiter=None, xtype=None, M1=None, M2=None, callback=None)¶

Use Quasi-Minimal Residual iteration to solve A x = b

Parameters:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

See also

LinearOperator

scipy.sparse.linalg.splu(A, permc_spec=2, diag_pivot_thresh=1.0, drop_tol=0.0, relax=1, panel_size=10)¶

A linear solver, for a sparse, square matrix A, using LU decomposition where L is a lower triangular matrix and U is an upper triagular matrix.

Returns a factored_lu object. (scipy.sparse.linalg.dsolve._superlu.SciPyLUType)

See scipy.sparse.linalg.dsolve._superlu.dgstrf for more info.

scipy.sparse.linalg.spsolve(A, b, permc_spec=2)¶: Solve the sparse linear system Ax=b

scipy.sparse.linalg.test¶

Run tests for module using nose.

Parameters:

label : {‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the tests to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

‘fast’ - the default - which corresponds to the nosetests -A

option of ‘not slow’.

‘full’ - fast (as above) and slow tests as in the

‘no -A’ option to nosetests - this is the same as ‘’.

None or ‘’ - run all tests. attribute_identifier - string passed directly to nosetests as ‘-A’.

verbose : int, optional

Verbosity value for test outputs, in the range 1-10. Default is 1.

extra_argv : list, optional

List with any extra arguments to pass to nosetests.

doctests : bool, optional

If True, run doctests in module. Default is False.

coverage : bool, optional

If True, report coverage of NumPy code. Default is False. (This requires the `coverage module:

<http://nedbatchelder.com/code/modules/coverage.html>`_).

Returns:

result : object

Returns the result of running the tests as a nose.result.TextTestResult object.

Notes

Each NumPy module exposes test in its namespace to run all tests for it. For example, to run all tests for numpy.lib:

>>> np.lib.test()

Examples

>>> result = np.lib.test()
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s

OK

>>> result.errors
[]
>>> result.knownfail
[]

scipy.sparse.linalg.use_solver(**kwargs)¶

Valid keyword arguments with defaults (other ignored):: useUmfpack = True assumeSortedIndices = False

The default sparse solver is umfpack when available. This can be changed by passing useUmfpack = False, which then causes the always present SuperLU based solver to be used.

Umfpack requires a CSR/CSC matrix to have sorted column/row indices. If sure that the matrix fulfills this, pass assumeSortedIndices=True to gain some speed.

Sparse linear algebra (`scipy.sparse.linalg`)¶