scico.linop

Linear operator functions and classes.

Functions

jacobian(F, u[, include_eval])	Construct Jacobian linear operator for a general operator.
linop_from_function(f, classname[, f_name])	Make a LinearOperator from a function.
operator_norm(A[, maxiter, key])	Estimate the norm of a LinearOperator.
power_iteration(A[, maxiter, key])	Compute largest eigenvalue of a diagonalizable LinearOperator.
valid_adjoint(A, AT[, eps, x, y, key])	Check whether LinearOperator AT is the adjoint of A.

Classes

CircularConvolve(h, input_shape[, ndims, ...])	A circular convolution linear operator.
ComposedLinearOperator(A, B[, jit])	A composition of two LinearOperator objects.
Convolve(h, input_shape[, input_dtype, ...])	A convolution linear operator.
Crop(crop_width, input_shape[, input_dtype, jit])	A linear operator for cropping an array.
DFT(input_shape[, axes, axes_shape, norm, jit])	Multi-dimensional discrete Fourier transform.
Diagonal(diagonal[, input_shape, input_dtype])	Diagonal linear operator.
DiagonalStack(ops[, collapse_input, ...])	A diagonal stack of linear operators.
FiniteDifference(input_shape[, input_dtype, ...])	Finite difference operator.
Identity(input_shape[, input_dtype])	Identity operator.
LinearOperator(input_shape[, output_shape, ...])	Generic linear operator base class
MatrixOperator(A[, input_cols])	Linear operator implementing matrix multiplication.
Pad(input_shape, *args[, input_dtype, ...])	Linear operator version of scico.numpy.pad.
Parallel2dProjector(im_shape, angles[, x0, ...])	Parallel ray, single axis, 2D X-ray projector.
Reshape(input_shape, *args[, input_dtype, ...])	Linear operator version of scico.numpy.reshape.
SingleAxisFiniteDifference(input_shape[, ...])	Finite difference operator acting along a single axis.
Slice(idx, input_shape[, input_dtype, jit])	A linear operator for slicing an array.
Sum(input_shape, *args[, input_dtype, ...])	Linear operator version of scico.numpy.sum.
Transpose(input_shape, *args[, input_dtype, ...])	Linear operator version of scico.numpy.transpose.
VerticalStack(ops[, collapse_output, jit])	A vertical stack of linear operators.
XRayTransform(projector)	X-ray transform operator.

class scico.linop.CircularConvolve(h, input_shape, ndims=None, input_dtype=<class 'jax.numpy.float32'>, h_is_dft=False, h_center=None, jit=True, **kwargs)

Bases: LinearOperator

A circular convolution linear operator.

This linear operator implements circular, multi-dimensional convolution via pointwise multiplication in the DFT domain. In its simplest form, it implements a single convolution and can be represented by linear operator \(H\) such that

\[H \mb{x} = \mb{h} \ast \mb{x} \;,\]

where \(\mb{h}\) is a user-defined filter.

More complex forms, corresponding to the case where either the input (as represented by parameter input_shape) or filter (parameter h) have additional axes that are not involved in the convolution are also supported. These follow numpy broadcasting rules. For example:

Additional axes in the input \(\mb{x}\) and not in \(\mb{h}\) corresponds to the operation

\[\begin{split}H \mb{x} = \left( \begin{array}{ccc} H' & 0 & \ldots\\ 0 & H' & \ldots\\ \vdots & \vdots & \ddots \end{array} \right) \left( \begin{array}{c} \mb{x}_0\\ \mb{x}_1\\ \vdots \end{array} \right) \;.\end{split}\]

Additional axes in \(\mb{h}\) corresponds to multiple filters, which will be denoted by \(\{\mb{h}_m\}\), with corresponding individual linear operations being denoted by \(h_m \mb{x}_m = \mb{h}_m \ast \mb{x}_m\). The full linear operator can then be represented as

\[\begin{split}H \mb{x} = \left( \begin{array}{c} H_0\\ H_1\\ \vdots \end{array} \right) \mb{x} \;.\end{split}\]

if the input is singleton, and as

\[\begin{split}H \mb{x} = \left( \begin{array}{ccc} H_0 & 0 & \ldots\\ 0 & H_1 & \ldots\\ \vdots & \vdots & \ddots \end{array} \right) \left( \begin{array}{c} \mb{x}_0\\ \mb{x}_1\\ \vdots \end{array} \right)\end{split}\]

