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

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

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