# 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)**

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

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.