otherwise.

Parameters:

• h (Array) – Array of filters.

• input_shape (Tuple[int, ...]) – Shape of input array.

• ndims (Optional[int]) – Number of (trailing) dimensions of the input and h involved in the convolution. Defaults to the number of dimensions in the input.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• h_is_dft (bool) – Flag indicating whether h is in the DFT domain.

• h_center (Union[Array, Sequence, float, int, None]) – Array of length ndims specifying the center of the filter. Defaults to the upper left corner, i.e., h_center = [0, 0, …, 0], may be noninteger. May be a float or int if h is one-dimensional.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

static from_operator(H, ndims=None, center=None, jit=True)

Construct a CircularConvolve version of a given operator.

Construct a CircularConvolve version of a given operator, which is assumed to be linear and shift invariant (LSI).

Parameters:

• H (Operator) – Input operator.

• ndims (Optional[int]) – Number of trailing dims over which the H acts.

• center (Optional[Tuple[int, ...]]) – Location at which to place the Kronecker delta. For LSI inputs, this will not matter. Defaults to the center of H.input_shape, i.e., (n_1 // 2, n_2 // 2, …).

• jit (bool) – If True, jit the resulting CircularConvolve.

class scico.linop.Convolve(h, input_shape, input_dtype=<class 'numpy.float32'>, mode='full', jit=True, **kwargs)

Bases: LinearOperator

A convolution linear operator.

Wrap jax.scipy.signal.convolve as a LinearOperator.

Wrap jax.scipy.signal.convolve as a LinearOperator.

Parameters:

• h (Array) – Convolutional filter. Must have same number of dimensions as len(input_shape).

• input_shape (Tuple[int, ...]) – Shape of input array.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• mode (str) – A string indicating the size of the output. One of “full”, “valid”, “same”. Defaults to “full”.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

For more details on mode, see jax.scipy.signal.convolve.

class scico.linop.DFT(input_shape, axes=None, axes_shape=None, norm=None, jit=True, **kwargs)

Bases: LinearOperator

Multi-dimensional discrete Fourier transform.

Parameters:

• input_shape (Tuple[int, ...]) – Shape of input array.

• axes (Optional[Sequence]) – Axes over which to compute the DFT. If None, the DFT is computed over all axes.

• axes_shape (Optional[Tuple[int, ...]]) – Output shape on the subset of array axes selected by axes. This parameter has the same behavior as the s parameter of numpy.fft.fftn.

• norm (Optional[str]) – DFT normalization mode. See the norm parameter of numpy.fft.fftn.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

inv(z)

Compute the inverse of this LinearOperator.

Compute the inverse of this LinearOperator applied to z.

Parameters:

z (Array) – Input array to inverse DFT.

Return type:

Array

class scico.linop.Diagonal(diagonal, input_shape=None, input_dtype=None, **kwargs)

Bases: LinearOperator

Diagonal linear operator.

Parameters:

• diagonal (Union[Array, BlockArray]) – Diagonal elements of this LinearOperator.

• input_shape (Optional[Tuple[int, ...]]) – Shape of input array. By default, equal to diagonal.shape, but may also be set to a shape that is broadcast-compatible with diagonal.shape.

• input_dtype (Optional[DType]) – dtype of input argument. The default, None, means diagonal.dtype.

property H: Diagonal

Hermitian transpose of this Diagonal.

property T: Diagonal

Transpose of this Diagonal.

conj()

Complex conjugate of this Diagonal.

Return type:

Diagonal

property gram_op: Diagonal

Gram operator of this Diagonal.

Return a new Diagonal G such that G(x) = A.adj(A(x))).

norm(ord=None)

Compute the matrix norm of the diagonal operator.

Valid values of ord and the corresponding norm definition are those listed under “norm for matrices” in the scico.numpy.linalg.norm documentation.

class scico.linop.Identity(input_shape, input_dtype=<class 'jax.numpy.float32'>, **kwargs)

Bases: Diagonal

Identity operator.

Parameters:

input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

class scico.linop.FiniteDifference(input_shape, input_dtype=<class 'numpy.float32'>, axes=None, prepend=None, append=None, circular=False, jit=True, **kwargs)

Bases: VerticalStack

Finite difference operator.

