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


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


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


Returns an operator which flips the order of entries in input upside down


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.


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


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.


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))\)


Returns the transpose of a given operator \(T = A^T\)


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


Returns an operator which represents the DFT orthonormal basis


Returns an operator which represents the DCT-II orthonormal basis


Returns an operator which represents the Walsh Hadamard Transform Basis

Unions of bases


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\)


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


Computes the upper frame bound for a linear operator using full SVD



Returns the same linear operator with compiled times and trans functions


Converts a linear operator to a matrix


Converts the adjoint of a linear operator to a matrix


Converts a linear operator to a matrix in complex numbers