import numpy as np
from scipy.sparse.sputils import isshape
from scipy.sparse import isspmatrix
__all__ = ['LinearOperator', 'aslinearoperator']
class LinearOperator:
"""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.
Optional Parameters
-------------------
rmatvec : callable f(v)
Returns A^H * v, where A^H is the conjugate transpose of A.
matmat : callable f(V)
Returns A * V, where V is a dense matrix with dimensions (N,K).
dtype : dtype
Data type of the matrix.
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.])
"""
def __init__(self, shape, matvec, rmatvec=None, matmat=None, dtype=None):
shape = tuple(shape)
if not isshape(shape):
raise ValueError('invalid shape')
self.shape = shape
self._matvec = matvec
if rmatvec is None:
def rmatvec(v):
raise NotImplementedError('rmatvec is not defined')
self.rmatvec = rmatvec
else:
self.rmatvec = rmatvec
if matmat is not None:
# matvec each column of V
self._matmat = matmat
if dtype is not None:
self.dtype = np.dtype(dtype)
def _matmat(self, X):
"""Default matrix-matrix multiplication handler. Falls back on
the user-defined matvec() routine, which is always provided.
"""
return np.hstack( [ self.matvec(col.reshape(-1,1)) for col in X.T ] )
def matvec(self, 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.
"""
x = np.asanyarray(x)
M,N = self.shape
if x.shape != (N,) and x.shape != (N,1):
raise ValueError('dimension mismatch')
y = self._matvec(x)
if isinstance(x, np.matrix):
y = np.asmatrix(y)
else:
y = np.asarray(y)
if x.ndim == 1:
y = y.reshape(M)
elif x.ndim == 2:
y = y.reshape(M,1)
else:
raise ValueError('invalid shape returned by user-defined matvec()')
return y
def matmat(self, 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.
"""
X = np.asanyarray(X)
if X.ndim != 2:
raise ValueError('expected rank-2 ndarray or matrix')
M,N = self.shape
if X.shape[0] != N:
raise ValueError('dimension mismatch')
Y = self._matmat(X)
if isinstance(Y, np.matrix):
Y = np.asmatrix(Y)
return Y
def __mul__(self,x):
x = np.asarray(x)
if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1:
return self.matvec(x)
elif x.ndim == 2:
return self.matmat(x)
else:
raise ValueError('expected rank-1 or rank-2 array or matrix')
def __repr__(self):
M,N = self.shape
if hasattr(self,'dtype'):
dt = 'dtype=' + str(self.dtype)
else:
dt = 'unspecified dtype'
return '<%dx%d LinearOperator with %s>' % (M,N,dt)
def 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>
"""
if isinstance(A, LinearOperator):
return A
elif isinstance(A, np.ndarray) or isinstance(A, np.matrix):
if A.ndim > 2:
raise ValueError('array must have rank <= 2')
A = np.atleast_2d(np.asarray(A))
def matvec(v):
return np.dot(A, v)
def rmatvec(v):
return np.dot(A.conj().transpose(), v)
def matmat(V):
return np.dot(A, V)
return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
matmat=matmat, dtype=A.dtype)
elif isspmatrix(A):
def matvec(v):
return A * v
def rmatvec(v):
return A.conj().transpose() * v
def matmat(V):
return A * V
return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
matmat=matmat, dtype=A.dtype)
else:
if hasattr(A, 'shape') and hasattr(A, 'matvec'):
rmatvec = None
dtype = None
if hasattr(A, 'rmatvec'):
rmatvec = A.rmatvec
if hasattr(A, 'dtype'):
dtype = A.dtype
return LinearOperator(A.shape, A.matvec,
rmatvec=rmatvec, dtype=dtype)
else:
raise TypeError('type not understood')
|