Computes finite differences along the specified axes, returning the results in a jax.Array (when possible) or BlockArray. See VerticalStack for details on how this choice is made. See SingleAxisFiniteDifference for the mathematical implications of the different boundary handling options prepend, append, and circular.

Example

>>> A = FiniteDifference((2, 3))
>>> x = snp.array([[1, 2, 4],
...                [0, 4, 1]])
>>> (A @ x)[0]
Array([[-1,  2, -3]], dtype=int32)
>>> (A @ x)[1]
Array([[ 1,  2],
       [ 4, -3]], dtype=int32)

Parameters:

• input_shape (Tuple[int, ...]) – Shape of input array.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• axes (Union[Tuple[int, ...], int, None]) – Axis or axes over which to apply finite difference operator. If not specified, or None, differences are evaluated along all axes.

• prepend (Union[Literal[0], Literal[1], None]) – Flag indicating handling of the left/top/etc. boundary. If None, there is no boundary extension. Values of 0 or 1 indicate respectively that zeros or the initial value in the array are prepended to the difference array.

• append (Union[Literal[0], Literal[1], None]) – Flag indicating handling of the right/bottom/etc. boundary. If None, there is no boundary extension. Values of 0 or 1 indicate respectively that zeros or -1 times the final value in the array are appended to the difference array.

• circular (bool) – If True, perform circular differences, i.e., include x[-1] - x[0]. If True, prepend and append must both be None.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

class scico.linop.SingleAxisFiniteDifference(input_shape, input_dtype=<class 'numpy.float32'>, axis=-1, prepend=None, append=None, circular=False, jit=True, **kwargs)

Bases: LinearOperator

Finite difference operator acting along a single axis.

By default (i.e. prepend and append set to None and circular set to False), the difference operator corresponds to the matrix

\[\begin{split}\left(\begin{array}{rrrrr} -1 & 1 & 0 & \ldots & 0\\ 0 & -1 & 1 & \ldots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \ldots & -1 & 1 \end{array}\right) \;,\end{split}\]

mapping \(\mbb{R}^N \rightarrow \mbb{R}^{N-1}\), while if circular is True, it corresponds to the \(\mbb{R}^N \rightarrow \mbb{R}^N\) mapping

\[\begin{split}\left(\begin{array}{rrrrr} -1 & 1 & 0 & \ldots & 0\\ 0 & -1 & 1 & \ldots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \ldots & -1 & 1\\ 1 & 0 & \dots & 0 & -1 \end{array}\right) \;.\end{split}\]

Other possible choices include prepend set to None and append set to 0, giving the \(\mbb{R}^N \rightarrow \mbb{R}^N\) mapping

\[\begin{split}\left(\begin{array}{rrrrr} -1 & 1 & 0 & \ldots & 0\\ 0 & -1 & 1 & \ldots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \ldots & -1 & 1\\ 0 & 0 & \dots & 0 & 0 \end{array}\right) \;,\end{split}\]

and both prepend and append set to 1, giving the \(\mbb{R}^N \rightarrow \mbb{R}^{N+1}\) mapping

\[\begin{split}\left(\begin{array}{rrrrr} 1 & 0 & 0 & \ldots & 0\\ -1 & 1 & 0 & \ldots & 0\\ 0 & -1 & 1 & \ldots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \ldots & -1 & 1\\ 0 & 0 & \dots & 0 & -1 \end{array}\right) \;.\end{split}\]

Parameters:

• input_shape (Tuple[int, ...]) – Shape of input array.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• axis (int) – Axis over which to apply finite difference operator.

• prepend (Union[Literal[0], Literal[1], None]) – Flag indicating handling of the left/top/etc. boundary. If None, there is no boundary extension. Values of 0 or 1 indicate respectively that zeros or the initial value in the array are prepended to the difference array.

• append (Union[Literal[0], Literal[1], None]) – Flag indicating handling of the right/bottom/etc. boundary. If None, there is no boundary extension. Values of 0 or 1 indicate respectively that zeros or -1 times the final value in the array are appended to the difference array.

• circular (bool) – If True, perform circular differences, i.e., include x[-1] - x[0]. If True, prepend and append must both be None.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

class scico.linop.Crop(crop_width, input_shape, input_dtype=<class 'jax.numpy.float32'>, jit=True, **kwargs)

Bases: LinearOperator

A linear operator for cropping an array.

Parameters:

