matrix_representation

odl.matrix_representation(op)

Return a matrix representation of a linear operator.

Parameters

opOperator: The linear operator of which one wants a matrix representation. If the domain or range is a ProductSpace, it must be a power-space.

Returns

matrixnumpy.ndarray: The matrix representation of the operator. The shape will be op.domain.shape + op.range.shape and the dtype is the promoted (greatest) dtype of the domain and range.

Examples

Approximate a matrix on its own:

>>> mat = np.array([[1, 2, 3],
...                 [4, 5, 6],
...                 [7, 8, 9]])
>>> op = odl.MatrixOperator(mat)
>>> matrix_representation(op)
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

It also works with ProductSpace’s and higher dimensional TensorSpace’s. In this case, the returned “matrix” will also be higher dimensional:

>>> space = odl.uniform_discr([0, 0], [2, 2], (2, 2))
>>> grad = odl.Gradient(space)
>>> tensor = odl.matrix_representation(grad)
>>> tensor.shape == (2, 2, 2, 2, 2)
True

Since the “matrix” is now higher dimensional, we need to use e.g. numpy.tensordot if we want to compute with the matrix representation:

>>> x = space.element(lambda x: x[0] ** 2 + 2 * x[1] ** 2)
>>> grad(x)
ProductSpace(uniform_discr([ 0.,  0.], [ 2.,  2.], (2, 2)), 2).element([

        [[ 2.  ,  2.  ],
         [-2.75, -6.75]],

        [[ 4.  , -4.75],
         [ 4.  , -6.75]]
])
>>> np.tensordot(tensor, x.data, axes=grad.domain.ndim)
array([[[ 2.  ,  2.  ],
        [-2.75, -6.75]],

       [[ 4.  , -4.75],
        [ 4.  , -6.75]]])

Notes

The algorithm works by letting the operator act on all unit vectors, and stacking the output as a matrix.