matrix_representation¶

odl.operator.oputils.matrix_representation(op)[source]¶

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.

Notes

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

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, axes=grad.domain.ndim)
array([[[ 2.  ,  2.  ],
        [-2.75, -6.75]],

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