• crop_width (Union[int, Sequence]) – Specify the crop width using the same format as the pad_width parameter of snp.pad.

• input_shape (Tuple[int, ...]) – Shape of input jax.Array.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

class scico.linop.Pad(input_shape, *args, input_dtype=<class 'jax.numpy.float32'>, output_shape=None, output_dtype=None, jit=True, **kwargs)

Bases: LinearOperator

Linear operator version of scico.numpy.pad.

Parameters:

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

• args (Any) – Positional arguments passed to scico.numpy.pad.

• input_dtype (DType) – dtype for input argument. Defaults to float32. If the LinearOperator implements complex-valued operations, this must be a complex dtype (typically complex64) for correct adjoint and gradient calculation.

• output_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...], None]) – Shape of output array. Defaults to None. If None, output_shape is determined by evaluating self.__call__ on an input array of zeros.

• output_dtype (Optional[DType]) – dtype for output argument. Defaults to None. If None, output_dtype is determined by evaluating self.__call__ on an input array of zeros.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

• kwargs (Any) – Keyword arguments passed to scico.numpy.pad.

class scico.linop.Reshape(input_shape, *args, input_dtype=<class 'jax.numpy.float32'>, output_shape=None, output_dtype=None, jit=True, **kwargs)

Bases: LinearOperator

Linear operator version of scico.numpy.reshape.

Parameters:

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

• args (Any) – Positional arguments passed to scico.numpy.reshape.

• input_dtype (DType) – dtype for input argument. Defaults to float32. If the LinearOperator implements complex-valued operations, this must be a complex dtype (typically complex64) for correct adjoint and gradient calculation.

• output_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...], None]) – Shape of output array. Defaults to None. If None, output_shape is determined by evaluating self.__call__ on an input array of zeros.

• output_dtype (Optional[DType]) – dtype for output argument. Defaults to None. If None, output_dtype is determined by evaluating self.__call__ on an input array of zeros.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

• kwargs (Any) – Keyword arguments passed to scico.numpy.reshape.

class scico.linop.Slice(idx, input_shape, input_dtype=<class 'jax.numpy.float32'>, jit=True, **kwargs)

Bases: LinearOperator

A linear operator for slicing an array.

This operator may be applied to either a jax.Array or a BlockArray. In the latter case, parameter idx must conform to the BlockArray indexing requirements.

Parameters:

• idx (ArrayIndex) – An array indexing expression, as generated by numpy.s_, for example.

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input jax.Array or BlockArray.

• input_dtype (DType) – dtype for input argument. Defaults to float32.

• jit (bool) – If True, jit the evaluation, adjoint, and gram functions of the LinearOperator.

class scico.linop.Sum(input_shape, *args, input_dtype=<class 'jax.numpy.float32'>, output_shape=None, output_dtype=None, jit=True, **kwargs)

Bases: LinearOperator

Linear operator version of scico.numpy.sum.

Parameters:

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

• args (Any) – Positional arguments passed to scico.numpy.sum.

• input_dtype (DType) – dtype for input argument. Defaults to float32. If the LinearOperator implements complex-valued operations, this must be a complex dtype (typically complex64) for correct adjoint and gradient calculation.

• output_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...], None]) – Shape of output array. Defaults to None. If None, output_shape is determined by evaluating self.__call__ on an input array of zeros.

• output_dtype (Optional[DType]) – dtype for output argument. Defaults to None. If None, output_dtype is determined by evaluating self.__call__ on an input array of zeros.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

• kwargs (Any) – Keyword arguments passed to scico.numpy.sum.

class scico.linop.Transpose(input_shape, *args, input_dtype=<class 'jax.numpy.float32'>, output_shape=None, output_dtype=None, jit=True, **kwargs)

Bases: LinearOperator

Linear operator version of scico.numpy.transpose.

Parameters:

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

• args (Any) – Positional arguments passed to scico.numpy.transpose.

• input_dtype (DType) – dtype for input argument. Defaults to float32. If the LinearOperator implements complex-valued operations, this must be a complex dtype (typically complex64) for correct adjoint and gradient calculation.

• output_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...], None]) – Shape of output array. Defaults to None. If None, output_shape is determined by evaluating self.__call__ on an input array of zeros.

• output_dtype (Optional[DType]) – dtype for output argument. Defaults to None. If None, output_dtype is determined by evaluating self.__call__ on an input array of zeros.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

