from __future__ import annotations
__all__ = [
"LinearOperator",
"aslinearoperator",
]
import warnings
from abc import ABC, abstractmethod
from collections.abc import Callable
from concurrent.futures import ThreadPoolExecutor
from multiprocessing import Pool
from multiprocessing.pool import Pool as PoolClass
from typing import Any
import numpy as np
import scipy as sp
from numpy.linalg import solve as np_solve
from scipy.linalg import lstsq
from scipy.linalg import solve as sp_solve
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import LinearOperator as spLinearOperator
from scipy.sparse.linalg import eigs as sp_eigs
from scipy.sparse.linalg import lsqr, spsolve
# need to check scipy version since the interface submodule changed into
# _interface from scipy>=1.8.0
sp_version = sp.__version__.split(".")
if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8:
from scipy.sparse.sputils import isintlike, isshape
else:
from scipy.sparse._sputils import isintlike, isshape
from pylops import get_ndarray_multiplication
from pylops._multioperator import _matvec_rmatvec_map
from pylops.optimization.basic import cgls
from pylops.utils._internal import _get_dtype
from pylops.utils.backend import (
get_array_module,
get_eigsh,
get_eigvals,
get_lobpcg,
get_module,
get_sparse_eye,
)
from pylops.utils.decorators import count
from pylops.utils.estimators import trace_hutchinson, trace_hutchpp, trace_nahutchpp
from pylops.utils.typing import (
DTypeLike,
InputDimsLike,
NDArray,
ShapeLike,
Tbackend,
Tpool,
)
class _LinearOperator(ABC):
"""Meta-class for Linear operator"""
@abstractmethod
def _matvec(self, x: NDArray) -> NDArray:
"""Matrix-vector multiplication handler."""
pass
@abstractmethod
def _rmatvec(self, x: NDArray) -> NDArray:
"""Matrix-vector adjoint multiplication handler."""
pass
[docs]
class LinearOperator(_LinearOperator):
"""Common interface for performing matrix-vector products.
This class acts as an abstract interface between matrix-like
objects and iterative solvers, providing methods to perform
matrix-vector and adjoint matrix-vector products as as
well as convenience methods such as ``eigs``, ``cond``, and
``conj``.
.. note:: End users of PyLops should not use this class directly but simply
use operators that are already implemented. This class is meant for
developers and it has to be used as the parent class of any new operator
developed within PyLops. Find more details regarding implementation of
new operators at :ref:`addingoperator`.
Parameters
----------
Op : :obj:`scipy.sparse.linalg.LinearOperator` or :obj:`pylops.linearoperator.LinearOperator`
Operator. If other arguments are provided, they will overwrite those obtained from ``Op``.
dtype : :obj:`str` or :obj:`numpy.dtype`, optional
Type of elements in input array. Any ``numpy.dtype``-like input
(e.g., ``"float64"`` or ``np.float64``) is always cast to a
:obj:`numpy.dtype` via ``np.dtype(dtype)``, both at initialization
and when ``dtype`` is set afterwards. If ``None``, ``dtype`` is taken
from ``Op`` or, if ``Op`` is not provided, defaults to ``numpy.float64``
(following ``np.dtype(None)``); the same applies when setting ``dtype``
to ``None`` afterwards.
shape : :obj:`tuple(int, int)`, optional
Shape of operator. If not provided, obtained from ``dims`` and ``dimsd``.
dims : :obj:`tuple(int, ..., int)`, optional
.. versionadded:: 2.0.0
Dimensions of model. If not provided, ``(self.shape[1],)`` is used.
dimsd : :obj:`tuple(int, ..., int)`, optional
.. versionadded:: 2.0.0
Dimensions of data. If not provided, ``(self.shape[0],)`` is used.
clinear : :obj:`bool`, optional
.. versionadded:: 1.17.0
Operator is complex-linear. Defaults to ``True``.
explicit : :obj:`bool`, optional
Operator contains a matrix that can be solved explicitly
(``True``) or not (``False``). Defaults to ``False``.
forceflat : :obj:`bool`, optional
.. versionadded:: 2.2.0
Force an array to be flattened after matvec/rmatvec if the input is ambiguous
(i.e., is a 1D array both when operating with ND arrays and with 1D arrays).
Defaults to ``None`` for operators that have no ambiguity or to ``True``
for those with ambiguity.
name : :obj:`str`, optional
.. versionadded:: 2.0.0
Name of operator (to be used by :func:`pylops.utils.describe.describe`)
Attributes
----------
matvec_count : :obj:`int`
Counter tracking the number of ``matvec``
operations performed since the operator was instantiated.
rmatvec_count : :obj:`int`
Counter tracking the number of ``rmatvec``
operations performed since the operator was instantiated.
matmat_count : :obj:`int`
Counter tracking the number of ``matmat``
operations performed since the operator was instantiated.
rmatmat_count : :obj:`int`
Counter tracking the number of ``rmatmat``
operations performed since the operator was instantiated.
"""
def __init__(
self,
Op: spLinearOperator | LinearOperator | None = None,
dtype: DTypeLike | None = None,
shape: ShapeLike | None = None,
dims: ShapeLike | None = None,
dimsd: ShapeLike | None = None,
clinear: bool | None = None,
explicit: bool | None = None,
forceflat: bool | None = None,
name: str | None = None,
) -> None:
if Op is not None:
self.Op = Op
# All Operators must have shape and dtype
dtype = self.Op.dtype if dtype is None else dtype
shape = self.Op.shape if shape is None else shape
# Optionally, some operators have other attributes
dims = getattr(Op, "dims", (Op.shape[1],)) if dims is None else dims
dimsd = getattr(Op, "dimsd", (Op.shape[0],)) if dimsd is None else dimsd
clinear = getattr(Op, "clinear", True) if clinear is None else clinear
explicit = (
getattr(self.Op, "explicit", False) if explicit is None else explicit
)
if explicit and hasattr(Op, "A"):
self.A = Op.A
forceflat = (
getattr(self.Op, "forceflat", None) if forceflat is None else forceflat
)
name = getattr(Op, "name", None) if name is None else name
# set dtype to np.float64 when neither Op not dtype is provided
# (consistent with np.dtype(None))
self.dtype = np.float64 if dtype is None else dtype
# set other params
if shape is not None:
self.shape = shape
if dims is not None:
self.dims = dims
if dimsd is not None:
self.dimsd = dimsd
if clinear is not None:
self.clinear = clinear
if explicit is not None:
self.explicit = explicit
if forceflat is not None:
self.forceflat = forceflat
self.name = name
# ndim: required by some eigenvalue solvers
self.ndim = len(self.shape)
# counters
self.matvec_count = 0
self.rmatvec_count = 0
self.matmat_count = 0
self.rmatmat_count = 0
@property
def dtype(self) -> np.dtype:
return self._dtype
@dtype.setter
def dtype(self, new_dtype: DTypeLike) -> None:
# always store a numpy dtype
self._dtype = np.dtype(new_dtype)
@dtype.deleter
def dtype(self) -> None:
del self._dtype
@property
def shape(self):
_shape = getattr(self, "_shape", None)
if _shape is None: # Cannot find shape, falling back on dims and dimsd
dims = getattr(self, "_dims", None)
dimsd = getattr(self, "_dimsd", None)
if dims is None or dimsd is None: # Cannot find both dims and dimsd, error
msg = (
f"'{self.__class__.__name__}' object has no attribute 'shape' "
"nor both fallback attributes ('dims', 'dimsd')"
)
raise AttributeError(msg)
_shape = (int(np.prod(dimsd)), int(np.prod(dims)))
self._shape = _shape # Update to not redo everything above on next call
return _shape
@shape.setter
def shape(self, new_shape: ShapeLike) -> None:
new_shape = tuple(new_shape)
if not isshape(new_shape):
msg = f"Invalid shape; must be 2-d tuple of integers, got {new_shape}"
raise ValueError(msg)
dims = getattr(self, "_dims", None)
dimsd = getattr(self, "_dimsd", None)
if dims is not None and dimsd is not None: # Found dims and dimsd
if np.prod(dimsd) != new_shape[0] and np.prod(dims) != new_shape[1]:
msg = "New shape incompatible with dims and dimsd"
raise ValueError(msg)
elif np.prod(dimsd) != new_shape[0]:
msg = "New shape incompatible with dimsd"
raise ValueError(msg)
elif np.prod(dims) != new_shape[1]:
msg = "New shape incompatible with dims"
raise ValueError(msg)
self._shape = new_shape
@shape.deleter
def shape(self):
del self._shape
@property
def dims(self):
_dims = getattr(self, "_dims", None)
if _dims is None:
shape = getattr(self, "_shape", None)
if shape is None:
msg = (
f"'{self.__class__.__name__}' object has no "
"attributes 'dims' or 'shape'"
)
raise AttributeError(msg)
_dims = (shape[1],)
return _dims
@dims.setter
def dims(self, new_dims: ShapeLike) -> None:
new_dims = tuple(new_dims)
shape = getattr(self, "_shape", None)
if shape is None: # shape not set yet
self._dims = new_dims
else:
if np.prod(new_dims) == self.shape[1]:
self._dims = new_dims
else:
msg = "dims incompatible with shape[1]"
raise ValueError(msg)
@dims.deleter
def dims(self):
del self._dims
@property
def dimsd(self):
_dimsd = getattr(self, "_dimsd", None)
if _dimsd is None:
shape = getattr(self, "_shape", None)
if shape is None:
msg = (
f"'{self.__class__.__name__}' object has "
"no attributes 'dimsd' or 'shape'"
)
raise AttributeError(msg)
_dimsd = (shape[0],)
return _dimsd
@dimsd.setter
def dimsd(self, new_dimsd: ShapeLike) -> None:
new_dimsd = tuple(new_dimsd)
shape = getattr(self, "_shape", None)
if shape is None: # shape not set yet
self._dimsd = new_dimsd
else:
if np.prod(new_dimsd) == self.shape[0]:
self._dimsd = new_dimsd
else:
msg = "dimsd incompatible with shape[0]"
raise ValueError(msg)
@dimsd.deleter
def dimsd(self):
del self._dimsd
@property
def clinear(self):
return getattr(self, "_clinear", True)
@clinear.setter
def clinear(self, new_clinear: bool) -> None:
self._clinear = bool(new_clinear)
@clinear.deleter
def clinear(self):
del self._clinear
@property
def explicit(self):
return getattr(self, "_explicit", False)
@explicit.setter
def explicit(self, new_explicit: bool) -> None:
self._explicit = bool(new_explicit)
@explicit.deleter
def explicit(self):
del self._explicit
@property
def forceflat(self):
return getattr(self, "_forceflat", None)
@forceflat.setter
def forceflat(self, new_forceflat: bool) -> None:
# note that this can also be None so we check before forcing bool
self._forceflat = (
new_forceflat if new_forceflat is None else bool(new_forceflat)
)
@forceflat.deleter
def forceflat(self):
del self._forceflat
@property
def name(self):
return getattr(self, "_name", None)
@name.setter
def name(self, new_name: str) -> None:
self._name = new_name
@name.deleter
def name(self):
del self._name
def __mul__(self, x: float | LinearOperator) -> LinearOperator:
return self.dot(x)
def __matmul__(self, other):
if np.isscalar(other):
msg = "Scalar not allowed, use * instead"
raise ValueError(msg)
return self.__mul__(other)
def __rmul__(self, x: float) -> LinearOperator:
if np.isscalar(x):
Op = _ScaledLinearOperator(self, x)
self._copy_attributes(
Op,
exclude=[
"explicit",
"name",
],
)
Op.explicit = False
return Op
else:
return NotImplemented
def __rmatmul__(self, other):
if np.isscalar(other):
msg = "Scalar not allowed, use * instead"
raise ValueError(msg)
return self.__rmul__(other)
def __pow__(self, p: int) -> LinearOperator:
if np.isscalar(p):
Op = _PowerLinearOperator(self, p)
self._copy_attributes(
Op,
exclude=[
"explicit",
"name",
],
)
Op.explicit = False
return Op
else:
return NotImplemented
def __add__(self, x: LinearOperator) -> LinearOperator:
if isinstance(x, (LinearOperator, spLinearOperator)):
# cast x to pylops linear operator if not already (this is done
# to allow mixing pylops and scipy operators)
Opx = aslinearoperator(x)
Op = _SumLinearOperator(self, Opx)
self._copy_attributes(
Op,
exclude=[
"explicit",
"forceflat",
"name",
],
)
Op.clinear = Op.clinear and Opx.clinear
Op.explicit = False
if self.forceflat is None and Opx.forceflat is None:
Op.forceflat = None
elif self.forceflat is not None and Opx.forceflat is not None:
# Define forceflat only if differing, otherwise raise error
if self.forceflat != Opx.forceflat:
msg = f"Operators have conflicting forceflat {Op.forceflat} != {Opx.forceflat}"
raise ValueError(msg)
Op.forceflat = self.forceflat
else: # Only one of them is None
Op.forceflat = (
self.forceflat if self.forceflat is not None else Opx.forceflat
)
# Replace if shape-like
if len(self.dims) == 1:
Op.dims = Opx.dims
if len(self.dimsd) == 1:
Op.dimsd = Opx.dimsd
return Op
else:
return NotImplemented
def __neg__(self) -> LinearOperator:
Op = _ScaledLinearOperator(self, -1)
self._copy_attributes(
Op,
exclude=[
"explicit",
"name",
],
)
Op.explicit = False
return Op
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 _copy_attributes(
self,
dest: LinearOperator,
exclude: list[str] | None = None,
) -> None:
"""Copy attributes from one LinearOperator to another"""
if exclude is None:
exclude = ["name"]
attrs = ["dims", "dimsd", "clinear", "explicit", "forceflat", "name"]
if exclude is not None:
for item in exclude:
attrs.remove(item)
for attr in attrs:
if hasattr(self, attr):
setattr(dest, attr, getattr(self, attr))
def _matvec(self, x: NDArray) -> NDArray:
"""Matrix-vector multiplication handler."""
if self.Op is not None:
return self.Op._matvec(x)
def _rmatvec(self, x: NDArray) -> NDArray:
"""Matrix-vector adjoint multiplication handler."""
if self.Op is not None:
return self.Op._rmatvec(x)
def _matmat_serial(self, X: NDArray) -> NDArray:
"""Matrix-matrix multiplication (serial version)"""
ncp = get_array_module(X)
if sp.sparse.issparse(X):
y = ncp.vstack([self._matvec(col.toarray().reshape(-1)) for col in X.T]).T
else:
y = ncp.vstack([self._matvec(col.reshape(-1)) for col in X.T]).T
return y
def _matmat_multithread(self, X: NDArray, pool: ThreadPoolExecutor) -> NDArray:
"""Matrix-matrix multiplication (multithreaded version)"""
if sp.sparse.issparse(X):
ys = list(
pool.map(
lambda args: _matvec_rmatvec_map(*args),
[(self._matvec, col.toarray().reshape(-1)) for col in X.T],
)
)
else:
ys = list(
pool.map(
lambda args: _matvec_rmatvec_map(*args),
[(self._matvec, col.reshape(-1)) for col in X.T],
)
)
y = np.vstack(ys).T
return y
def _matmat_multiproc(self, X: NDArray, pool: Pool) -> NDArray:
"""Matrix-matrix multiplication (multiprocess version)"""
if sp.sparse.issparse(X):
ys = pool.starmap(
_matvec_rmatvec_map,
[(self._matvec, col.toarray().reshape(-1)) for col in X.T],
)
else:
ys = pool.starmap(
_matvec_rmatvec_map,
[(self._matvec, col.reshape(-1)) for col in X.T],
)
y = np.vstack(ys).T
return y
def _matmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
"""Matrix-matrix multiplication handler.
Modified version of scipy _matmat to avoid having trailing dimension
in col when provided to matvec
"""
if pool is None:
# serial
return self._matmat_serial(X)
elif isinstance(pool, ThreadPoolExecutor):
# multithread
return self._matmat_multithread(X, pool)
elif isinstance(pool, PoolClass):
# multiprocess
return self._matmat_multiproc(X, pool)
else:
msg = f"Received pool of unsupported type ({type(pool)})"
raise NotImplementedError(msg)
def _rmatmat_serial(self, X: NDArray) -> NDArray:
"""Matrix-matrix adjoint multiplication (serial version)"""
ncp = get_array_module(X)
if sp.sparse.issparse(X):
y = ncp.vstack([self._rmatvec(col.toarray().reshape(-1)) for col in X.T]).T
else:
y = ncp.vstack([self._rmatvec(col.reshape(-1)) for col in X.T]).T
return y
def _rmatmat_multithread(self, X: NDArray, pool: ThreadPoolExecutor) -> NDArray:
"""Matrix-matrix adjoint multiplication (multithreaded version)"""
if sp.sparse.issparse(X):
ys = list(
pool.map(
lambda args: _matvec_rmatvec_map(*args),
[(self._rmatvec, col.toarray().reshape(-1)) for col in X.T],
)
)
else:
ys = list(
pool.map(
lambda args: _matvec_rmatvec_map(*args),
[(self._rmatvec, col.reshape(-1)) for col in X.T],
)
)
y = np.vstack(ys).T
return y
def _rmatmat_multiproc(self, X: NDArray, pool: Pool) -> NDArray:
"""Matrix-matrix adjoint multiplication (multiprocess version)"""
if sp.sparse.issparse(X):
ys = pool.starmap(
_matvec_rmatvec_map,
[(self._rmatvec, col.toarray().reshape(-1)) for col in X.T],
)
else:
ys = pool.starmap(
_matvec_rmatvec_map,
[(self._rmatvec, col.reshape(-1)) for col in X.T],
)
y = np.vstack(ys).T
return y
def _rmatmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
"""Matrix-matrix adjoint multiplication handler.
Modified version of scipy _rmatmat to avoid having trailing dimension
in col when provided to rmatvec
"""
if pool is None:
# serial
return self._rmatmat_serial(X)
elif isinstance(pool, ThreadPoolExecutor):
# multithread
return self._rmatmat_multithread(X, pool)
elif isinstance(pool, PoolClass):
# multiprocess
return self._rmatmat_multiproc(X, pool)
else:
msg = f"Received pool of unsupported type ({type(pool)})"
raise NotImplementedError(msg)
def _adjoint(self) -> LinearOperator:
Op = _AdjointLinearOperator(self)
self._copy_attributes(Op, exclude=["dims", "dimsd", "explicit", "name"])
Op.explicit = False
Op.dims = self.dimsd
Op.dimsd = self.dims
return Op
def _transpose(self) -> LinearOperator:
Op = _TransposedLinearOperator(self)
self._copy_attributes(Op, exclude=["dims", "dimsd", "explicit", "name"])
Op.explicit = False
Op.dims = self.dimsd
Op.dimsd = self.dims
return Op
[docs]
def adjoint(self):
return self._adjoint()
H: Callable[[LinearOperator], LinearOperator] = property(adjoint)
[docs]
def transpose(self):
return self._transpose()
T: Callable[[LinearOperator], LinearOperator] = property(transpose)
[docs]
@count(forward=True)
def matvec(self, x: NDArray) -> NDArray:
"""Matrix-vector multiplication.
Modified version of scipy matvec which does not consider the case
where the input vector is ``np.matrix`` (the use ``np.matrix`` is now
discouraged in numpy's documentation).
Parameters
----------
x : :obj:`numpy.ndarray`
Input array of shape (N,) or (N,1)
Returns
-------
y : :obj:`numpy.ndarray`
Output array of shape (M,) or (M,1)
"""
M, N = self.shape
if x.shape != (N,) and x.shape != (N, 1):
msg = f"Dimension mismatch. Got {x.shape}, but expected ({N},) or ({N}, 1)."
raise ValueError(msg)
y = self._matvec(x)
if x.ndim == 1:
y = y.reshape(M)
elif x.ndim == 2:
y = y.reshape(M, 1)
else:
msg = "Invalid shape returned by user-defined matvec()"
raise ValueError(msg)
return y
[docs]
@count(forward=False)
def rmatvec(self, x: NDArray) -> NDArray:
"""Adjoint matrix-vector multiplication.
Modified version of scipy rmatvec which does not consider the case
where the input vector is ``np.matrix`` (the use ``np.matrix`` is now
discouraged in numpy's documentation).
Parameters
----------
y : :obj:`numpy.ndarray`
Input array of shape (M,) or (M,1)
Returns
-------
x : :obj:`numpy.ndarray`
Output array of shape (N,) or (N,1)
"""
M, N = self.shape
if x.shape != (M,) and x.shape != (M, 1):
msg = f"Dimension mismatch. Got {x.shape}, but expected ({M},) or ({M}, 1)."
raise ValueError(msg)
y = self._rmatvec(x)
if x.ndim == 1:
y = y.reshape(N)
elif x.ndim == 2:
y = y.reshape(N, 1)
else:
msg = "Invalid shape returned by user-defined rmatvec()"
raise ValueError(msg)
return y
[docs]
@count(forward=True, matmat=True)
def matmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
"""Matrix-matrix multiplication.
Modified version of scipy matmat which does not consider the case
where the input vector is ``np.matrix`` (the use ``np.matrix`` is now
discouraged in numpy's documentation).
Parameters
----------
x : :obj:`numpy.ndarray`
Input array of shape (N,K)
pool : :obj:`multiprocessing.Pool` or :obj:`concurrent.futures.ThreadPoolExecutor` or :obj:`None`
Pool of workers used to evaluate the operator on different columns of ``X`` in parallel.
Returns
-------
y : :obj:`numpy.ndarray`
Output array of shape (M,K)
"""
if X.ndim != 2:
msg = f"Expected 2-d ndarray or matrix, not {X.ndim}-d ndarray"
raise ValueError(msg)
if X.shape[0] != self.shape[1]:
msg = f"Dimension mismatch: {self.shape}, {X.shape}"
raise ValueError(msg)
Y = self._matmat(X, pool=pool)
return Y
[docs]
@count(forward=False, matmat=True)
def rmatmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
"""Matrix-matrix multiplication.
Modified version of scipy rmatmat which does not consider the case
where the input vector is ``np.matrix`` (the use ``np.matrix`` is now
discouraged in numpy's documentation).
Parameters
----------
x : :obj:`numpy.ndarray`
Input array of shape (M,K)
pool : :obj:`multiprocessing.Pool` or :obj:`concurrent.futures.ThreadPoolExecutor` or :obj:`None`
Pool of workers used to evaluate the operator on different columns of ``X`` in parallel.
Returns
-------
y : :obj:`numpy.ndarray`
Output array of shape (N,K)
"""
if X.ndim != 2:
msg = f"Expected 2-d ndarray or matrix, not {X.ndim}-d ndarray"
raise ValueError(msg)
if X.shape[0] != self.shape[0]:
msg = f"Dimension mismatch: {self.shape}, {X.shape}"
raise ValueError(msg)
Y = self._rmatmat(X, pool=pool)
return Y
[docs]
def dot(self, x: float | NDArray | LinearOperator) -> NDArray | LinearOperator:
"""Dot-product
Dot-product between a PyLops operator and a:
- scalar: returns a scaled operator, i.e., ``alpha * Op``;
- 1d- or nd-array: applies the operator to the input and returns
an array. A 1d-array of size ``shape[1]`` (or an nd-array of shape
``dims`` when ``dims`` is defined and ``forceflat`` is not set) is
treated as a vector and passed to :func:`pylops.LinearOperator.matvec`;
a 2d-array of shape ``(shape[1], P)`` (or an nd-array of shape
``(*dims, P)``) is treated as a matrix of ``P`` vectors and passed to
:func:`pylops.LinearOperator.matmat`. Arrays with more than 2
dimensions are accepted only if nd-array multiplication is enabled
via :func:`pylops.set_ndarray_multiplication`, in which case the
output is reshaped to ``dimsd`` (or ``(*dimsd, P)``);
- linear operator: returns the product of the two
operators, i.e., ``Op1 @ Op2``.
Parameters
----------
x : :obj:`float`, :obj:`numpy.ndarray` or :obj:`pylops.LinearOperator`
Input scalar, array (or matrix) or linear operator
Returns
-------
y : :obj:`numpy.ndarray` or :obj:`pylops.LinearOperator`
Scaled operator (scalar input), array (or matrix) that represents
the result of applying the linear operator on ``x`` (array input),
or product operator (operator input).
Raises
------
ValueError
If ``x`` is an array with an unsupported shape, or if the two
operators have conflicting ``forceflat``.
"""
if isinstance(x, (LinearOperator, spLinearOperator)):
# cast x to pylops linear operator if not already (this is done
# to allow mixing pylops and scipy operators)
Opx = aslinearoperator(x)
Op = _ProductLinearOperator(self, Opx)
self._copy_attributes(Op, exclude=["dims", "explicit", "forceflat", "name"])
Op.clinear = Op.clinear and Opx.clinear
Op.explicit = False
if self.forceflat is None and Opx.forceflat is None:
Op.forceflat = None
elif self.forceflat is not None and Opx.forceflat is not None:
# Define forceflat only if differing, otherwise raise error
if self.forceflat != Opx.forceflat:
msg = f"Operators have conflicting forceflat: {Op.forceflat} != {Opx.forceflat}"
raise ValueError(msg)
Op.forceflat = self.forceflat
else: # Only one of them is None
Op.forceflat = (
self.forceflat if self.forceflat is not None else Opx.forceflat
)
Op.dims = Opx.dims
return Op
elif np.isscalar(x):
Op = _ScaledLinearOperator(self, x)
self._copy_attributes(
Op,
exclude=["explicit", "name"],
)
Op.explicit = False
return Op
else:
if not get_ndarray_multiplication() and (
x.ndim > 2 or (x.ndim == 2 and x.shape[0] != self.shape[1])
):
msg = (
"Operator can only be applied to 1D vectors or 2D matrices. "
"Enable ndarray multiplication with pylops.set_ndarray_multiplication(True)."
)
raise ValueError(msg)
is_dims_shaped = x.shape == self.dims
is_dims_shaped_matrix = len(x.shape) > 1 and x.shape[:-1] == self.dims
if is_dims_shaped:
# (dims1, ..., dimsK) => (dims1 * ... * dimsK,) == self.shape
x = x.ravel()
if is_dims_shaped_matrix and not self.forceflat:
# (dims1, ..., dimsK, P) => (dims1 * ... * dimsK, P)
x = x.reshape((-1, x.shape[-1]))
if x.ndim == 1:
y = self.matvec(x)
if (
is_dims_shaped
and not self.forceflat
and get_ndarray_multiplication()
):
y = y.reshape(self.dimsd)
return y
elif x.ndim == 2:
y = self.matmat(x)
if (
is_dims_shaped_matrix
and not self.forceflat
and get_ndarray_multiplication()
):
y = y.reshape((*self.dimsd, -1))
return y
else:
msg = (
"Wrong shape.\nFor vector multiplication, expects either a 1d "
"array or, an ndarray of size `dims` when `dims` and `dimsd` "
"both are available.\n For matrix multiplication, expects a 2d "
f"array with its first dimension is equal to {self.shape[1]}.\n"
f"Instead, received an array of shape {x.shape}."
)
raise ValueError(msg)
[docs]
def div(
self,
y: float | NDArray,
niter: int = 100,
densesolver: str = "scipy",
) -> NDArray | LinearOperator:
r"""Solve the linear problem or divide by a scalar
Overloading of operator ``/`` to improve expressivity of `Pylops`
when solving inverse problems.
Division of a PyLops operator by a:
- scalar: returns a scaled operator, so ``Op / a`` is the same as
``Op * (1 / a)``.
- 1d-array: solves the linear problem
:math:`\mathbf{y}=\mathbf{A}\mathbf{x}` and returns the estimated
model. Explicit operators are solved with a sparse or dense direct
solver, whilst any other operator is solved iteratively;
Parameters
----------
y : :obj:`float` or :obj:`numpy.ndarray`
Scalar or data of size ``shape[0]``.
niter : :obj:`int`, optional
Number of iterations (to be used only when ``explicit=False``)
densesolver : :obj:`str`, optional
Use scipy (``scipy``) or numpy (``numpy``) dense solver (to be
used only for square explicit operators with dense matrix)
Returns
-------
xest : :obj:`pylops.LinearOperator` or :obj:`numpy.ndarray`
Scaled operator or estimated model.
"""
xest = self.__truediv__(y, niter=niter, densesolver=densesolver)
return xest
def __truediv__(
self,
y: NDArray | float,
niter: int = 100,
densesolver: str = "scipy",
) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator, symmetric with ``Op * a``
return self * (1 / y)
if self.explicit is True:
if sp.sparse.issparse(self.A):
# use scipy solver for sparse matrices
xest = spsolve(self.A, y)
elif isinstance(self.A, np.ndarray):
# use scipy solvers for dense matrices (used for backward
# compatibility, could be switched to numpy equivalents)
if self.A.shape[0] == self.A.shape[1]:
if densesolver == "scipy":
xest = sp_solve(self.A, y)
else:
xest = np_solve(self.A, y)
else:
xest = lstsq(self.A, y)[0]
else:
# use numpy/cupy solvers for dense matrices
ncp = get_array_module(y)
if self.A.shape[0] == self.A.shape[1]:
xest = ncp.linalg.solve(self.A, y)
else:
xest = ncp.linalg.lstsq(self.A, y)[0]
else:
if isinstance(y, np.ndarray):
# numpy backend
xest = lsqr(self, y, iter_lim=niter, atol=1e-8, btol=1e-8)[0]
else:
# cupy backend
ncp = get_array_module(y)
xest = cgls(
self,
y,
x0=ncp.zeros(int(self.shape[1]), dtype=self.dtype),
niter=niter,
)[0]
return xest
[docs]
def todense(
self,
backend: str = "numpy",
) -> NDArray:
r"""Return dense matrix.
The operator is converted into its dense matrix equivalent. In order
to do so, square or tall operators are applied to an identity matrix
whose number of rows and columns is equivalent to the number of
columns of the operator. Conversely, for skinny operators, the
transpose operator is applied to an identity matrix
whose number of rows and columns is equivalent to the number of
rows of the operator and the resulting matrix is transposed
(and complex conjugated).
Note that this operation may be costly for operators with large number
of rows and columns and it should be used mostly as a way to inspect
the structure of the matricial equivalent of the operator.
Parameters
----------
backend : :obj:`str`, optional
Backend used to densify matrix (``numpy`` or ``cupy`` or ``jax``). Note that
this must be consistent with how the operator has been created.
Returns
-------
matrix : :obj:`numpy.ndarray` or :obj:`cupy.ndarray`
Dense matrix.
"""
ncp = get_module(backend)
Op = self
# Create identity matrix
shapemin = min(Op.shape)
if shapemin <= 1e3:
# use numpy for small matrices (faster but heavier on memory)
identity = ncp.eye(shapemin, dtype=self.dtype)
else:
# use scipy for large matrices (slower but lighter on memory)
identity = get_sparse_eye(ncp.ones(1))(shapemin, dtype=self.dtype).tocsc()
# Apply operator
if Op.shape[1] == shapemin:
matrix = Op.matmat(identity)
else:
matrix = ncp.conj(Op.rmatmat(identity)).T
return matrix
[docs]
def tosparse(self) -> NDArray:
r"""Return sparse matrix.
The operator in converted into its sparse (CSR) matrix equivalent. In order
to do so, the operator is applied to series of unit vectors with length equal
to the number of coloumns in the original operator.
Returns
-------
matrix : :obj:`scipy.sparse.csr_matrix`
Sparse matrix.
"""
Op = self
(_, n) = self.shape
# stores non-zero data for the sparse matrix creation
entries = []
indices = []
# loop through columns of self
for i in range(n):
# make a unit vector for the ith column
unit_i = np.zeros(n)
unit_i[i] = 1
# multiply unit vector to self and find the non-zeros
res_i: NDArray = Op * unit_i
rows_nz = np.where(res_i != 0)[0]
# append the non-zero values and indices to the lists
for j in rows_nz:
indices.append([i, j])
entries_i = res_i[rows_nz]
for e in entries_i:
entries.append(e)
# post process the entries / indices for scipy.sparse.csr_matrix
entries = np.array(entries)
indices = np.array(indices)
i, j = indices[:, 0], indices[:, 1]
# construct a sparse, CSR matrix from the entries / indices data.
matrix = csr_matrix((entries, (j, i)), shape=self.shape, dtype=self.dtype)
return matrix
[docs]
def eigs(
self,
neigs: int | None = None,
symmetric: bool = False,
niter: int | None = None,
uselobpcg: bool = False,
backend: Tbackend = "numpy",
**kwargs_eig: dict[str, Any],
) -> NDArray:
r"""N-most significant eigenvalues of a linear operator.
Return an estimate of the ``neigs`` most significant eigenvalues
of the linear operator. If the operator has rectangular
shape (``shape[0] != shape[1]``), eigenvalues are first
computed for the square operator :math:`\mathbf{A^H}\mathbf{A}`
and the square-root values are returned.
.. note::
If ``backend="cupy"``, the eigenvalues of a square operator
can be computed only when in symmetric (``symmetric=True``),
since CuPy does not provide an equivalent
of :func:`scipy.sparse.linalg.eigs`. A
:obj:`NotImplementedError` is raised.
Parameters
----------
neigs : :obj:`int`
Number of eigenvalues to compute. Note that for ``explicit=False``,
only :math:`N-1` eigenvalues can be computed for real-valued
symmetric or non-square operators, and :math:`N-2` eigenvalues
otherwise, where :math:`N` is the size of the operator in the
model space. If ``None``, return :math:`N` for explicit case and
:math:`N-2` for implicit operators.
symmetric : :obj:`bool`, optional
Operator is real-symmetric or complex-hermitian (``True``)
or not (``False``). User should set this parameter to ``True``
only when it is guaranteed that the operator is real-symmetric
or complex-hermitian
niter : :obj:`int`, optional
Number of iterations for eigenvalue estimation
uselobpcg : :obj:`bool`, optional
Use :func:`scipy.sparse.linalg.lobpcg`. Note that this method
is only available for symmetric operators.
backend : :obj:`str`, optional
.. versionadded:: 2.9.0
Backend to use (`numpy` or `cupy`)
**kwargs_eig
Arbitrary keyword arguments for
:func:`scipy.linalg.eigvals`,
:func:`scipy.sparse.linalg.eigs`,
:func:`scipy.sparse.linalg.eigsh`, or
:func:`scipy.sparse.linalg.lobpcg` and
CuPy equivalents.
Returns
-------
eigenvalues : :obj:`numpy.ndarray`
Operator eigenvalues.
Raises
------
ValueError
If ``uselobpcg=True`` for a non-symmetric square matrix with
complex type, or if ``uselobpcg=True`` and the size of the operator
minus the number of constraints ``Y`` is smaller than 5 times ``neigs``
(in which case :func:`scipy.sparse.linalg.lobpcg` would switch to a
dense solver)
NotImplementedError
If ``backend="cupy"`` and ``symmetric=False`` for a square
operator for which only a limited number of eigenvalues is
requested (``neigs`` smaller than the size of the operator, or
``explicit=False``), as CuPy does not provide an equivalent of
:func:`scipy.sparse.linalg.eigs`
Notes
-----
Depending on the size of the operator, whether it is explicit or not
and the number of eigenvalues requested, different algorithms are
used by this routine.
More precisely, when only a limited number of eigenvalues is requested
the :func:`scipy.sparse.linalg.eigsh` method is used in case of
``symmetric=True`` and the :func:`scipy.sparse.linalg.eigs` method
is used ``symmetric=False``. However, when the matrix is represented
explicitly within the linear operator (``explicit=True``) and all the
eigenvalues are requested (or ``neigs`` is equal to the number of
columns of the operator minus one for a square, non-symmetric operator) the
:func:`scipy.linalg.eigvals` is used instead. Note that for
or a complex-valued operator, when ``neigs`` is equal to the number of
columns of the operator minus one, the eigenvalues are selected based on
``which`` and ``sigma`` is ignored.
Finally, when only a limited number of eigenvalues is required,
it is also possible to explicitly choose to use the
``scipy.sparse.linalg.lobpcg`` method via the ``uselobpcg`` input
parameter flag.
Most of these algorithms are a port of ARPACK [1]_, a Fortran package
which provides routines for quickly finding eigenvalues/eigenvectors
of a matrix. As ARPACK requires only left-multiplication
by the matrix in question, eigenvalues/eigenvectors can also be
estimated for linear operators when the dense matrix is not available.
.. [1] http://www.caam.rice.edu/software/ARPACK/
"""
def _check_lobpcg_size(n: int, neigs: int, Y: NDArray | None = None) -> None:
"""Check that scipy's lobpcg will iterate rather than use a dense solver
:func:`scipy.sparse.linalg.lobpcg` falls back to a dense eigensolver
(which cannot be applied to PyLops operators) when
``n - sizeY < 5 * neigs``, where ``sizeY`` is the number of
constraints (columns of ``Y``, which can be passed via ``kwargs_eig``).
"""
sizeY = Y.shape[1] if Y is not None and Y.ndim == 2 else 0
if n - sizeY < 5 * neigs:
msg = (
f"uselobpcg=True requires the operator size ({n}) minus the "
f"number of constraints ({sizeY}) to be at least 5 times the "
f"number of requested eigenvalues ({neigs}); reduce neigs to "
f"at most {(n - sizeY) // 5} or use uselobpcg=False"
)
raise ValueError(msg)
# Check if lobpcgc is feasible, raise error otherwise
square = self.shape[0] == self.shape[1]
if square and not symmetric and uselobpcg:
msg = (
"Cannot use scipy.sparse.linalg.lobpcg "
"for non real-symmetric or complex-hermitian operators..."
)
raise ValueError(msg)
ncp = get_module(backend)
explicit = self.explicit and isinstance(self.A, ncp.ndarray)
if neigs is None:
neigs = self.shape[1] if explicit else self.shape[1] - 2
# scipy.sparse.linalg.eigs (also used by eigsh for complex operators)
# requires neigs < N - 1, use dense solver instead
useeigs = (square and not symmetric) or np.iscomplexobj(
np.ones(1, dtype=self.dtype)
)
# Check if all eigenvalues should be computed
alleigs = explicit and (
neigs >= self.shape[1] or (useeigs and neigs == self.shape[1] - 1)
)
# Check that the number of eigenvalues requested is not larger
# than the number of columns in the operator minus two (eigs) or
# minus one (eigsh)
neigsmax = self.shape[1] - 2 if useeigs else self.shape[1] - 1
if not explicit and neigs > neigsmax:
warnings.warn(
f"The provided neigs ({neigs}) exceeds the number of columns "
f"in the operator minus {self.shape[1] - neigsmax}, "
f"reset to {neigsmax}",
stacklevel=2,
)
neigs = neigsmax
if square:
# Square operator
Op = self.A if explicit else self
if alleigs:
# All eigenvalues (only reached for explicit=True)
eigenvalues = get_eigvals(ncp.ones(1))(Op)
else:
# Limited number of eigenvalues
if symmetric and uselobpcg:
if not explicit:
_check_lobpcg_size(self.shape[0], neigs, kwargs_eig.get("Y"))
X = ncp.random.rand(self.shape[0], neigs).astype(self.dtype)
eigenvalues = get_lobpcg(ncp.ones(1))(
Op, X=X, maxiter=niter, **kwargs_eig
)[0]
elif symmetric:
eigenvalues = get_eigsh(ncp.ones(1))(
Op, k=neigs, maxiter=niter, **kwargs_eig
)[0]
else:
if backend == "cupy":
msg = (
"Eigenvalues of non-symmetric square "
"operators are not available with CuPy backend"
)
raise NotImplementedError(msg)
eigenvalues = sp_eigs(Op, k=neigs, maxiter=niter, **kwargs_eig)[0]
else:
# Non-square operator (square-root of eigenvalues of A^H A)
Op = ncp.dot(ncp.conj(self.A.T), self.A) if explicit else self.H * self
if alleigs:
# All eigenvalues (only reached for explicit=True)
eigenvalues = get_eigvals(ncp.ones(1))(Op)
elif uselobpcg:
# Limited number of eigenvalues
if not explicit:
_check_lobpcg_size(self.shape[1], neigs, kwargs_eig.get("Y"))
X = ncp.random.rand(self.shape[1], neigs).astype(self.dtype)
eigenvalues = get_lobpcg(ncp.ones(1))(
Op, X=X, maxiter=niter, **kwargs_eig
)[0]
else:
# Limited number of eigenvalues
eigenvalues = get_eigsh(ncp.ones(1))(
Op, k=neigs, maxiter=niter, **kwargs_eig
)[0]
eigenvalues = ncp.sqrt(eigenvalues)
if alleigs and neigs == self.shape[1] - 1:
# dense solver returns all eigenvalues: select neigs based on which
which = kwargs_eig.get("which", "LM")
key = {"M": ncp.abs, "R": ncp.real, "A": ncp.real, "I": ncp.imag}[which[1]]
key = key(eigenvalues)
idx = ncp.argsort(key if which[0] == "S" else -key)[:neigs]
eigenvalues = eigenvalues[idx]
return -ncp.sort(-eigenvalues)
[docs]
def cond(
self,
uselobpcg: bool = False,
backend: Tbackend = "numpy",
**kwargs_eig: dict[str, Any],
) -> NDArray:
r"""Condition number of linear operator.
Return an estimate of the condition number of the linear operator as
the ratio of the largest and lowest estimated eigenvalues.
.. note::
If ``backend="cupy"``, the condition number can only be computed
with the ``uselobpcg=True`` option (and therefore for square,
symmetric operators only).
Parameters
----------
uselobpcg : :obj:`bool`, optional
Use :func:`scipy.sparse.linalg.lobpcg` to compute eigenvalues
backend : :obj:`str`, optional
.. versionadded:: 2.9.0
Backend to use (`numpy` or `cupy`).
**kwargs_eig
Arbitrary keyword arguments for :func:`scipy.sparse.linalg.eigs`,
:func:`scipy.sparse.linalg.eigsh`, or
:func:`scipy.sparse.linalg.lobpcg`
Returns
-------
eigenvalues : :obj:`numpy.ndarray`
Operator eigenvalues.
Notes
-----
The condition number of a matrix (or linear operator) can be estimated
as the ratio of the largest and lowest estimated eigenvalues:
.. math::
k= \frac{\lambda_{max}}{\lambda_{min}}
The condition number provides an indication of the rate at which the
solution of the inversion of the linear operator :math:`A` will
change with respect to a change in the data :math:`y`.
Thus, if the condition number is large, even a small error in :math:`y`
may cause a large error in :math:`x`. On the other hand, if the
condition number is small then the error in :math:`x` is not much
bigger than the error in :math:`y`. A problem with a low condition
number is said to be *well-conditioned*, while a problem with a high
condition number is said to be *ill-conditioned*.
"""
if not uselobpcg:
lmax = self.eigs(neigs=1, which="LM", backend=backend, **kwargs_eig).item()
lmin = self.eigs(neigs=1, which="SM", backend=backend, **kwargs_eig).item()
cond = lmax / lmin
else:
lmax = self.eigs(
neigs=1,
symmetric=True,
uselobpcg=True,
largest=True,
backend=backend,
**kwargs_eig,
).item()
lmin = self.eigs(
neigs=1,
symmetric=True,
uselobpcg=True,
largest=False,
backend=backend,
**kwargs_eig,
).item()
cond = lmax / lmin
return cond
[docs]
def conj(self) -> LinearOperator:
"""Complex conjugate operator
Returns
-------
conjop : :obj:`pylops.LinearOperator`
Complex conjugate operator
"""
conjop = _ConjLinearOperator(self)
return conjop
[docs]
def apply_columns(self, cols: InputDimsLike) -> LinearOperator:
"""Apply subset of columns of operator
This method can be used to wrap a LinearOperator and mimic the action
of a subset of columns of the operator on a reduced model in forward
mode, and retrieve only the result of a subset of rows in adjoint mode.
Note that unless the operator has ``explicit=True``, this is not
optimal as the entire forward and adjoint passes of the original
operator will have to be perfomed. It can however be useful for the
implementation of solvers such as Orthogonal Matching Pursuit (OMP)
that iteratively build a solution by evaluate only a subset of the
columns of the operator.
Parameters
----------
cols : :obj:`list`
Columns to be selected
Returns
-------
colop : :obj:`pylops.LinearOperator`
Apply column operator
"""
colop = _ColumnLinearOperator(self, cols)
return colop
[docs]
def toreal(self, forw: bool = True, adj: bool = True) -> LinearOperator:
"""Real operator
Parameters
----------
forw : :obj:`bool`, optional
Apply real to output of forward pass
adj : :obj:`bool`, optional
Apply real to output of adjoint pass
Returns
-------
realop : :obj:`pylops.LinearOperator`
Real operator
"""
realop = _RealImagLinearOperator(self, forw, adj, True)
return realop
[docs]
def toimag(self, forw: bool = True, adj: bool = True) -> LinearOperator:
"""Imag operator
Parameters
----------
forw : :obj:`bool`, optional
Apply imag to output of forward pass
adj : :obj:`bool`, optional
Apply imag to output of adjoint pass
Returns
-------
imagop : :obj:`pylops.LinearOperator`
Imag operator
"""
imagop = _RealImagLinearOperator(self, forw, adj, False)
return imagop
[docs]
def trace(
self,
neval: int | None = None,
method: str | None = None,
backend: str = "numpy",
**kwargs_trace: dict[str, Any],
) -> float:
r"""Trace of linear operator.
Returns the trace (or its estimate) of the linear operator.
Parameters
----------
neval : :obj:`int`, optional
Maximum number of matrix-vector products compute. Default depends
``method``.
method : :obj:`str`, optional
Should be one of the following:
- **explicit**: If the operator is not explicit, will convert to
dense first.
- **hutchinson**: see :obj:`pylops.utils.trace_hutchinson`
- **hutch++**: see :obj:`pylops.utils.trace_hutchpp`
- **na-hutch++**: see :obj:`pylops.utils.trace_nahutchpp`
Defaults to 'explicit' for explicit operators, and 'Hutch++' for
the rest.
backend : :obj:`str`, optional
Backend used to densify matrix (``numpy`` or ``cupy``). Note that
this must be consistent with how the operator has been created.
**kwargs_trace
Arbitrary keyword arguments passed to
:obj:`pylops.utils.trace_hutchinson`,
:obj:`pylops.utils.trace_hutchpp`, or
:obj:`pylops.utils.trace_nahutchpp`
Returns
-------
trace : :obj:`float`
Operator trace.
Raises
------
ValueError
If the operator has rectangular shape (``shape[0] != shape[1]``)
NotImplementedError
If the ``method`` is not one of the available methods.
"""
if self.shape[0] != self.shape[1]:
msg = f"Operator must be square, got {self.shape}"
raise ValueError(msg)
ncp = get_module(backend)
if method is None:
method = "explicit" if self.explicit else "hutch++"
method_l = method.lower()
if method_l == "explicit":
A = self.A if self.explicit else self.todense(backend=backend)
return ncp.trace(A)
elif method_l == "hutchinson":
return trace_hutchinson(self, neval=neval, backend=backend, **kwargs_trace)
elif method_l == "hutch++":
return trace_hutchpp(self, neval=neval, backend=backend, **kwargs_trace)
elif method_l == "na-hutch++":
return trace_nahutchpp(self, neval=neval, backend=backend, **kwargs_trace)
else:
msg = f"`method={method}` not available."
raise NotImplementedError(msg)
[docs]
def reset_count(self) -> None:
"""Reset counters
When invoked all counters are set back to 0.
"""
self.matvec_count = 0
self.rmatvec_count = 0
self.matmat_count = 0
self.rmatmat_count = 0
class _ScaledLinearOperator(LinearOperator):
"""Scaled Linear Operator"""
def __init__(
self,
A: LinearOperator,
alpha: float,
) -> None:
if not isinstance(A, LinearOperator):
msg = "LinearOperator expected as A"
raise ValueError(msg)
if not np.isscalar(alpha):
msg = "Scalar expected as alpha"
raise ValueError(msg)
if isinstance(alpha, complex) and not np.iscomplexobj(
np.ones(1, dtype=A.dtype)
):
# if the scalar is of complex type but not the operator, find out type
dtype = _get_dtype([A], [type(alpha)])
else:
# if both the scalar and operator are of real or complex type, use type
# of the operator
dtype = A.dtype
super().__init__(dtype=dtype, shape=A.shape)
self.args = (A, alpha)
def _matvec(self, x: NDArray) -> NDArray:
return self.args[1] * self.args[0].matvec(x)
def _rmatvec(self, x: NDArray) -> NDArray:
return np.conj(self.args[1]) * self.args[0].rmatvec(x)
def _matmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return self.args[1] * self.args[0]._matmat(x, pool=pool)
def _rmatmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return np.conj(self.args[1]) * self.args[0]._rmatmat(x, pool=pool)
def _adjoint(self) -> LinearOperator:
A, alpha = self.args
return A.H * np.conj(alpha)
class _ConjLinearOperator(LinearOperator):
"""Complex conjugate linear operator"""
def __init__(self, Op: LinearOperator) -> None:
if not isinstance(Op, LinearOperator):
msg = "Op must be a LinearOperator"
raise TypeError(msg)
super().__init__(Op, shape=Op.shape)
self.Op = Op
def _matvec(self, x: NDArray) -> NDArray:
return (self.Op._matvec(x.conj())).conj()
def _rmatvec(self, x: NDArray) -> NDArray:
return (self.Op._rmatvec(x.conj())).conj()
def _adjoint(self) -> LinearOperator:
return _ConjLinearOperator(self.Op.H)
class _ColumnLinearOperator(LinearOperator):
"""Column selector linear operator
Produces the forward and adjoint passes with a subset of columns of the
original operator
"""
def __init__(
self,
Op: LinearOperator,
cols: InputDimsLike,
) -> None:
if not isinstance(Op, LinearOperator):
msg = "Op must be a LinearOperator"
raise TypeError(msg)
super().__init__(Op)
self.Op = Op
self.cols = cols
self._shape = (Op.shape[0], len(cols))
self._dims = len(cols)
if self.explicit:
self.Opcol = Op.A[:, cols]
def _matvec(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
if self.explicit:
y = self.Opcol @ x
else:
y = ncp.zeros(int(self.Op.shape[1]), dtype=self.dtype)
y[self.cols] = x
y = self.Op._matvec(y)
return y
def _rmatvec(self, x: NDArray) -> NDArray:
if self.explicit:
y = self.Opcol.T.conj() @ x
else:
y = self.Op._rmatvec(x)
y = y[self.cols]
return y
class _AdjointLinearOperator(LinearOperator):
"""Adjoint of Linear Operator"""
def __init__(self, A: LinearOperator):
shape = (A.shape[1], A.shape[0])
super().__init__(shape=shape, dtype=A.dtype)
self.A = A
self.args = (A,)
def _matvec(self, x: NDArray) -> NDArray:
return self.A._rmatvec(x)
def _rmatvec(self, x: NDArray) -> NDArray:
return self.A._matvec(x)
def _matmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return self.A._rmatmat(X, pool=pool)
def _rmatmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return self.A._matmat(X, pool=pool)
class _TransposedLinearOperator(LinearOperator):
"""Transposition of Linear Operator"""
def __init__(self, A: LinearOperator):
shape = (A.shape[1], A.shape[0])
super().__init__(shape=shape, dtype=A.dtype)
self.A = A
self.args = (A,)
def _matvec(self, x: NDArray) -> NDArray:
return np.conj(self.A._rmatvec(np.conj(x)))
def _rmatvec(self, x: NDArray) -> NDArray:
return np.conj(self.A._matvec(np.conj(x)))
def _matmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return np.conj(self.A._rmatmat(np.conj(X), pool=pool))
def _rmatmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return np.conj(self.A._matmat(np.conj(X), pool=pool))
class _ProductLinearOperator(LinearOperator):
"""Product of Linear Operators"""
def __init__(self, A: LinearOperator, B: LinearOperator):
if not isinstance(A, LinearOperator) or not isinstance(B, LinearOperator):
msg = (
f"Both operands have to be a LinearOperator - got {type(A)}, {type(B)}"
)
raise ValueError(msg)
if A.shape[1] != B.shape[0]:
msg = f"Cannot multiply {A} and {B}: shape mismatch"
raise ValueError(msg)
super().__init__(dtype=_get_dtype([A, B]), shape=(A.shape[0], B.shape[1]))
self.args = (A, B)
def _matvec(self, x: NDArray) -> NDArray:
return self.args[0]._matvec(self.args[1]._matvec(x))
def _rmatvec(self, x: NDArray) -> NDArray:
return self.args[1]._rmatvec(self.args[0]._rmatvec(x))
def _matmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return self.args[0]._matmat(self.args[1]._matmat(X), pool=pool)
def _rmatmat(self, X: NDArray, pool: Tpool | None = None) -> NDArray:
return self.args[1]._rmatmat(self.args[0]._rmatmat(X), pool=pool)
def _adjoint(self):
A, B = self.args
return B.H * A.H
class _SumLinearOperator(LinearOperator):
def __init__(
self,
A: LinearOperator,
B: LinearOperator,
) -> None:
if not isinstance(A, LinearOperator) or not isinstance(B, LinearOperator):
msg = "Both operands have to be a LinearOperator"
raise ValueError(msg)
if A.shape != B.shape:
msg = f"Cannot add {A} to {B}: shape mismatch"
raise ValueError(msg)
self.args = (A, B)
super().__init__(dtype=_get_dtype([A, B]), shape=A.shape)
def _matvec(self, x: NDArray) -> NDArray:
return self.args[0]._matvec(x) + self.args[1]._matvec(x)
def _rmatvec(self, x: NDArray) -> NDArray:
return self.args[0]._rmatvec(x) + self.args[1]._rmatvec(x)
def _matmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return self.args[0]._matmat(x, pool=pool) + self.args[1]._matmat(x, pool=pool)
def _rmatmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return self.args[0]._rmatmat(x, pool=pool) + self.args[1]._rmatmat(x, pool=pool)
def _adjoint(self) -> LinearOperator:
A, B = self.args
return A.H + B.H
class _PowerLinearOperator(LinearOperator):
def __init__(self, A: LinearOperator, p: int) -> None:
if not isinstance(A, LinearOperator):
msg = "LinearOperator expected as A"
raise ValueError(msg)
if A.shape[0] != A.shape[1]:
msg = f"Square LinearOperator expected, got {A}"
raise ValueError(msg)
if not isintlike(p) or p < 0:
msg = "Non-negative integer expected as p"
raise ValueError(msg)
super().__init__(dtype=A.dtype, shape=A.shape)
self.args = (A, p)
def _power(self, fun: Callable, x: NDArray, **kwargs) -> NDArray:
res = x.copy()
for _ in range(self.args[1]):
res = fun(res, **kwargs)
return res
def _matvec(self, x: NDArray) -> NDArray:
return self._power(self.args[0].matvec, x)
def _rmatvec(self, x: NDArray) -> NDArray:
return self._power(self.args[0].rmatvec, x)
def _matmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return self._power(self.args[0]._matmat, x, pool=pool)
def _rmatmat(self, x: NDArray, pool: Tpool | None = None) -> NDArray:
return self._power(self.args[0]._rmatmat, x, pool=pool)
def _adjoint(self) -> LinearOperator:
A, p = self.args
return A.H**p
class _RealImagLinearOperator(LinearOperator):
"""Real-Imag linear operator
Computes forward and adjoint passes of an operator Op and returns only
its real (or imaginary) component. Note that for the adjoint step the
output must be complex conjugated (i.e. opposite of the imaginary part is
returned)
"""
def __init__(
self,
Op: LinearOperator,
forw: bool = True,
adj: bool = True,
real: bool = True,
) -> None:
if not isinstance(Op, LinearOperator):
msg = "Op must be a LinearOperator"
raise TypeError(msg)
super().__init__(Op, shape=Op.shape)
self.Op = Op
self.real = real
self.forw = forw
self.adj = adj
self.dtype = np.array(0, dtype=self.Op.dtype).real.dtype
def _matvec(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
y = self.Op._matvec(x)
if self.forw:
if self.real:
y = ncp.real(y)
else:
y = ncp.imag(y)
return y
def _rmatvec(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
y = self.Op._rmatvec(x)
if self.adj:
if self.real:
y = ncp.real(y)
else:
y = -ncp.imag(y)
return y
def aslinearoperator(Op: spLinearOperator | LinearOperator) -> LinearOperator:
"""Return Op as a LinearOperator.
Converts any operator compatible with pylops definition of LinearOperator into a pylops
LinearOperator. This can be used for example when `Op` is a scipy operator to ensure
that the returned operator has all properties and methods of the pylops class.
Parameters
----------
Op : :obj:`pylops.LinearOperator` or any other Operator
Operator of any type
Returns
-------
Op : :obj:`pylops.LinearOperator`
Operator of type :obj:`pylops.LinearOperator`
"""
if isinstance(Op, LinearOperator):
return Op
else:
return LinearOperator(Op)
PyLops