pylops

name	pylops
description	Linear operators for large-scale inverse problems with matrix-free representations. Use when Claude needs to: (1) Define linear operators for forward/adjoint operations, (2) Solve inverse problems (deconvolution, imaging, tomography), (3) Apply signal processing transforms (FFT, convolution, derivatives), (4) Compose operators for complex workflows, (5) Perform regularized inversion with smoothness or sparsity constraints, (6) Process seismic or image data at scale.
version	1.0.0
author	Geoscience Skills
license	MIT
tags	["Linear Operators","Inverse Problems","Deconvolution","Signal Processing"]
dependencies	["pylops>=2.0.0","numpy","scipy"]
complements	["devito","simpeg","segyio"]
workflow_role	modelling

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

1. Never form full matrix - Use .matvec() and .rmatvec() for memory efficiency
2. Check shapes - Operators have .shape attribute like matrices
3. Verify adjoint - Always run pylops.utils.dottest() for custom operators
4. Start simple - Test on small problems before scaling up
5. 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