• kwargs (Any) – Keyword arguments passed to scico.numpy.transpose.

scico.linop.linop_from_function(f, classname, f_name=None)

Make a LinearOperator from a function.

Example

>>> Sum = linop_from_function(snp.sum, 'Sum')
>>> H = Sum((2, 10), axis=1)
>>> H @ snp.ones((2, 10))
Array([10., 10.], dtype=float32)

Parameters:

• f (Callable) – Function from which to create a LinearOperator.

• classname (str) – Name of the resulting class.

• f_name (Optional[str]) – Name of f for use in docstrings. Useful for getting the correct version of wrapped functions. Defaults to f”{f.__module__}.{f.__name__}”.

class scico.linop.ComposedLinearOperator(A, B, jit=False)

Bases: LinearOperator

A composition of two LinearOperator objects.

A new LinearOperator formed by the composition of two other LinearOperator objects.

A ComposedLinearOperator AB implements AB @ x == A @ B @ x. LinearOperator A and B are stored as attributes of the ComposedLinearOperator.

LinearOperator A and B must have compatible shapes and dtypes: A.input_shape == B.output_shape and A.input_dtype == B.input_dtype.

Parameters:

• A (LinearOperator) – First (left) LinearOperator.

• B (LinearOperator) – Second (right) LinearOperator.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

class scico.linop.LinearOperator(input_shape, output_shape=None, eval_fn=None, adj_fn=None, input_dtype=<class 'numpy.float32'>, output_dtype=None, jit=False)

Bases: Operator

Generic linear operator base class

Parameters:

• input_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...]]) – Shape of input array.

• output_shape (Union[Tuple[int, ...], Tuple[Tuple[int, ...], ...], None]) – Shape of output array. Defaults to None. If None, output_shape is determined by evaluating self.__call__ on an input array of zeros.

• eval_fn (Optional[Callable]) – Function used in evaluating this LinearOperator. Defaults to None. If None, then self.__call__ must be defined in any derived classes.

• adj_fn (Optional[Callable]) – Function used to evaluate the adjoint of this LinearOperator. Defaults to None. If None, the adjoint is not set, and the _set_adjoint will be called silently at the first adj call or can be called manually.

• input_dtype (DType) – dtype for input argument. Defaults to float32. If the LinearOperator implements complex-valued operations, this must be a complex dtype (typically complex64) for correct adjoint and gradient calculation.

• output_dtype (Optional[DType]) – dtype for output argument. Defaults to None. If None, output_dtype is determined by evaluating self.__call__ on an input array of zeros.

• jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

property H: LinearOperator

Hermitian transpose of this LinearOperator.

Return a new LinearOperator that is the Hermitian transpose of this LinearOperator. For a real-valued LinearOperator A (A.input_dtype is float32 or float64), the LinearOperator A.H is equivalent to A.T. For a complex-valued LinearOperator A (A.input_dtype is complex64 or complex128), the LinearOperator A.H implements the adjoint of A : A.H @ y == A.adj(y) == A.conj().T @ y).

For the non-conjugate transpose, see T.

property T: LinearOperator

Transpose of this LinearOperator.

Return a new LinearOperator that implements the transpose of this LinearOperator. For a real-valued LinearOperator A (A.input_dtype is float32 or float64), the LinearOperator A.T implements the adjoint: A.T(y) == A.adj(y). For a complex-valued LinearOperator A (A.input_dtype is complex64 or complex128), the LinearOperator A.T is not the adjoint. For the conjugate transpose, use .conj().T or H.

__call__(x)

Evaluate this LinearOperator at the point \(\mb{x}\).

Parameters:

x (Union[LinearOperator, Array, BlockArray]) – Point at which to evaluate this LinearOperator. If x is a jax.Array or BlockArray, must have shape == self.input_shape. If x is a LinearOperator, must have x.output_shape == self.input_shape.

Return type:

Union[LinearOperator, Array, BlockArray]

adj(y)

Adjoint of this LinearOperator.

Compute the adjoint of this LinearOperator applied to input y.

Parameters:

y (Union[LinearOperator, Array, BlockArray]) – Point at which to compute adjoint. If y is jax.Array or BlockArray, must have shape == self.output_shape. If y is a LinearOperator, must have y.output_shape == self.output_shape.

Return type:

Union[LinearOperator, Array, BlockArray]

