sparse.linalg.matfuncs

Sparse matrix functions

Module Contents

Classes

MatrixPowerOperator(self,A,p,structure=None)
ProductOperator(self,*args,**kwargs) For now, this is limited to products of multiple square matrices.
_ExpmPadeHelper(self,A,structure=None,use_exact_onenorm=False) Help lazily evaluate a matrix exponential.

Functions

inv(A) Compute the inverse of a sparse matrix
_onenorm_matrix_power_nnm(A,p) Compute the 1-norm of a non-negative integer power of a non-negative matrix.
_onenorm(A)
_ident_like(A)
_is_upper_triangular(A)
_smart_matrix_product(A,B,alpha=None,structure=None) A matrix product that knows about sparse and structured matrices.
_onenormest_matrix_power(A,p,t=2,itmax=5,compute_v=False,compute_w=False,structure=None) Efficiently estimate the 1-norm of A^p.
_onenormest_product(operator_seq,t=2,itmax=5,compute_v=False,compute_w=False,structure=None) Efficiently estimate the 1-norm of the matrix product of the args.
expm(A) Compute the matrix exponential using Pade approximation.
_expm(A,use_exact_onenorm)
_solve_P_Q(U,V,structure=None) A helper function for expm_2009.
_sinch(x) Stably evaluate sinch.
_eq_10_42(lam_1,lam_2,t_12) Equation (10.42) of Functions of Matrices: Theory and Computation.
_fragment_2_1(X,T,s) A helper function for expm_2009.
_ell(A,m) A helper function for expm_2009.
inv(A)

Compute the inverse of a sparse matrix

A : (M,M) ndarray or sparse matrix
square matrix to be inverted
Ainv : (M,M) ndarray or sparse matrix
inverse of A

This computes the sparse inverse of A. If the inverse of A is expected to be non-sparse, it will likely be faster to convert A to dense and use scipy.linalg.inv.

>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import inv
>>> A = csc_matrix([[1., 0.], [1., 2.]])
>>> Ainv = inv(A)
>>> Ainv
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Column format>
>>> A.dot(Ainv)
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in Compressed Sparse Column format>
>>> A.dot(Ainv).todense()
matrix([[ 1.,  0.],
        [ 0.,  1.]])

New in version 0.12.0.

_onenorm_matrix_power_nnm(A, p)

Compute the 1-norm of a non-negative integer power of a non-negative matrix.

A : a square ndarray or matrix or sparse matrix
Input matrix with non-negative entries.
p : non-negative integer
The power to which the matrix is to be raised.
out : float
The 1-norm of the matrix power p of A.
_onenorm(A)
_ident_like(A)
_is_upper_triangular(A)
_smart_matrix_product(A, B, alpha=None, structure=None)

A matrix product that knows about sparse and structured matrices.

A : 2d ndarray
First matrix.
B : 2d ndarray
Second matrix.
alpha : float
The matrix product will be scaled by this constant.
structure : str, optional
A string describing the structure of both matrices A and B. Only upper_triangular is currently supported.
M : 2d ndarray
Matrix product of A and B.
class MatrixPowerOperator(A, p, structure=None)
__init__(A, p, structure=None)
_matvec(x)
_rmatvec(x)
_matmat(X)
T()
class ProductOperator(*args, **kwargs)

For now, this is limited to products of multiple square matrices.

__init__(*args, **kwargs)
_matvec(x)
_rmatvec(x)
_matmat(X)
T()
_onenormest_matrix_power(A, p, t=2, itmax=5, compute_v=False, compute_w=False, structure=None)

Efficiently estimate the 1-norm of A^p.

A : ndarray
Matrix whose 1-norm of a power is to be computed.
p : int
Non-negative integer power.
t : int, optional
A positive parameter controlling the tradeoff between accuracy versus time and memory usage. Larger values take longer and use more memory but give more accurate output.
itmax : int, optional
Use at most this many iterations.
compute_v : bool, optional
Request a norm-maximizing linear operator input vector if True.
compute_w : bool, optional
Request a norm-maximizing linear operator output vector if True.
est : float
An underestimate of the 1-norm of the sparse matrix.
v : ndarray, optional
The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm.
w : ndarray, optional
The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input.
_onenormest_product(operator_seq, t=2, itmax=5, compute_v=False, compute_w=False, structure=None)

Efficiently estimate the 1-norm of the matrix product of the args.

