pylops

PyLops - Linear Operators Library

Quick Reference

import numpy as np
import pylops

# Create operator and apply forward/adjoint
A = pylops.FirstDerivative(n=100, dtype='float64')
y = A @ x        # Forward: y = A @ x
x_adj = A.H @ y  # Adjoint: x = A.H @ y
x_est = A / y    # Solve inverse problem

Key Classes

Class	Purpose
`LinearOperator`	Base class for all operators
`VStack/HStack`	Vertical/horizontal operator stacking
`BlockDiag`	Block diagonal operator composition

Essential Operations

Basic Operators

# Diagonal operator
D = pylops.Diagonal(np.array([1., 2., 3.]))
y = D @ x; x_adj = D.H @ y

# Derivatives
D1 = pylops.FirstDerivative(n, dtype='float64')
D2 = pylops.SecondDerivative(n, dtype='float64')
G = pylops.Gradient(dims=(64, 64), dtype='float64')

Convolution

wavelet = np.sin(np.linspace(0, 2*np.pi, 21)) * np.hanning(21)
C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=10)
y = C @ x      # Convolve
x_adj = C.H @ y  # Correlation (adjoint)

Compose and Stack Operators

# Chain: y = C @ B @ A @ x
composed = pylops.Smoothing1D(5, n) @ pylops.FirstDerivative(n) @ pylops.Identity(n)

# Stack operators
V = pylops.VStack([A, B])     # Vertical: (2n, n)
H = pylops.HStack([A, B])     # Horizontal: (n, 2n)
BD = pylops.BlockDiag([A, B]) # Block diagonal: (2n, 2n)

Solve Inverse Problems

# Simple least squares
x_est = A / y

# Normal equations
x_est = pylops.optimization.leastsquares.NormalEquationsInversion(A, None, y)

# Regularized inversion with smoothness
Reg = pylops.SecondDerivative(n)
x_est = pylops.optimization.leastsquares.RegularizedInversion(
    A, [Reg], y, epsRs=[0.1]
)

Iterative Solvers

x_lsqr = pylops.optimization.solver.lsqr(A, y, iter_lim=100)[0]
x_cgls = pylops.optimization.solver.cgls(A, y, niter=100)[0]

Sparsity-Promoting Inversion

x_l1 = pylops.optimization.sparsity.fista(A, y, niter=100, eps=0.1)[0]

Verify Adjoint (Dot Test)

A = pylops.FirstDerivative(100)
pylops.utils.dottest(A, 100, 100, verb=True)  # Dot test passed!

Common Patterns

Seismic Deconvolution

C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=len(wavelet)//2)
seismic = C @ reflectivity

# Deconvolve (regularized)
Reg = pylops.SecondDerivative(n)
reflectivity_est = pylops.optimization.leastsquares.RegularizedInversion(
    C, [Reg], seismic, epsRs=[0.01]
)

Image Denoising with TV

ny, nx = image.shape
G = pylops.Gradient(dims=(ny, nx))
x_tv = pylops.optimization.sparsity.splitbregman(
    pylops.Identity(ny*nx), G, image.ravel(),
    niter_inner=5, niter_outer=10, mu=1.0, epsRL1s=[0.1]
)[0].reshape(ny, nx)

When to Use vs Alternatives

Scenario	Recommendation
Matrix-free linear operators for large inverse problems	PyLops - purpose-built, memory efficient
Sparse matrix operations with known structure	scipy.sparse - standard, well-documented
Simple convolution/deconvolution	PyLops - clean API with `Convolve1D`
Custom operators for small problems	Custom NumPy/SciPy - no extra dependency
GPU-accelerated linear algebra	PyLops - pass CuPy arrays for automatic GPU
Seismic deconvolution or imaging operators	PyLops - rich signal processing operator library

Choose PyLops when: You need matrix-free linear operators that scale to large problems without forming explicit matrices. Its operator algebra (@, VStack, BlockDiag) and built-in solvers (LSQR, FISTA, Split Bregman) make inverse problem workflows concise.

Avoid PyLops when: Your problem is small enough for explicit matrices (use NumPy/SciPy), or you need nonlinear operators (PyLops is strictly linear).

Common Workflows

Regularized seismic deconvolution

Define wavelet array and create Convolve1D operator
Generate or load seismic trace data
Run dot test to verify operator adjoint: pylops.utils.dottest()
Set up regularization operator (e.g., SecondDerivative for smoothness)
Run RegularizedInversion(C, [Reg], data, epsRs=[eps])
Compare estimated reflectivity against true (if available)
Tune epsRs parameter: higher = smoother, lower = sharper
For sparse solutions, use pylops.optimization.sparsity.fista() instead

Tips

Never form full matrix - Use .matvec() and .rmatvec() for memory efficiency
Check shapes - Operators have .shape attribute like matrices
Verify adjoint - Always run pylops.utils.dottest() for custom operators
Start simple - Test on small problems before scaling up
Use GPU - Pass CuPy arrays for automatic GPU acceleration

References

Operators Reference - Complete list of available operators
Solvers Reference - Solver methods and configuration options

Scripts

scripts/deconvolution.py - Seismic deconvolution example