concrete — Concrete linear operators

Basic operators

class Identity(shape: tuple[int, ...], name: str = 'I')[source]

Bases: LinOp

Identity operator of specific shape.

Parameters:

  • shape (tuple of int) – The shape of the input and output.

  • name (str, optional) – Name of the operator.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

asmatrix(like: Array | None = None) → Array[source]

Return the matrix corresponding to the linear operator.

Applies forward to N unit vectors where N = linop.isize.

Parameters:

like (Array, optional) – If provided, use its array namespace; otherwise use float64 numpy array. See guide.

Returns:

2D array of shape (osize, isize).

Return type:

Array

Notes

Can be very heavy depending on the size of the operator.

class Diag(diag: Array, name: str = 'D')[source]

Bases: LinOp

Diagonal operator.

Parameters:

  • diag (Array) – The diagonal values. Input and output have the same shape as diag. The array namespace is inferred from this array.

  • name (str, optional) – Name of the operator.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

fwadj(point: Array) → Array[source]

Apply Aᴴ·A to point.

Parameters:

point (Array) – Input array of shape ishape.

Returns:

Output array of shape ishape.

Return type:

Array

asmatrix(like: Array | None = None) → Array[source]

Return the matrix corresponding to the linear operator.

Applies forward to N unit vectors where N = linop.isize.

Parameters:

like (Array, optional) – If provided, use its array namespace; otherwise use float64 numpy array. See guide.

Returns:

2D array of shape (osize, isize).

Return type:

Array

Notes

Can be very heavy depending on the size of the operator.

Fourier transforms

class DFT(shape: tuple[int, ...], ndim: int, name: str = 'DFT')[source]

Bases: LinOp

Discrete Fourier Transform on the last N axes.

Parameters:

  • shape (tuple of int) – The shape of the input.

  • ndim (int) – The number of last axes over which to compute the DFT.

  • name (str, optional) – Name of the operator.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

fwadj(point: Array) → Array[source]

Apply Aᴴ·A to point.

Parameters:

point (Array) – Input array of shape ishape.

Returns:

Output array of shape ishape.

Return type:

Array

class RealDFT(shape: tuple[int, ...], ndim: int, name: str = 'rDFT')[source]

Bases: LinOp

Real Discrete Fourier Transform on the last N axes.

Parameters:

  • shape (tuple of int) – The shape of the input.

  • ndim (int) – The number of last axes over which to compute the DFT.

  • name (str, optional) – Name of the operator.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

fwadj(point: Array) → Array[source]

Apply Aᴴ·A to point.

Parameters:

point (Array) – Input array of shape ishape.

Returns:

Output array of shape ishape.

Return type:

Array

Convolutions

class Conv(ir: Array, ishape: tuple[int, ...], dim: int, name: str = 'Conv')[source]

Bases: LinOp

ND convolution on the last N axes.

Does not assume a periodic or circular boundary condition.

Parameters:

  • ir (Array) – The impulse response. Must have at least dim dimensions. The array namespace is inferred from this array.

  • ishape (tuple of int) – The shape of the input. Images are on the last dim axes.

  • dim (int) – Number of last axes over which convolution applies.

  • name (str, optional) – Name of the operator.

imp_resp

The impulse response.

Type:

Array

freq_resp

The frequency response.

Type:

Array

dim

The last dim axes where convolution applies.

Type:

int

Notes

Uses FFT internally for fast computation. The forward method is equivalent to “valid” boundary condition and adjoint is equivalent to “full” boundary condition with zero filling.

Uses the array namespace of the impulse response.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

class DirectConv(ir: Array, ishape: tuple[int, ...], name: str = 'DConv')[source]

Bases: LinOp

Direct convolution.

The convolution is performed on the last N axes where N = ir.ndim.

Parameters:

  • ir (Array) – The impulse response (numpy array).

  • ishape (tuple of int) – The shape of the input array.

  • name (str, optional) – Name of the operator.

Notes

Numpy-only. Uses scipy.signal.oaconvolve (Overlap-Add method), which is generally faster than FFT-based convolution when one array is much larger than the other. Requires scipy.

property ir: Array

The impulse response.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

class FreqFilter(ir: Array, ishape: tuple[int, ...], name: str = 'Filter')[source]

Bases: Diag

Frequency filter in Fourier space.

Parameters:

  • ir (Array) – The impulse response.

  • ishape (tuple of int) – The shape of the input array (used to compute the frequency response).

  • name (str, optional) – Name of the operator.

diag

The frequency response of the filter.

Type:

Array

Notes

Almost like diagonal but assumes a complex Fourier space and is defined by an impulse response. If you have the frequency response, just use Diag.

class CircConv(imp_resp: Array, shape: tuple[int, ...], name: str = 'CConv')[source]

Bases: LinOp

Circulant (periodic) convolution.

Parameters:

  • imp_resp (Array) – The impulse response.

  • shape (tuple of int) – Shape of the input and output arrays.

  • name (str, optional) – Name of the operator.

imp_resp

The impulse response.

Type:

Array

property freq_resp: Array

The frequency response.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

fwadj(point: Array) → Array[source]

Apply Aᴴ·A to point.

Parameters:

point (Array) – Input array of shape ishape.

Returns:

Output array of shape ishape.

Return type:

Array

class Diff(axis: int, ishape: tuple[int, ...], name: str = 'Diff')[source]

