897 lines
27 KiB
Python
897 lines
27 KiB
Python
"""Abstract linear algebra library.
|
|
|
|
This module defines a class hierarchy that implements a kind of "lazy"
|
|
matrix representation, called the ``LinearOperator``. It can be used to do
|
|
linear algebra with extremely large sparse or structured matrices, without
|
|
representing those explicitly in memory. Such matrices can be added,
|
|
multiplied, transposed, etc.
|
|
|
|
As a motivating example, suppose you want have a matrix where almost all of
|
|
the elements have the value one. The standard sparse matrix representation
|
|
skips the storage of zeros, but not ones. By contrast, a LinearOperator is
|
|
able to represent such matrices efficiently. First, we need a compact way to
|
|
represent an all-ones matrix::
|
|
|
|
>>> import numpy as np
|
|
>>> from scipy.sparse.linalg._interface import LinearOperator
|
|
>>> class Ones(LinearOperator):
|
|
... def __init__(self, shape):
|
|
... super().__init__(dtype=None, shape=shape)
|
|
... def _matvec(self, x):
|
|
... return np.repeat(x.sum(), self.shape[0])
|
|
|
|
Instances of this class emulate ``np.ones(shape)``, but using a constant
|
|
amount of storage, independent of ``shape``. The ``_matvec`` method specifies
|
|
how this linear operator multiplies with (operates on) a vector. We can now
|
|
add this operator to a sparse matrix that stores only offsets from one::
|
|
|
|
>>> from scipy.sparse.linalg._interface import aslinearoperator
|
|
>>> from scipy.sparse import csr_matrix
|
|
>>> offsets = csr_matrix([[1, 0, 2], [0, -1, 0], [0, 0, 3]])
|
|
>>> A = aslinearoperator(offsets) + Ones(offsets.shape)
|
|
>>> A.dot([1, 2, 3])
|
|
array([13, 4, 15])
|
|
|
|
The result is the same as that given by its dense, explicitly-stored
|
|
counterpart::
|
|
|
|
>>> (np.ones(A.shape, A.dtype) + offsets.toarray()).dot([1, 2, 3])
|
|
array([13, 4, 15])
|
|
|
|
Several algorithms in the ``scipy.sparse`` library are able to operate on
|
|
``LinearOperator`` instances.
|
|
"""
|
|
|
|
import warnings
|
|
|
|
import numpy as np
|
|
|
|
from scipy.sparse import issparse
|
|
from scipy.sparse._sputils import isshape, isintlike, asmatrix, is_pydata_spmatrix
|
|
|
|
__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.
|
|
|
|
To construct a concrete LinearOperator, either pass appropriate
|
|
callables to the constructor of this class, or subclass it.
|
|
|
|
A subclass must implement either one of the methods ``_matvec``
|
|
and ``_matmat``, and the attributes/properties ``shape`` (pair of
|
|
integers) and ``dtype`` (may be None). It may call the ``__init__``
|
|
on this class to have these attributes validated. Implementing
|
|
``_matvec`` automatically implements ``_matmat`` (using a naive
|
|
algorithm) and vice-versa.
|
|
|
|
Optionally, a subclass may implement ``_rmatvec`` or ``_adjoint``
|
|
to implement the Hermitian adjoint (conjugate transpose). As with
|
|
``_matvec`` and ``_matmat``, implementing either ``_rmatvec`` or
|
|
``_adjoint`` implements the other automatically. Implementing
|
|
``_adjoint`` is preferable; ``_rmatvec`` is mostly there for
|
|
backwards compatibility.
|
|
|
|
Parameters
|
|
----------
|
|
shape : tuple
|
|
Matrix dimensions (M, N).
|
|
matvec : callable f(v)
|
|
Returns returns A * v.
|
|
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.
|
|
rmatmat : callable f(V)
|
|
Returns A^H * V, where V is a dense matrix with dimensions (M, K).
|
|
|
|
Attributes
|
|
----------
|
|
args : tuple
|
|
For linear operators describing products etc. of other linear
|
|
operators, the operands of the binary operation.
|
|
ndim : int
|
|
Number of dimensions (this is always 2)
|
|
|
|
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.
|
|
|
|
LinearOperator instances can also be multiplied, added with each
|
|
other and exponentiated, all lazily: the result of these operations
|
|
is always a new, composite LinearOperator, that defers linear
|
|
operations to the original operators and combines the results.
|
|
|
|
More details regarding how to subclass a LinearOperator and several
|
|
examples of concrete LinearOperator instances can be found in the
|
|
external project `PyLops <https://pylops.readthedocs.io>`_.
|
|
|
|
|
|
Examples
|
|
--------
|
|
>>> import numpy as np
|
|
>>> from scipy.sparse.linalg import LinearOperator
|
|
>>> def mv(v):
|
|
... return np.array([2*v[0], 3*v[1]])
|
|
...
|
|
>>> A = LinearOperator((2,2), matvec=mv)
|
|
>>> A
|
|
<2x2 _CustomLinearOperator with dtype=float64>
|
|
>>> A.matvec(np.ones(2))
|
|
array([ 2., 3.])
|
|
>>> A * np.ones(2)
|
|
array([ 2., 3.])
|
|
|
|
"""
|
|
|
|
ndim = 2
|
|
# Necessary for right matmul with numpy arrays.
|
|
__array_ufunc__ = None
|
|
|
|
def __new__(cls, *args, **kwargs):
|
|
if cls is LinearOperator:
|
|
# Operate as _CustomLinearOperator factory.
|
|
return super().__new__(_CustomLinearOperator)
|
|
else:
|
|
obj = super().__new__(cls)
|
|
|
|
if (type(obj)._matvec == LinearOperator._matvec
|
|
and type(obj)._matmat == LinearOperator._matmat):
|
|
warnings.warn("LinearOperator subclass should implement"
|
|
" at least one of _matvec and _matmat.",
|
|
category=RuntimeWarning, stacklevel=2)
|
|
|
|
return obj
|
|
|
|
def __init__(self, dtype, shape):
|
|
"""Initialize this LinearOperator.
|
|
|
|
To be called by subclasses. ``dtype`` may be None; ``shape`` should
|
|
be convertible to a length-2 tuple.
|
|
"""
|
|
if dtype is not None:
|
|
dtype = np.dtype(dtype)
|
|
|
|
shape = tuple(shape)
|
|
if not isshape(shape):
|
|
raise ValueError(f"invalid shape {shape!r} (must be 2-d)")
|
|
|
|
self.dtype = dtype
|
|
self.shape = shape
|
|
|
|
def _init_dtype(self):
|
|
"""Called from subclasses at the end of the __init__ routine.
|
|
"""
|
|
if self.dtype is None:
|
|
v = np.zeros(self.shape[-1])
|
|
self.dtype = np.asarray(self.matvec(v)).dtype
|
|
|
|
def _matmat(self, X):
|
|
"""Default matrix-matrix multiplication handler.
|
|
|
|
Falls back on the user-defined _matvec method, so defining that will
|
|
define matrix multiplication (though in a very suboptimal way).
|
|
"""
|
|
|
|
return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])
|
|
|
|
def _matvec(self, x):
|
|
"""Default matrix-vector multiplication handler.
|
|
|
|
If self is a linear operator of shape (M, N), then this method will
|
|
be called on a shape (N,) or (N, 1) ndarray, and should return a
|
|
shape (M,) or (M, 1) ndarray.
|
|
|
|
This default implementation falls back on _matmat, so defining that
|
|
will define matrix-vector multiplication as well.
|
|
"""
|
|
return self.matmat(x.reshape(-1, 1))
|
|
|
|
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 1-d 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 or overridden
|
|
_matvec method 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 = 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 rmatvec(self, x):
|
|
"""Adjoint matrix-vector multiplication.
|
|
|
|
Performs the operation y = A^H * x where A is an MxN linear
|
|
operator and x is a column vector or 1-d array.
|
|
|
|
Parameters
|
|
----------
|
|
x : {matrix, ndarray}
|
|
An array with shape (M,) or (M,1).
|
|
|
|
Returns
|
|
-------
|
|
y : {matrix, ndarray}
|
|
A matrix or ndarray with shape (N,) or (N,1) depending
|
|
on the type and shape of the x argument.
|
|
|
|
Notes
|
|
-----
|
|
This rmatvec wraps the user-specified rmatvec routine or overridden
|
|
_rmatvec method to ensure that y has the correct shape and type.
|
|
|
|
"""
|
|
|
|
x = np.asanyarray(x)
|
|
|
|
M,N = self.shape
|
|
|
|
if x.shape != (M,) and x.shape != (M,1):
|
|
raise ValueError('dimension mismatch')
|
|
|
|
y = self._rmatvec(x)
|
|
|
|
if isinstance(x, np.matrix):
|
|
y = asmatrix(y)
|
|
else:
|
|
y = np.asarray(y)
|
|
|
|
if x.ndim == 1:
|
|
y = y.reshape(N)
|
|
elif x.ndim == 2:
|
|
y = y.reshape(N,1)
|
|
else:
|
|
raise ValueError('invalid shape returned by user-defined rmatvec()')
|
|
|
|
return y
|
|
|
|
def _rmatvec(self, x):
|
|
"""Default implementation of _rmatvec; defers to adjoint."""
|
|
if type(self)._adjoint == LinearOperator._adjoint:
|
|
# _adjoint not overridden, prevent infinite recursion
|
|
raise NotImplementedError
|
|
else:
|
|
return self.H.matvec(x)
|
|
|
|
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 or overridden
|
|
_matmat method to ensure that y has the correct type.
|
|
|
|
"""
|
|
if not (issparse(X) or is_pydata_spmatrix(X)):
|
|
X = np.asanyarray(X)
|
|
|
|
if X.ndim != 2:
|
|
raise ValueError(f'expected 2-d ndarray or matrix, not {X.ndim}-d')
|
|
|
|
if X.shape[0] != self.shape[1]:
|
|
raise ValueError(f'dimension mismatch: {self.shape}, {X.shape}')
|
|
|
|
try:
|
|
Y = self._matmat(X)
|
|
except Exception as e:
|
|
if issparse(X) or is_pydata_spmatrix(X):
|
|
raise TypeError(
|
|
"Unable to multiply a LinearOperator with a sparse matrix."
|
|
" Wrap the matrix in aslinearoperator first."
|
|
) from e
|
|
raise
|
|
|
|
if isinstance(Y, np.matrix):
|
|
Y = asmatrix(Y)
|
|
|
|
return Y
|
|
|
|
def rmatmat(self, X):
|
|
"""Adjoint matrix-matrix multiplication.
|
|
|
|
Performs the operation y = A^H * x where A is an MxN linear
|
|
operator and x is a column vector or 1-d array, or 2-d array.
|
|
The default implementation defers to the adjoint.
|
|
|
|
Parameters
|
|
----------
|
|
X : {matrix, ndarray}
|
|
A matrix or 2D array.
|
|
|
|
Returns
|
|
-------
|
|
Y : {matrix, ndarray}
|
|
A matrix or 2D array depending on the type of the input.
|
|
|
|
Notes
|
|
-----
|
|
This rmatmat wraps the user-specified rmatmat routine.
|
|
|
|
"""
|
|
if not (issparse(X) or is_pydata_spmatrix(X)):
|
|
X = np.asanyarray(X)
|
|
|
|
if X.ndim != 2:
|
|
raise ValueError('expected 2-d ndarray or matrix, not %d-d'
|
|
% X.ndim)
|
|
|
|
if X.shape[0] != self.shape[0]:
|
|
raise ValueError(f'dimension mismatch: {self.shape}, {X.shape}')
|
|
|
|
try:
|
|
Y = self._rmatmat(X)
|
|
except Exception as e:
|
|
if issparse(X) or is_pydata_spmatrix(X):
|
|
raise TypeError(
|
|
"Unable to multiply a LinearOperator with a sparse matrix."
|
|
" Wrap the matrix in aslinearoperator() first."
|
|
) from e
|
|
raise
|
|
|
|
if isinstance(Y, np.matrix):
|
|
Y = asmatrix(Y)
|
|
return Y
|
|
|
|
def _rmatmat(self, X):
|
|
"""Default implementation of _rmatmat defers to rmatvec or adjoint."""
|
|
if type(self)._adjoint == LinearOperator._adjoint:
|
|
return np.hstack([self.rmatvec(col.reshape(-1, 1)) for col in X.T])
|
|
else:
|
|
return self.H.matmat(X)
|
|
|
|
def __call__(self, x):
|
|
return self*x
|
|
|
|
def __mul__(self, x):
|
|
return self.dot(x)
|
|
|
|
def __truediv__(self, other):
|
|
if not np.isscalar(other):
|
|
raise ValueError("Can only divide a linear operator by a scalar.")
|
|
|
|
return _ScaledLinearOperator(self, 1.0/other)
|
|
|
|
def dot(self, x):
|
|
"""Matrix-matrix or matrix-vector multiplication.
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like
|
|
1-d or 2-d array, representing a vector or matrix.
|
|
|
|
Returns
|
|
-------
|
|
Ax : array
|
|
1-d or 2-d array (depending on the shape of x) that represents
|
|
the result of applying this linear operator on x.
|
|
|
|
"""
|
|
if isinstance(x, LinearOperator):
|
|
return _ProductLinearOperator(self, x)
|
|
elif np.isscalar(x):
|
|
return _ScaledLinearOperator(self, x)
|
|
else:
|
|
if not issparse(x) and not is_pydata_spmatrix(x):
|
|
# Sparse matrices shouldn't be converted to numpy arrays.
|
|
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 1-d or 2-d array or matrix, got %r'
|
|
% x)
|
|
|
|
def __matmul__(self, other):
|
|
if np.isscalar(other):
|
|
raise ValueError("Scalar operands are not allowed, "
|
|
"use '*' instead")
|
|
return self.__mul__(other)
|
|
|
|
def __rmatmul__(self, other):
|
|
if np.isscalar(other):
|
|
raise ValueError("Scalar operands are not allowed, "
|
|
"use '*' instead")
|
|
return self.__rmul__(other)
|
|
|
|
def __rmul__(self, x):
|
|
if np.isscalar(x):
|
|
return _ScaledLinearOperator(self, x)
|
|
else:
|
|
return self._rdot(x)
|
|
|
|
def _rdot(self, x):
|
|
"""Matrix-matrix or matrix-vector multiplication from the right.
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like
|
|
1-d or 2-d array, representing a vector or matrix.
|
|
|
|
Returns
|
|
-------
|
|
xA : array
|
|
1-d or 2-d array (depending on the shape of x) that represents
|
|
the result of applying this linear operator on x from the right.
|
|
|
|
Notes
|
|
-----
|
|
This is copied from dot to implement right multiplication.
|
|
"""
|
|
if isinstance(x, LinearOperator):
|
|
return _ProductLinearOperator(x, self)
|
|
elif np.isscalar(x):
|
|
return _ScaledLinearOperator(self, x)
|
|
else:
|
|
if not issparse(x) and not is_pydata_spmatrix(x):
|
|
# Sparse matrices shouldn't be converted to numpy arrays.
|
|
x = np.asarray(x)
|
|
|
|
# We use transpose instead of rmatvec/rmatmat to avoid
|
|
# unnecessary complex conjugation if possible.
|
|
if x.ndim == 1 or x.ndim == 2 and x.shape[0] == 1:
|
|
return self.T.matvec(x.T).T
|
|
elif x.ndim == 2:
|
|
return self.T.matmat(x.T).T
|
|
else:
|
|
raise ValueError('expected 1-d or 2-d array or matrix, got %r'
|
|
% x)
|
|
|
|
def __pow__(self, p):
|
|
if np.isscalar(p):
|
|
return _PowerLinearOperator(self, p)
|
|
else:
|
|
return NotImplemented
|
|
|
|
def __add__(self, x):
|
|
if isinstance(x, LinearOperator):
|
|
return _SumLinearOperator(self, x)
|
|
else:
|
|
return NotImplemented
|
|
|
|
def __neg__(self):
|
|
return _ScaledLinearOperator(self, -1)
|
|
|
|
def __sub__(self, x):
|
|
return self.__add__(-x)
|
|
|
|
def __repr__(self):
|
|
M,N = self.shape
|
|
if self.dtype is None:
|
|
dt = 'unspecified dtype'
|
|
else:
|
|
dt = 'dtype=' + str(self.dtype)
|
|
|
|
return '<%dx%d %s with %s>' % (M, N, self.__class__.__name__, dt)
|
|
|
|
def adjoint(self):
|
|
"""Hermitian adjoint.
|
|
|
|
Returns the Hermitian adjoint of self, aka the Hermitian
|
|
conjugate or Hermitian transpose. For a complex matrix, the
|
|
Hermitian adjoint is equal to the conjugate transpose.
|
|
|
|
Can be abbreviated self.H instead of self.adjoint().
|
|
|
|
Returns
|
|
-------
|
|
A_H : LinearOperator
|
|
Hermitian adjoint of self.
|
|
"""
|
|
return self._adjoint()
|
|
|
|
H = property(adjoint)
|
|
|
|
def transpose(self):
|
|
"""Transpose this linear operator.
|
|
|
|
Returns a LinearOperator that represents the transpose of this one.
|
|
Can be abbreviated self.T instead of self.transpose().
|
|
"""
|
|
return self._transpose()
|
|
|
|
T = property(transpose)
|
|
|
|
def _adjoint(self):
|
|
"""Default implementation of _adjoint; defers to rmatvec."""
|
|
return _AdjointLinearOperator(self)
|
|
|
|
def _transpose(self):
|
|
""" Default implementation of _transpose; defers to rmatvec + conj"""
|
|
return _TransposedLinearOperator(self)
|
|
|
|
|
|
class _CustomLinearOperator(LinearOperator):
|
|
"""Linear operator defined in terms of user-specified operations."""
|
|
|
|
def __init__(self, shape, matvec, rmatvec=None, matmat=None,
|
|
dtype=None, rmatmat=None):
|
|
super().__init__(dtype, shape)
|
|
|
|
self.args = ()
|
|
|
|
self.__matvec_impl = matvec
|
|
self.__rmatvec_impl = rmatvec
|
|
self.__rmatmat_impl = rmatmat
|
|
self.__matmat_impl = matmat
|
|
|
|
self._init_dtype()
|
|
|
|
def _matmat(self, X):
|
|
if self.__matmat_impl is not None:
|
|
return self.__matmat_impl(X)
|
|
else:
|
|
return super()._matmat(X)
|
|
|
|
def _matvec(self, x):
|
|
return self.__matvec_impl(x)
|
|
|
|
def _rmatvec(self, x):
|
|
func = self.__rmatvec_impl
|
|
if func is None:
|
|
raise NotImplementedError("rmatvec is not defined")
|
|
return self.__rmatvec_impl(x)
|
|
|
|
def _rmatmat(self, X):
|
|
if self.__rmatmat_impl is not None:
|
|
return self.__rmatmat_impl(X)
|
|
else:
|
|
return super()._rmatmat(X)
|
|
|
|
def _adjoint(self):
|
|
return _CustomLinearOperator(shape=(self.shape[1], self.shape[0]),
|
|
matvec=self.__rmatvec_impl,
|
|
rmatvec=self.__matvec_impl,
|
|
matmat=self.__rmatmat_impl,
|
|
rmatmat=self.__matmat_impl,
|
|
dtype=self.dtype)
|
|
|
|
|
|
class _AdjointLinearOperator(LinearOperator):
|
|
"""Adjoint of arbitrary Linear Operator"""
|
|
|
|
def __init__(self, A):
|
|
shape = (A.shape[1], A.shape[0])
|
|
super().__init__(dtype=A.dtype, shape=shape)
|
|
self.A = A
|
|
self.args = (A,)
|
|
|
|
def _matvec(self, x):
|
|
return self.A._rmatvec(x)
|
|
|
|
def _rmatvec(self, x):
|
|
return self.A._matvec(x)
|
|
|
|
def _matmat(self, x):
|
|
return self.A._rmatmat(x)
|
|
|
|
def _rmatmat(self, x):
|
|
return self.A._matmat(x)
|
|
|
|
class _TransposedLinearOperator(LinearOperator):
|
|
"""Transposition of arbitrary Linear Operator"""
|
|
|
|
def __init__(self, A):
|
|
shape = (A.shape[1], A.shape[0])
|
|
super().__init__(dtype=A.dtype, shape=shape)
|
|
self.A = A
|
|
self.args = (A,)
|
|
|
|
def _matvec(self, x):
|
|
# NB. np.conj works also on sparse matrices
|
|
return np.conj(self.A._rmatvec(np.conj(x)))
|
|
|
|
def _rmatvec(self, x):
|
|
return np.conj(self.A._matvec(np.conj(x)))
|
|
|
|
def _matmat(self, x):
|
|
# NB. np.conj works also on sparse matrices
|
|
return np.conj(self.A._rmatmat(np.conj(x)))
|
|
|
|
def _rmatmat(self, x):
|
|
return np.conj(self.A._matmat(np.conj(x)))
|
|
|
|
def _get_dtype(operators, dtypes=None):
|
|
if dtypes is None:
|
|
dtypes = []
|
|
for obj in operators:
|
|
if obj is not None and hasattr(obj, 'dtype'):
|
|
dtypes.append(obj.dtype)
|
|
return np.result_type(*dtypes)
|
|
|
|
|
|
class _SumLinearOperator(LinearOperator):
|
|
def __init__(self, A, B):
|
|
if not isinstance(A, LinearOperator) or \
|
|
not isinstance(B, LinearOperator):
|
|
raise ValueError('both operands have to be a LinearOperator')
|
|
if A.shape != B.shape:
|
|
raise ValueError(f'cannot add {A} and {B}: shape mismatch')
|
|
self.args = (A, B)
|
|
super().__init__(_get_dtype([A, B]), A.shape)
|
|
|
|
def _matvec(self, x):
|
|
return self.args[0].matvec(x) + self.args[1].matvec(x)
|
|
|
|
def _rmatvec(self, x):
|
|
return self.args[0].rmatvec(x) + self.args[1].rmatvec(x)
|
|
|
|
def _rmatmat(self, x):
|
|
return self.args[0].rmatmat(x) + self.args[1].rmatmat(x)
|
|
|
|
def _matmat(self, x):
|
|
return self.args[0].matmat(x) + self.args[1].matmat(x)
|
|
|
|
def _adjoint(self):
|
|
A, B = self.args
|
|
return A.H + B.H
|
|
|
|
|
|
class _ProductLinearOperator(LinearOperator):
|
|
def __init__(self, A, B):
|
|
if not isinstance(A, LinearOperator) or \
|
|
not isinstance(B, LinearOperator):
|
|
raise ValueError('both operands have to be a LinearOperator')
|
|
if A.shape[1] != B.shape[0]:
|
|
raise ValueError(f'cannot multiply {A} and {B}: shape mismatch')
|
|
super().__init__(_get_dtype([A, B]),
|
|
(A.shape[0], B.shape[1]))
|
|
self.args = (A, B)
|
|
|
|
def _matvec(self, x):
|
|
return self.args[0].matvec(self.args[1].matvec(x))
|
|
|
|
def _rmatvec(self, x):
|
|
return self.args[1].rmatvec(self.args[0].rmatvec(x))
|
|
|
|
def _rmatmat(self, x):
|
|
return self.args[1].rmatmat(self.args[0].rmatmat(x))
|
|
|
|
def _matmat(self, x):
|
|
return self.args[0].matmat(self.args[1].matmat(x))
|
|
|
|
def _adjoint(self):
|
|
A, B = self.args
|
|
return B.H * A.H
|
|
|
|
|
|
class _ScaledLinearOperator(LinearOperator):
|
|
def __init__(self, A, alpha):
|
|
if not isinstance(A, LinearOperator):
|
|
raise ValueError('LinearOperator expected as A')
|
|
if not np.isscalar(alpha):
|
|
raise ValueError('scalar expected as alpha')
|
|
if isinstance(A, _ScaledLinearOperator):
|
|
A, alpha_original = A.args
|
|
# Avoid in-place multiplication so that we don't accidentally mutate
|
|
# the original prefactor.
|
|
alpha = alpha * alpha_original
|
|
|
|
dtype = _get_dtype([A], [type(alpha)])
|
|
super().__init__(dtype, A.shape)
|
|
self.args = (A, alpha)
|
|
|
|
def _matvec(self, x):
|
|
return self.args[1] * self.args[0].matvec(x)
|
|
|
|
def _rmatvec(self, x):
|
|
return np.conj(self.args[1]) * self.args[0].rmatvec(x)
|
|
|
|
def _rmatmat(self, x):
|
|
return np.conj(self.args[1]) * self.args[0].rmatmat(x)
|
|
|
|
def _matmat(self, x):
|
|
return self.args[1] * self.args[0].matmat(x)
|
|
|
|
def _adjoint(self):
|
|
A, alpha = self.args
|
|
return A.H * np.conj(alpha)
|
|
|
|
|
|
class _PowerLinearOperator(LinearOperator):
|
|
def __init__(self, A, p):
|
|
if not isinstance(A, LinearOperator):
|
|
raise ValueError('LinearOperator expected as A')
|
|
if A.shape[0] != A.shape[1]:
|
|
raise ValueError('square LinearOperator expected, got %r' % A)
|
|
if not isintlike(p) or p < 0:
|
|
raise ValueError('non-negative integer expected as p')
|
|
|
|
super().__init__(_get_dtype([A]), A.shape)
|
|
self.args = (A, p)
|
|
|
|
def _power(self, fun, x):
|
|
res = np.array(x, copy=True)
|
|
for i in range(self.args[1]):
|
|
res = fun(res)
|
|
return res
|
|
|
|
def _matvec(self, x):
|
|
return self._power(self.args[0].matvec, x)
|
|
|
|
def _rmatvec(self, x):
|
|
return self._power(self.args[0].rmatvec, x)
|
|
|
|
def _rmatmat(self, x):
|
|
return self._power(self.args[0].rmatmat, x)
|
|
|
|
def _matmat(self, x):
|
|
return self._power(self.args[0].matmat, x)
|
|
|
|
def _adjoint(self):
|
|
A, p = self.args
|
|
return A.H ** p
|
|
|
|
|
|
class MatrixLinearOperator(LinearOperator):
|
|
def __init__(self, A):
|
|
super().__init__(A.dtype, A.shape)
|
|
self.A = A
|
|
self.__adj = None
|
|
self.args = (A,)
|
|
|
|
def _matmat(self, X):
|
|
return self.A.dot(X)
|
|
|
|
def _adjoint(self):
|
|
if self.__adj is None:
|
|
self.__adj = _AdjointMatrixOperator(self)
|
|
return self.__adj
|
|
|
|
class _AdjointMatrixOperator(MatrixLinearOperator):
|
|
def __init__(self, adjoint):
|
|
self.A = adjoint.A.T.conj()
|
|
self.__adjoint = adjoint
|
|
self.args = (adjoint,)
|
|
self.shape = adjoint.shape[1], adjoint.shape[0]
|
|
|
|
@property
|
|
def dtype(self):
|
|
return self.__adjoint.dtype
|
|
|
|
def _adjoint(self):
|
|
return self.__adjoint
|
|
|
|
|
|
class IdentityOperator(LinearOperator):
|
|
def __init__(self, shape, dtype=None):
|
|
super().__init__(dtype, shape)
|
|
|
|
def _matvec(self, x):
|
|
return x
|
|
|
|
def _rmatvec(self, x):
|
|
return x
|
|
|
|
def _rmatmat(self, x):
|
|
return x
|
|
|
|
def _matmat(self, x):
|
|
return x
|
|
|
|
def _adjoint(self):
|
|
return self
|
|
|
|
|
|
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 additional information.
|
|
|
|
Notes
|
|
-----
|
|
If 'A' has no .dtype attribute, the data type is determined by calling
|
|
:func:`LinearOperator.matvec()` - set the .dtype attribute to prevent this
|
|
call upon the linear operator creation.
|
|
|
|
Examples
|
|
--------
|
|
>>> import numpy as np
|
|
>>> from scipy.sparse.linalg import aslinearoperator
|
|
>>> M = np.array([[1,2,3],[4,5,6]], dtype=np.int32)
|
|
>>> aslinearoperator(M)
|
|
<2x3 MatrixLinearOperator 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 ndim <= 2')
|
|
A = np.atleast_2d(np.asarray(A))
|
|
return MatrixLinearOperator(A)
|
|
|
|
elif issparse(A) or is_pydata_spmatrix(A):
|
|
return MatrixLinearOperator(A)
|
|
|
|
else:
|
|
if hasattr(A, 'shape') and hasattr(A, 'matvec'):
|
|
rmatvec = None
|
|
rmatmat = None
|
|
dtype = None
|
|
|
|
if hasattr(A, 'rmatvec'):
|
|
rmatvec = A.rmatvec
|
|
if hasattr(A, 'rmatmat'):
|
|
rmatmat = A.rmatmat
|
|
if hasattr(A, 'dtype'):
|
|
dtype = A.dtype
|
|
return LinearOperator(A.shape, A.matvec, rmatvec=rmatvec,
|
|
rmatmat=rmatmat, dtype=dtype)
|
|
|
|
else:
|
|
raise TypeError('type not understood')
|