PyLops

This Python library is inspired by the MATLAB Spot – A Linear-Operator Toolbox project.

Linear operators and inverse problems are at the core of many of the most used algorithms in signal processing, image processing, and remote sensing. When dealing with small-scale problems, the Python numerical scientific libraries numpy and scipy allow to perform most of the underlying matrix operations (e.g., computation of matrix-vector products and manipulation of matrices) in a simple and expressive way.

Many useful operators, however, do not lend themselves to an explicit matrix representation when used to solve large-scale problems. PyLops operators, on the other hand, still represent a matrix and can be treated in a similar way, but do not rely on the explicit creation of a dense (or sparse) matrix itself. Conversely, the forward and adjoint operators are represented by small pieces of codes that mimic the effect of the matrix on a vector or another matrix.

Luckily, many iterative methods (e.g. cg, lsqr) do not need to know the individual entries of a matrix to solve a linear system. Such solvers only require the computation of forward and adjoint matrix-vector products as done for any of the PyLops operators.

Here is a simple example showing how a dense first-order first derivative operator can be created, applied and inverted using numpy/scipy commands:

import numpy as np
from scipy.linalg import lstsq

nx = 7
x = np.arange(nx) - (nx-1)/2

D = np.diag(0.5*np.ones(nx-1), k=1) - \
    np.diag(0.5*np.ones(nx-1), k=-1)
D[0] = D[-1] = 0 # take away edge effects

# y = Dx
y = np.dot(D, x)
# x = D'y
xadj = np.dot(D.T, y)
# xinv = D^-1 y
xinv = lstsq(D, y)[0]

and similarly using PyLops commands:

from pylops import FirstDerivative

nx = 7
x = range(-(nx // 2), nx // 2 + (1 if nx % 2 else 0))

Dlop = FirstDerivative(nx, dtype='float64')

# y = Dx
y = Dlop*x
# x = D'y
xadj = Dlop.H*y
# xinv = D^-1 y
xinv = Dlop / y

Note how this second approach does not require creating a dense matrix, reducing both the memory load and the computational cost of applying a derivative to an input vector \(\mathbf{x}\). Moreover, the code becomes even more compact and espressive than in the previous case letting the user focus on the formulation of equations of the forward problem to be solved by inversion.

Terminology

A common terminology is used within the entire documentation of PyLops. Every linear operator and its application to a model will be referred to as forward model (or operation)

\[\mathbf{y} = \mathbf{A} \mathbf{x}\]

while its application to a data is referred to as adjoint modelling (or operation)

\[\mathbf{x} = \mathbf{A}^H \mathbf{y}\]

where \(\mathbf{x}\) is called model and \(\mathbf{y}\) is called data. The operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{F}^n\) effectively maps a vector of size \(m\) in the model space to a vector of size \(n\) in the data space, conversely the adjoint operator \(\mathbf{A}^H:\mathbb{F}^n \to \mathbb{F}^m\) maps a vector of size \(n\) in the data space to a vector of size \(m\) in the model space. As linear operators mimics the effect a matrix on a vector we can also loosely refer to \(m\) as the number of columns and \(n\) as the number of rows of the operator.

Ultimately, solving an inverse problems accounts to removing the effect of \(\mathbf{A}\) from the data \(\mathbf{y}\) to retrieve the model \(\mathbf{x}\).

For a more detailed description of the concepts of linear operators, adjoints and inverse problems in general, you can head over to one of Jon Claerbout’s books such as Basic Earth Imaging.

Implementation

PyLops is build on top of the scipy class scipy.sparse.linalg.LinearOperator.

This class allows in fact for the creation of objects (or interfaces) for matrix-vector and matrix-matrix products that can ultimately be used to solve any inverse problem of the form \(\mathbf{y}=\mathbf{A}\mathbf{x}\).

As explained in the scipy LinearOperator official documentation, to construct a scipy.sparse.linalg.LinearOperator, a user is required to pass appropriate callables to the constructor of this class, or subclass it. More specifically one of the methods _matvec and _matmat must be implemented for the forward operator and one of the methods _rmatvec or _adjoint may be implemented to apply the Hermitian adjoint. The attributes/properties shape (pair of integers) and dtype (may be None) must also be provided during __init__ of this class.

Any linear operator developed within the PyLops library follows this philosophy. As explained more in details in Implementing new operators section, a linear operator is created by subclassing the scipy.sparse.linalg.LinearOperator class and _matvec and _rmatvec are implemented.

History

PyLops was initially written by Equinor It is a flexible and scalable python library for large-scale optimization with linear operators that can be tailored to our needs, and as contribution to the free software community.