Bases: LinOp

Difference operator.

Compute the first-order differences along an axis.

Parameters:

  • axis (int) – The axis along which to perform the diff.

  • ishape (tuple of int) – The shape of the input array.

  • name (str, optional) – Name of the operator.

axis

The axis along which the differences are performed.

Type:

int

forward(point: Array) → Array[source]

The forward application A·x.

This corresponds to the application of the following matrix in 1D:

-1  1  0  0
 0 -1  1  0
 0  0 -1  1
adjoint(point: Array) → Array[source]

The adjoint application Aᴴ·y.

This corresponds to the application of the following matrix in 1D:

-1  0  0
 1 -1  0
 0  1 -1
 0  0  1

Other operators

class Sampling(ishape: tuple[int, ...], index: tuple)[source]

Bases: LinOp

Sampling operator using numpy fancy indexing.

Numpy-only. Index is a tuple of index arrays as in numpy fancy indexing.

Parameters:

  • ishape (tuple of int) – The shape of the input array.

  • index (tuple of array of int) – Tuple of index arrays, one per dimension, as in numpy fancy indexing. All arrays must have the same shape, which becomes the output shape.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

class Slice(ishape: tuple[int, ...], idx: tuple)[source]

Bases: LinOp

Equivalent to obj[::2, 1, …] etc.

Parameters:

  • ishape (tuple of int) – The shape of the input array.

  • idx (index expression) – The index expression to apply (use np.index_exp to build it). The output shape is inferred as np.empty(ishape)[idx].shape.

See also

Sampling

when you have an array of indices that can handle multiple sampling of the same value.

Notes

Use np.index_exp to build the idx argument.

Examples

>>> s = Slice((10, 10), idx=np.index_exp[::2, 1])
>>> y = s.forward(np.empty((10, 10)))  # shape (5,)
>>> x = s.adjoint(y)                   # shape (10, 10)
property idx: tuple

The index expression.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

Wavelet operators

class DWT(shape: tuple[int, ...], level: int | None = None, wavelet: str = 'haar', name: str = 'DWT')[source]

Bases: LinOp

Unitary Discrete Wavelet Transform.

Parameters:

  • shape (tuple of int) – The input shape.

  • level (int, optional) – The decomposition level.

  • wavelet (str, optional) – The wavelet to use.

  • name (str, optional) – Name of the operator.

wlt

The wavelet.

Type:

str

lvl

The decomposition level.

Type:

int

Notes

NumPy-only. Uses pywt internally.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

fwadj(point: Array) → Array[source]

Apply Aᴴ·A to point.

Parameters:

point (Array) – Input array of shape ishape.

Returns:

Output array of shape ishape.

Return type:

Array

class Analysis2(shape: tuple[int, int], level: int, wavelet: str = 'haar', name: str = 'A')[source]

Bases: LinOp

2D analysis operator with stationary wavelet decomposition.

Parameters:

  • shape (tuple of (int, int)) – The input shape.

  • level (int) – The decomposition level.

  • wavelet (str, optional) – The wavelet to use.

  • name (str, optional) – Name of the operator.

Notes

NumPy-only. Uses pywt internally. The output is a 3D array where the first axis is the coefficient axis, with the approximation coefficients at index 0 and the detail coefficients at indices 1 to 3*level. The second and third axes are the spatial axes. See pywt.swt2 documentation for more details on the output format.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

cube2coeffs(point: Array) → list[source]

Return pywt coefficients from 3D array.

static coeffs2cube(coeffs: list) → Array[source]

Return 3D array from pywt coefficients.

im2coeffs(point: Array) → list[source]

Return pywt coefficients from an image array.

static coeffs2im(coeffs: list) → Array[source]

Return an image array from pywt coefficients.

cube2im(cube: Array) → Array[source]

Convert 3D coefficient cube to image-stacked array.

im2cube(im: Array) → Array[source]

Convert image-stacked array to 3D coefficient cube.

get_irs() → Array[source]

Return the impulse response of the filter bank.

get_frs() → Array[source]

Return the frequency response of the filter bank.

class Synthesis2(shape: tuple[int, int], level: int, wavelet: str = 'haar', name: str = 'S')[source]

Bases: LinOp

2D synthesis operator with stationary wavelet decomposition.

Parameters:

  • shape (tuple of (int, int)) – The input shape.

  • level (int) – The decomposition level.

  • wavelet (str, optional) – The wavelet to use.

  • name (str, optional) – Name of the operator.

forward(point: Array) → Array[source]

Returns the forward application A·x.

adjoint(point: Array) → Array[source]

Returns the adjoint application Aᴴ·y.

cube2coeffs(point: Array) → list[source]

Return pywt coefficients from 3D array.

coeffs2cube(coeffs: list) → Array[source]

Return 3D array from pywt coefficients.

im2coeffs(point: Array) → list[source]

Return pywt coefficients from image.

coeffs2im(coeffs: list) → Array[source]

Return image from pywt coefficients.

cube2im(cube: Array) → Array[source]

Convert 3D coefficient cube to image-stacked array.

im2cube(im: Array) → Array[source]

Convert image-stacked array to 3D coefficient cube.

get_irs() → Array[source]

Return the impulse response of the filter bank.

get_frs() → Array[source]

Return the frequency response of the filter bank.