Returns:

Adjoint evaluated at y.

conj()

Complex conjugate of this LinearOperator.

Return a new LinearOperator Ac such that Ac(x) = conj(A)(x).

Return type:

LinearOperator

gram(x)

Compute A.adj(A(x)).

Parameters:

x (Union[LinearOperator, Array, BlockArray]) – Point at which to evaluate the gram operator. If x is a jax.Array or BlockArray, must have shape == self.input_shape. If x is a LinearOperator, must have x.output_shape == self.input_shape.

Return type:

Union[LinearOperator, Array, BlockArray]

Returns:

Result of A.adj(A(x)).

property gram_op: LinearOperator

Gram operator of this LinearOperator.

Return a new LinearOperator G such that G(x) = A.adj(A(x))).

jit()

Replace the private functions _eval, _adj, _gram with jitted versions.

class scico.linop.MatrixOperator(A, input_cols=0)

Bases: LinearOperator

Linear operator implementing matrix multiplication.

Parameters:

• A (Union[Array, ndarray, bool_, number, bool, int, float, complex]) – Dense array. The action of the created LinearOperator will implement matrix multiplication with A.

• input_cols (int) – If this parameter is set to the default of 0, the MatrixOperator takes a vector (one-dimensional array) input. If the input is intended to be a matrix (two-dimensional array), this parameter should specify number of columns in the matrix.

property H

Hermitian (conjugate) transpose of this MatrixOperator.

Return a MatrixOperator corresponding to the Hermitian (conjugate) transpose of this matrix.

property T

Transpose of this MatrixOperator.

Return a MatrixOperator corresponding to the transpose of this matrix.

__call__(other)

Evaluate this LinearOperator at the point \(\mb{x}\).

Parameters:

x – Point at which to evaluate this LinearOperator. If x is a jax.Array or BlockArray, must have shape == self.input_shape. If x is a LinearOperator, must have x.output_shape == self.input_shape.

adj(y)

Adjoint of this LinearOperator.

Compute the adjoint of this LinearOperator applied to input y.

Parameters:

y – Point at which to compute adjoint. If y is jax.Array or BlockArray, must have shape == self.output_shape. If y is a LinearOperator, must have y.output_shape == self.output_shape.

Returns:

Adjoint evaluated at y.

conj()

Complex conjugate of this MatrixOperator.

Return a MatrixOperator with complex conjugated elements.

gram(other)

Compute A.adj(A(x)).

Parameters:

x – Point at which to evaluate the gram operator. If x is a jax.Array or BlockArray, must have shape == self.input_shape. If x is a LinearOperator, must have x.output_shape == self.input_shape.

Returns:

Result of A.adj(A(x)).

property gram_op

Gram operator of this MatrixOperator.

Return a new LinearOperator G such that G(x) = A.adj(A(x))).

norm(ord=None, axis=None, keepdims=False)

Compute the norm of the dense matrix self.A.

Call scico.numpy.linalg.norm on the dense matrix self.A.

to_array()

Return a numpy.ndarray containing self.A.

class scico.linop.DiagonalStack(ops, collapse_input=True, collapse_output=True, jit=True, **kwargs)

Bases: DiagonalStack, LinearOperator

A diagonal stack of linear operators.

Given linear operators \(A_1, A_2, \dots, A_N\), create the linear operator

\[\begin{split}H = \begin{pmatrix} A_1 & 0 & \ldots & 0\\ 0 & A_2 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & A_N \\ \end{pmatrix} \qquad \text{such that} \qquad H \begin{pmatrix} \mb{x}_1 \\ \mb{x}_2 \\ \vdots \\ \mb{x}_N \\ \end{pmatrix} = \begin{pmatrix} A_1(\mb{x}_1) \\ A_2(\mb{x}_2) \\ \vdots \\ A_N(\mb{x}_N) \\ \end{pmatrix} \;.\end{split}\]

By default, if the inputs \(\mb{x}_1, \mb{x}_2, \dots, \mb{x}_N\) all have the same (possibly nested) shape, S, this operator will work on the stack, i.e., have an input shape of (N, *S). If the inputs have distinct shapes, S1, S2, …, SN, this operator will work on the block concatenation, i.e., have an input shape of (S1, S2, …, SN). The same holds for the output shape.

Parameters:

• ops (Sequence[LinearOperator]) – Operators to stack.

