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

Warning

This documentation is work-in-progress and unorganized.

## Sparse Linear Algebra¶

The submodules of sparse.linalg:
1. eigen: sparse eigenvalue problem solvers
2. isolve: iterative methods for solving linear systems
3. 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.

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). 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). 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 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]

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 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]

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).

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. 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). 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).

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)

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: 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.

#### Previous topic

scipy.sparse.isspmatrix_dia