Linear Operators

We provide a collection of linear operators with efficient JAX based implementations that are relevant in standard signal/image processing problems. We also provide a bunch of utilities to combine and convert linear operators.

This module is inspired by pylops although the implementation approach is different.

A linear operator \(T : X \to Y\) connects a model space \(X\) to a data space \(Y\).

A linear operator satisfies following laws:

(1)\[T (x + y) = T (x) + T (y)\]

(2)\[T (\alpha x) = \alpha T(x)\]

Thus, for a general linear combination:

(3)\[T (\alpha x + \beta y) = \alpha T (x) + \beta T (y)\]

We are concerned with linear operators \(T : \mathbb{F}^n \to \mathbb{F}^m\) where \(\mathbb{F}\) is either the field of real numbers or complex numbers. \(X = \mathbb{F}^n\) is the model space and \(Y = \mathbb{F}^m\) is the data space.

Such a linear operator can be represented by a two dimensional matrix \(A\).

The forward operation is given by:

(4)\[y = A x.\]

The corresponding adjoint operation is given by:

(5)\[\hat{x} = A^H y\]

We represent a linear operator by a pair of functions times and trans. The times function implements the forward operation while the trans function implements the adjoint operation.

An inverse problem consists of computing \(x\) given \(y\) and \(A\).

1D, 2D, ND operators

• A 1D operator takes a 1D array as input and returns a 1D array as output. E.g. identity, pad, etc.

• A 2D operator takes a 2D array as input and returns a 2D array as output. E.g. conv2D, dwt2D, etc.

• An ND operator takes an ND array as input and returns an ND array as output.

The vectors may be stored using multi-dimensional arrays in memory. E.g., images are usually stored in 2-3 dimensional arrays. The operators themselves may work directly on multi-dimensions arrays. E.g. a 2D convolution operator can be applied directly to an image to result in another image.

In other words, the vectors from the model space as well as data space may be stored in memory using 1D,2D,…,ND array. They should still be treated as vectors for the purposes of this module.

This is a departure from pylops convention where input to a linear operator must be flattened into a 1D array and output needs to be reshaped again. In this library, the input and output to a 2D linear operator would be a 2D array.

axis parameter in a 1D linear operator

• A 1D linear operator may get an ND array as input.

• In this case, the axis parameter to the operator specifies the axis along which the linear operator is to be applied.

• The input ND array will be broken into slices of 1D arrays along the axis and the linear operator will be applied separately to each slice.

• Then the slices will be combined to generate the output ND array.

• E.g. if the input is a matrix then:

  • axis=0 means apply the linear operator over each column (along axis=0)

  • axis=1 means apply the linear operator over each row (along axis=1)

This is based on the convention followed by numpy.apply_along_axis.

Data types

Operator

Represents a finite linear operator \(T : A -> B\) where \(A\) and \(B\) are finite vector spaces.

Basic operators

identity(in_dim[, out_dim, axis])	Returns an identity linear operator from model space to data space
matrix(A[, axis])	Converts a matrix into a linear operator
real_matrix(A)	Converts a real matrix into a linear operator
sparse_real_matrix(A[, axis])	Converts a sparse real matrix into a linear operator
scalar_mult(alpha, n)	Returns a linear operator T such that \(T v = \alpha v\)
diagonal(d[, axis])	Returns a linear operator which mimics multiplication by a diagonal matrix
zero(in_dim[, out_dim, axis])	Returns a linear operator which maps everything to 0 vector in data space
flipud(n)	Returns an operator which flips the order of entries in input upside down
sum(n)	Returns an operator which computes the sum of a vector
dot(v[, adjoint, axis])	Returns a linear operator T such that \(T x = \langle v , x \rangle = v^H x\)
pad_zeros(n, before, after)	Adds zeros before and after a vector.
symmetrize(n)	An operator which constructs a symmetric vector by pre-pending the input in reversed order
restriction(n, I[, axis])	An operator which computes y = x[I] over an index set I
reshape(in_shape, out_shape[, order])	Returns a linear operator which reshapes vectors from model space to data space
arr2vec(shape)	Returns a linear operator which reshapes arrays to vectors
heaviside(n[, axis, normalized])	Returns a linear operator implements the Heaviside step function
inv_heaviside(n[, axis, normalized])	Returns a linear operator that computes the inverse of Heaviside/cumsum on input