operator_seq : linear operator sequence
Matrices whose 1-norm of product is to be computed.
t : int, optional
A positive parameter controlling the tradeoff between accuracy versus time and memory usage. Larger values take longer and use more memory but give more accurate output.
itmax : int, optional
Use at most this many iterations.
compute_v : bool, optional
Request a norm-maximizing linear operator input vector if True.
compute_w : bool, optional
Request a norm-maximizing linear operator output vector if True.
structure : str, optional
A string describing the structure of all operators. Only upper_triangular is currently supported.
est : float
An underestimate of the 1-norm of the sparse matrix.
v : ndarray, optional
The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm.
w : ndarray, optional
The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input.
class _ExpmPadeHelper(A, structure=None, use_exact_onenorm=False)

Help lazily evaluate a matrix exponential.

The idea is to not do more work than we need for high expm precision, so we lazily compute matrix powers and store or precompute other properties of the matrix.

__init__(A, structure=None, use_exact_onenorm=False)

Initialize the object.

A : a dense or sparse square numpy matrix or ndarray
The matrix to be exponentiated.
structure : str, optional
A string describing the structure of matrix A. Only upper_triangular is currently supported.
use_exact_onenorm : bool, optional
If True then only the exact one-norm of matrix powers and products will be used. Otherwise, the one-norm of powers and products may initially be estimated.
A2()
A4()
A6()
A8()
A10()
d4_tight()
d6_tight()
d8_tight()
d10_tight()
d4_loose()
d6_loose()
d8_loose()
d10_loose()
pade3()
pade5()
pade7()
pade9()
pade13_scaled(s)
expm(A)

Compute the matrix exponential using Pade approximation.

A : (M,M) array_like or sparse matrix
2D Array or Matrix (sparse or dense) to be exponentiated
expA : (M,M) ndarray
Matrix exponential of A

This is algorithm (6.1) which is a simplification of algorithm (5.1).

New in version 0.12.0.

[1]Awad H. Al-Mohy and Nicholas J. Higham (2009) “A New Scaling and Squaring Algorithm for the Matrix Exponential.” SIAM Journal on Matrix Analysis and Applications. 31 (3). pp. 970-989. ISSN 1095-7162
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import expm
>>> A = csc_matrix([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
>>> A.todense()
matrix([[1, 0, 0],
        [0, 2, 0],
        [0, 0, 3]], dtype=int64)
>>> Aexp = expm(A)
>>> Aexp
<3x3 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Column format>
>>> Aexp.todense()
matrix([[  2.71828183,   0.        ,   0.        ],
        [  0.        ,   7.3890561 ,   0.        ],
        [  0.        ,   0.        ,  20.08553692]])
_expm(A, use_exact_onenorm)
_solve_P_Q(U, V, structure=None)

A helper function for expm_2009.

U : ndarray
Pade numerator.
V : ndarray
Pade denominator.
structure : str, optional
A string describing the structure of both matrices U and V. Only upper_triangular is currently supported.

The structure argument is inspired by similar args for theano and cvxopt functions.

_sinch(x)

Stably evaluate sinch.

The strategy of falling back to a sixth order Taylor expansion was suggested by the Spallation Neutron Source docs which was found on the internet by google search. http://www.ornl.gov/~t6p/resources/xal/javadoc/gov/sns/tools/math/ElementaryFunction.html The details of the cutoff point and the Horner-like evaluation was picked without reference to anything in particular.

Note that sinch is not currently implemented in scipy.special, whereas the “engineer’s” definition of sinc is implemented. The implementation of sinc involves a scaling factor of pi that distinguishes it from the “mathematician’s” version of sinc.

_eq_10_42(lam_1, lam_2, t_12)

Equation (10.42) of Functions of Matrices: Theory and Computation.

This is a helper function for _fragment_2_1 of expm_2009. Equation (10.42) is on page 251 in the section on Schur algorithms. In particular, section 10.4.3 explains the Schur-Parlett algorithm. expm([[lam_1, t_12], [0, lam_1]) = [[exp(lam_1), t_12*exp((lam_1 + lam_2)/2)*sinch((lam_1 - lam_2)/2)], [0, exp(lam_2)]

_fragment_2_1(X, T, s)

A helper function for expm_2009.

The argument X is modified in-place, but this modification is not the same as the returned value of the function. This function also takes pains to do things in ways that are compatible with sparse matrices, for example by avoiding fancy indexing and by using methods of the matrices whenever possible instead of using functions of the numpy or scipy libraries themselves.

_ell(A, m)

A helper function for expm_2009.

A : linear operator
A linear operator whose norm of power we care about.
m : int
The power of the linear operator
value : int
A value related to a bound.