• collapse_input (Optional[bool]) – If True, inputs are expected to be stacked along the first dimension when possible.

• collapse_output (Optional[bool]) – If True, the output will be stacked along the first dimension when possible.

• jit (bool) – See jit in LinearOperator.

class scico.linop.VerticalStack(ops, collapse_output=True, jit=True, **kwargs)

Bases: VerticalStack, LinearOperator

A vertical stack of linear operators.

Given linear operators \(A_1, A_2, \dots, A_N\), create the linear operator

\[\begin{split}H = \begin{pmatrix} A_1 \\ A_2 \\ \vdots \\ A_N \\ \end{pmatrix} \qquad \text{such that} \qquad H \mb{x} = \begin{pmatrix} A_1(\mb{x}) \\ A_2(\mb{x}) \\ \vdots \\ A_N(\mb{x}) \\ \end{pmatrix} \;.\end{split}\]

Parameters:

• ops (Sequence[LinearOperator]) – Linear operators to stack.

• collapse_output (Optional[bool]) – If True and the output would be a BlockArray with shape ((m, n, …), (m, n, …), …), the output is instead a jax.Array with shape (S, m, n, …) where S is the length of ops.

• jit (bool) – See jit in LinearOperator.

scico.linop.jacobian(F, u, include_eval=False)

Construct Jacobian linear operator for a general operator.

For a specified Operator, construct a corresponding Jacobian LinearOperator, the application of which is equivalent to multiplication by the Jacobian of the Operator at a specified input value.

The implementation of this function is based on Operator.jvp and Operator.vjp, which are themselves based on jax.jvp and jax.vjp. For reasons of computational efficiency, these functions return the value of the Operator evaluated at the specified point in addition to the requested Jacobian-vector product. If the include_eval parameter of this function is True, the constructed LinearOperator returns a BlockArray output, the first component of which is the result of the Operator evaluation, and the second component of which is the requested Jacobian-vector product. If include_eval is False, then the Operator evaluation computed by jax.jvp and jax.vjp are discarded.

Parameters:

• F (Operator) – Operator of which the Jacobian is to be computed.

• u (Array) – Input value of the Operator at which the Jacobian is to be computed.

• include_eval (Optional[bool]) – Flag indicating whether the result of evaluating the Operator should be included (as the first component of a BlockArray) in the output of the Jacobian LinearOperator constructed by this function.

Return type:

LinearOperator

Returns:

A LinearOperator capable of computing Jacobian-vector products.

scico.linop.operator_norm(A, maxiter=100, key=None)

Estimate the norm of a LinearOperator.

Estimate the operator norm induced by the \(\ell_2\) vector norm, i.e. for LinearOperator \(A\),

\[\begin{split}\| A \|_2 &= \max \{ \| A \mb{x} \|_2 \, : \, \| \mb{x} \|_2 \leq 1 \} \\ &= \sqrt{ \lambda_{ \mathrm{max} }( A^H A ) } = \sigma_{\mathrm{max}}(A) \;,\end{split}\]

where \(\lambda_{\mathrm{max}}(B)\) and \(\sigma_{\mathrm{max}}(B)\) respectively denote the largest eigenvalue of \(B\) and the largest singular value of \(B\). The value is estimated via power iteration, using power_iteration, to estimate \(\lambda_{\mathrm{max}}(A^H A)\).

Parameters:

• A (LinearOperator) – LinearOperator for which operator norm is desired.

• maxiter (int) – Maximum number of power iterations to use. Default: 100

• key (Optional[Array]) – Jax PRNG key. Defaults to None, in which case a new key is created.

Returns:

float – Norm of operator \(A\).

scico.linop.power_iteration(A, maxiter=100, key=None)

Compute largest eigenvalue of a diagonalizable LinearOperator.

Compute largest eigenvalue of a diagonalizable LinearOperator using power iteration.

Parameters:

• A (LinearOperator) – LinearOperator used for computation. Must be diagonalizable.

• maxiter (int) – Maximum number of power iterations to use.

• key (Optional[Array]) – Jax PRNG key. Defaults to None, in which case a new key is created.

Returns:

tuple –

A tuple (mu, v) containing:

• mu: Estimate of largest eigenvalue of A.

• v: Eigenvector of A with eigenvalue mu.

scico.linop.valid_adjoint(A, AT, eps=1e-07, x=None, y=None, key=None)

