MatrixOperator

class odl.MatrixOperator(*args, **kwargs)

Bases: Operator

A matrix acting as a linear operator.

This operator uses a matrix to represent an operator, and get its adjoint and inverse by doing computations on the matrix. This is in general a rather slow and memory-inefficient approach, and users are recommended to use other alternatives if possible.

__init__(matrix, domain=None, range=None, impl: str | None = None, device: str | None = None, axis=0)[source]

Initialize a new instance.

Parameters

matrixarray-like or scipy.sparse.base.spmatrix: 2-dimensional array representing the linear operator. For Scipy sparse matrices only tensor spaces with ndim == 1 are allowed as domain. The matrix is copied to impl/device, if these are specified (only once, when the operator is initialized). If a plain Python list is supplied, it will first be converted to a NumPy array.
domainTensorSpace, optional: Space of elements on which the operator can act. Its dtype must be castable to range.dtype. For the default None, a space with 1 axis and size matrix.shape[1] is used, together with the matrix’ data type.
rangeTensorSpace, optional: Space of elements on to which the operator maps. Its shape and dtype attributes must match those of the result of the multiplication. For the default None, the range is inferred from matrix, domain and axis.
implArrayBackend-identifying str, optional: Which backend to use for the low-level matrix multiplication. If not explicitly provided, it will be inferred in the following order of preference, depending on what is available: 1. from domain 2. from range 3. from matrix
devicestr, optional: On which device to store the matrix. Same defaulting logic as for impl.
axisint, optional: Sum over this axis of an input tensor in the multiplication.

Examples

By default, domain and range are spaces of with one axis:

>>> m = np.ones((3, 4))
>>> op = MatrixOperator(m)
>>> op.domain
rn(4)
>>> op.range
rn(3)
>>> op([1, 2, 3, 4])
rn(3).element([ 10.,  10.,  10.])

For multi-dimensional arrays (tensors), the summation (contraction) can be performed along a specific axis. In this example, the number of matrix rows (4) must match the domain shape entry in the given axis:

>>> dom = odl.rn((5, 4, 4))  # can use axis=1 or axis=2
>>> op = MatrixOperator(m, domain=dom, axis=1)
>>> op(dom.one()).shape
(5, 3, 4)
>>> op = MatrixOperator(m, domain=dom, axis=2)
>>> op(dom.one()).shape
(5, 4, 3)

The operator also works on uniform_discr type spaces. Note, however, that the weighting of the domain is propagated to the range by default, in order to keep the correspondence between adjoint and transposed matrix:

>>> space = odl.uniform_discr(0, 1, 4)
>>> op = MatrixOperator(m, domain=space)
>>> op(space.one())
rn(3, weighting=0.25).element([ 4.,  4.,  4.])
>>> np.array_equal(op.adjoint.matrix, m.T)
True

Notes

For a matrix $A \in \mathbb{F}^{n \times m}$ , the operation on a tensor $T \in \mathbb{F}^{n_1 \times \dots \times n_d}$ is defined as the summation

$(A \cdot T)_{i_1, \dots, i_k, \dots, i_d} = \sum_{j=1}^m A_{i_k j} T_{i_1, \dots, j, \dots, i_d}.$

It produces a new tensor $A \cdot T \in \mathbb{F}^{ n_1 \times \dots \times n \times \dots \times n_d}$ .