Operator calculus

It is possible to combine one or more linear operators to create new linear operators. The functions in this section provide different ways to combine linear operators.

neg(A)	Returns the negative of a linear operator \(T = -A\)
scale(A, alpha)	Returns the linear operator \(T = \alpha A\) for the operator \(A\)
add(A, B)	Returns the sum of two linear operators \(T = A + B\)
subtract(A, B)	Returns a linear operator \(T = A - B\)
compose(A, B[, ignore_compatibility, shape])	Returns the composite linear operator \(T = AB\) such that \(T(x)= A(B(x))\)
transpose(A)	Returns the transpose of a given operator \(T = A^T\)
adjoint(A)	Returns the adjoint of a given operator \(T = A^H\)
hcat(A, B)	Returns the linear operator \(T = [A \, B]\)
power(A, p)	Returns the linear operator \(T = A^p\)
block_diag(operators[, axis])	Returns a block diagonal operator from 2 or more operators
windowed_op(m, T[, overlap, axis])	A wrapper to convert an operator into an overcomplete windowed operator

Signal processing operators

running_average(n, length[, axis])	Computes a running average of entries in x
fir_filter(n, h[, axis])	Implements an FIR filter defined by coeffs
convolve(n, h[, offset, axis])	Implements a convolution operator with the filter h
convolve2D(shape, h[, offset, axes])	Performs 2 dimensional convolution on the input array
convolveND(shape, h[, offset, axes])	Performs N dimensional convolution on input array

Orthonormal transforms and bases

fft(m[, n, axis, mode])	1D FFT operator
dwt(n[, wavelet, level, axis, basis])	Returns a 1D Discrete Wavelet Transform operator
dwt2D(shape[, wavelet, level, axes, basis])	Returns a 2D Discrete Wavelet Transform operator
fourier_basis(n)	Returns an operator which represents the DFT orthonormal basis
cosine_basis(n)	Returns an operator which represents the DCT-II orthonormal basis
walsh_hadamard_basis(n)	Returns an operator which represents the Walsh Hadamard Transform Basis

Unions of bases

dirac_fourier_basis(n)

Returns an operator for a two-ortho basis dictionary consisting of Dirac basis and Fourier basis

Random compressive sensing operators

gaussian_dict(key, m[, n, normalize_atoms, axis])	An operator which represents a Gaussian sensix matrix (with normalized columns)
rademacher_dict(key, m[, n, …])	An operator which represents a Rademacher sensing matrix
random_onb_dict(key, n[, axis])	An operator representing a random orthonormal basis
random_orthonormal_rows_dict(key, m[, n, axis])	An operator whose rows are orthonormal (sampled from a random orthonormal basis)
sparse_binary_dict(key, m[, n, d, …])	An operator which represents a sparse binary sensing matrix

Operators for special matrices

circulant(n, c[, axis])

Circulant matrix operator

Derivatives (finite differences)

first_derivative(n[, dx, kind])	Computes the first derivative
second_derivative(n[, dx])
tv(n[, kind, axis])	Returns a total variation linear operator for 1D signals
tv2D(shape[, kind])	Returns a total variation linear operator for 2D images

Convenience operators

These operators are technically not linear on \(\mathbb{F}^n \to \mathbb{F}^m\)

real(n)

Returns the real parts of a vector of complex numbers

Operator parts

column(T, i)	Returns the i-th column of the operator T
columns(T, indices)	Returns the i-th column of the operator T

Properties of a linear operator

These are still experimental and not efficient.

normest(operator[, max_iters, error_tolerance])	Estimates the norm of a linear operator by power method
normest_jit(operator[, max_iters, …])	Estimates the norm of a linear operator by power method
upper_frame_bound(T)	Computes the upper frame bound for a linear operator using full SVD

Utilities

jit(operator)	Returns the same linear operator with compiled times and trans functions
to_matrix(A)	Converts a linear operator to a matrix
to_adjoint_matrix(A)	Converts the adjoint of a linear operator to a matrix
to_complex_matrix(A)	Converts a linear operator to a matrix in complex numbers