Check whether LinearOperator AT is the adjoint of A.

Check whether LinearOperator \(\mathsf{AT}\) is the adjoint of \(\mathsf{A}\). The test exploits the identity

\[\mathbf{y}^T (A \mathbf{x}) = (\mathbf{y}^T A) \mathbf{x} = (A^T \mathbf{y})^T \mathbf{x}\]

by computing \(\mathbf{u} = \mathsf{A}(\mathbf{x})\) and \(\mathbf{v} = \mathsf{AT}(\mathbf{y})\) for random \(\mathbf{x}\) and \(\mathbf{y}\) and confirming that

\[\frac{| \mathbf{y}^T \mathbf{u} - \mathbf{v}^T \mathbf{x} |} {\max \left\{ | \mathbf{y}^T \mathbf{u} |, | \mathbf{v}^T \mathbf{x} | \right\}} < \epsilon \;.\]

If \(\mathsf{A}\) is a complex operator (with a complex input_dtype) then the test checks whether \(\mathsf{AT}\) is the Hermitian conjugate of \(\mathsf{A}\), with a test as above, but with all the \((\cdot)^T\) replaced with \((\cdot)^H\).

Parameters:

• A (LinearOperator) – Primary LinearOperator.

• AT (LinearOperator) – Adjoint LinearOperator.

• eps (Optional[float]) – Error threshold for validation of \(\mathsf{AT}\) as adjoint of \(\mathsf{AT}\). If None, the relative error is returned instead of a boolean value.

• x (Optional[Array]) – If not the default None, use the specified array instead of a random array as test vector \(\mb{x}\). If specified, the array must have shape A.input_shape.

• y (Optional[Array]) – If not the default None, use the specified array instead of a random array as test vector \(\mb{y}\). If specified, the array must have shape AT.input_shape.

• key (Optional[Array]) – Jax PRNG key. Defaults to None, in which case a new key is created.

Return type:

Union[bool, float]

Returns:

Boolean value indicating whether validation passed, or relative error of test, depending on type of parameter eps.

class scico.linop.Parallel2dProjector(im_shape, angles, x0=None, dx=None, y0=None, det_count=None)

Bases: object

Parallel ray, single axis, 2D X-ray projector.

This implementation approximates the projection of each rectangular pixel as a boxcar function (whereas the exact projection is a trapezoid). Detector pixels are modeled as bins (rather than points) and this approximation allows fast calculation of the contribution of each pixel to each bin because the integral of the boxcar is simple.

By requiring the side length of the pixels to be less than or equal to the bin width (which is assumed to be 1.0), we ensure that each pixel contributes to at most two bins, which accelerates the accumulation of pixel values into bins (equivalently, makes the linear operator sparse).

x0, dx, and y0 should be expressed in units such that the detector spacing dy is 1.0.

Parameters:

• im_shape (Tuple[int, ...]) – Shape of input array.

• angles (Union[Array, ndarray, bool_, number, bool, int, float, complex]) – (num_angles,) array of angles in radians. Viewing an (M, N) array as a matrix with M rows and N columns, an angle of 0 corresponds to summing rows, an angle of pi/2 corresponds to summing columns, and an angle of pi/4 corresponds to summing along antidiagonals.

• x0 (Union[Array, ndarray, bool_, number, bool, int, float, complex, None]) – (x, y) position of the corner of the pixel im[0,0]. By default, -im_shape / 2.

• dx (Union[Array, ndarray, bool_, number, bool, int, float, complex, None]) – Image pixel side length in x- and y-direction. Should be <= 1.0 in each dimension. By default, [1.0, 1.0].

• y0 (Optional[float]) – Location of the edge of the first detector bin. By default, -det_count / 2

• det_count (Optional[int]) – Number of elements in detector. If None, defaults to the size of the diagonal of im_shape.

project(im)

Compute X-ray projection.

class scico.linop.XRayTransform(projector)

Bases: LinearOperator

X-ray transform operator.

Wrap an X-ray projector object in a SCICO LinearOperator.

Parameters:

projector – instance of an X-ray projector object to wrap, currently the only option is Parallel2dProjector

Modules

scico.linop.abel	Abel transform LinearOperator wrapping the pyabel package.
scico.linop.optics	Optical propagator classes.
scico.linop.xray	X-ray transform classes.