Inverse problems, regularization and seismic imaging

Before this lesson: linear algebra 08 (SVD and pseudoinverse), linear algebra 12 (least squares), linear algebra 13 (conditioning) and lesson 03 (the adjoint method). Helpful: optimization 10 (Gauss–Newton) and the probability track (Bayes' rule with Gaussians). Time: about 100 minutes, plus the notebook.

---

Why this matters #

• Seismology. Earthquake location, travel-time tomography and full-waveform inversion (FWI) are inverse problems. FWI images the Earth's interior, from oil and gas reservoirs to the whole mantle, by fitting complete seismograms.
• Weather and climate. Data assimilation (3D-Var and 4D-Var, lesson 06) is a regularized inverse problem: find the atmospheric state from sparse observations, with a forecast as the prior.
• Imaging. Medical CT and MRI, deblurring of images and telescope data, and many more.
• Machine learning. Ridge regression is Tikhonov regularization. Gaussian-process regression (lesson 07) is the same Bayesian computation in function space. Learned regularizers, unrolled networks and diffusion-model priors are active research on inverse problems.

1. Predict first #

2. The picture #

Forward and inverse problems #

A forward map $F$ sends a model $m$ (what we want to know) to data $d$ (what we can measure): $d = F(m) + e$, with noise $e$. When $F$ is linear we write $F(m) = Gm$ with a matrix $G$.

problem	model $m$	data $d$	forward map
earthquake location	hypocentre $(x, y, z)$ and origin time $t_0$	arrival times at stations	$t_i = t_0 + \text{distance}_i / v$ (nonlinear)
travel-time tomography	slowness $1/c$ in each cell	travel times of rays	$G_{ij}$ = length of ray $i$ in cell $j$ (linear)
full-waveform inversion	velocity $c(x)$ everywhere	seismograms	solve the wave equation (nonlinear)
deblurring	sharp signal	blurred, noisy signal	convolution with a blur (linear)

Hadamard's conditions. A problem is well-posed if (1) a solution exists, (2) it is unique, and (3) it depends continuously on the data (stability). A problem that fails one of them is ill-posed. Inverse problems usually fail (3), often (2), and with noise also (1): noisy data are not exactly $F(m)$ for any $m$. After discretization the solution always depends continuously on $d$, but the dependence can be extremely sensitive: the matrix is ill-conditioned (linear algebra 13).

Why undoing a blur fails #

A blur averages neighbouring values, so it damps fine details (high frequencies) strongly and broad features (low frequencies) weakly. In the SVD $G = U\Sigma V^\top$ (linear algebra 08), the right singular vectors $v_i$ are patterns of increasing frequency, and the singular values $\sigma_i$ fall quickly. The noise has the same size in every direction $u_i$. The naive solution divides each data component by $\sigma_i$:

Where $\sigma_i$ is tiny, the noise part $u_i^\top e$ is divided by a tiny number and dominates. Regularization keeps the components with large $\sigma_i$ and damps the rest.

Interactive figure: loading…

3. Definitions and main results #

Linear inverse problems and the SVD #

Let $G \in \R^{n_d\times n_m}$ have the SVD $G = \sum_i \sigma_i u_iv_i^\top$ with $\sigma_1 \ge \sigma_2 \ge \dots > 0$ (for the nonzero ones). The minimum-norm least-squares solution (the pseudoinverse solution, linear algebra 08 and 12) is

Noise amplification. Let $d = Gm_\text{true} + e$ with noise $e \sim \mathcal N(0, \sigma^2 I)$. Here $\sigma$ (without index) is the noise level. Then

Because the $u_i$ are orthonormal, each $u_i^\top e$ is again $\mathcal N(0, \sigma^2)$. So the error along $v_i$ has standard deviation $\sigma/\sigma_i$. The condition number $\kappa = \sigma_1/\sigma_\text{min}$ bounds how much a relative error in $d$ can grow in $m$.

The discrete Picard condition (Hansen 1990). The naive coefficients $u_i^\top d/\sigma_i$ stay bounded only if the exact data coefficients $|u_i^\top d_\text{exact}|$ decay faster than $\sigma_i$, on average. Real data contain noise, and $|u_i^\top e| \approx \sigma$ for every $i$. So in the Picard plot, $|u_i^\top d|$ follows the decay of $\sigma_i$ and then levels off at the noise level. Every component after that point is mostly noise, and dividing it by $\sigma_i$ makes it explode.

Regularization: truncated SVD and Tikhonov #

Truncated SVD (TSVD) keeps the first $k$ components: $m_k = \sum_{i\le k}(u_i^\top d/\sigma_i)\,v_i$.

Tikhonov regularization (Tikhonov 1963) adds a penalty on the size of the model:

Theorem 1 (filter factors). In the SVD basis,

Proof.

1. The objective is a convex quadratic. Its gradient $2G^\top(Gm - d) + 2\alpha^2 m$ is zero at the minimum, so $(G^\top G + \alpha^2 I)m_\alpha = G^\top d$. The matrix is positive definite for $\alpha > 0$, so the solution is unique.
2. With $G = U\Sigma V^\top$ (full $V$, orthogonal): $G^\top G + \alpha^2 I = V(\Sigma^\top\Sigma + \alpha^2 I)V^\top$ and $G^\top d = V\Sigma^\top U^\top d$.
3. So $m_\alpha = V(\Sigma^\top\Sigma + \alpha^2 I)^{-1}\Sigma^\top U^\top d$. The diagonal matrix in the middle has entries $\sigma_i/(\sigma_i^2 + \alpha^2) = f_i/\sigma_i$. $\blacksquare$

What the filter factors do.

• $\sigma_i \gg \alpha$: $f_i \approx 1$, the component is kept.
• $\sigma_i \ll \alpha$: $f_i \approx \sigma_i^2/\alpha^2$, the component is removed.
• $\sigma_i = \alpha$: $f_i = \tfrac12$. So $\alpha$ is the singular value where the filter switches.
• The noise gain $f_i/\sigma_i = \sigma_i/(\sigma_i^2 + \alpha^2)$ is at most $1/(2\alpha)$ (its maximum, at $\sigma_i = \alpha$). So noise can be amplified by at most $1/(2\alpha)$, instead of $1/\sigma_\text{min}$.

TSVD uses $f_i \in \{0, 1\}$; Tikhonov switches smoothly. Larger $\alpha$ means less noise in $m_\alpha$ (less variance) but more lost detail (more bias).

Choosing $\alpha$.

• Discrepancy principle (Morozov 1966). A good model should fit the data only as well as the noise allows: choose $\alpha$ so that $\|Gm_\alpha - d\| \approx \|e\| \approx \sqrt{n_d}\,\sigma$. It needs the noise level, and it tends to choose an $\alpha$ slightly too large (slightly too smooth).
• L-curve (Hansen 1992). Plot $\log\|m_\alpha\|$ against $\log\|Gm_\alpha - d\|$ for many $\alpha$. The curve has the shape of an L. On the vertical part, a smaller $\alpha$ fits a little more noise and the model norm explodes; on the horizontal part, a larger $\alpha$ removes signal and the residual grows. Choose the corner. It does not need the noise level.
• Generalized cross-validation (GCV) is a third common choice; it predicts each datum from the others.

Interactive figure: loading…

General-form Tikhonov and sparsity

Instead of $\alpha^2\|m\|^2$ one often penalizes $\alpha^2\|Lm\|^2$ with a derivative matrix $L$, so smooth models are preferred, not small ones. The solution uses the generalized SVD of $(G, L)$. With an $\ell_1$ penalty, $\alpha\|Lm\|_1$ (total variation when $L$ is a first difference), sharp edges survive. The problem stays convex but is no longer quadratic; proximal methods solve it (optimization 09).

The Bayesian view: Tikhonov is MAP #

Model the noise and the prior as Gaussians: $e \sim \mathcal N(0, \sigma^2 I)$ and $m \sim \mathcal N(0, \tau^2 I)$. By Bayes' rule the posterior density is

Theorem 2 (Tikhonov = MAP). The maximum a posteriori (MAP) estimate is the Tikhonov solution with $\alpha^2 = \sigma^2/\tau^2$. Equivalently, Tikhonov with weight $\alpha$ is MAP for the prior $m \sim \mathcal N(0, \alpha^{-2}\sigma^2 I)$. The posterior is Gaussian with mean $m_\alpha$ and covariance

Proof.

1. Maximizing $p(m\mid d)$ is minimizing $-2\sigma^2\log p(m\mid d) = \|Gm - d\|^2 + (\sigma^2/\tau^2)\|m\|^2 + \text{const}$: the Tikhonov objective with $\alpha^2 = \sigma^2/\tau^2$.
2. The exponent is a quadratic in $m$ with Hessian $G^\top G/\sigma^2 + I/\tau^2$. A density of the form $\exp(-\tfrac12(m - \mu)^\top P(m - \mu))$ is Gaussian with covariance $P^{-1}$; its mean is the minimizer, $m_\alpha$. Multiply $P$ by $\sigma^2$ and use $\sigma^2/\tau^2 = \alpha^2$ for the second form. $\blacksquare$

In the SVD basis the posterior variance along $v_i$ is $\sigma^2/(\sigma_i^2 + \alpha^2)$. For large $\sigma_i$ it is about $\sigma^2/\sigma_i^2$: the data determine this direction. For small $\sigma_i$ it is about $\sigma^2/\alpha^2 = \tau^2$, the prior variance: the data say nothing about this direction, and we keep the prior.

General Gaussian prior and noise. With $m \sim \mathcal N(m_b, B)$ and $e \sim \mathcal N(0, R)$, MAP minimizes

This is exactly the 3D-Var cost of data assimilation (lesson 06), with background state $m_b$ and observation operator $G$. Its posterior mean, written with kernel matrices, is the posterior mean of Gaussian-process regression (lesson 07).

Nonlinear problems: Gauss–Newton and earthquake location #

For a nonlinear forward map we minimize the misfit $J(m) = \tfrac12\|F(m) - d\|^2$. Let $r = F(m) - d$ be the residual and $J_F$ the Jacobian of $F$ at $m$. Gauss–Newton linearizes $F(m + \delta) \approx F(m) + J_F\delta$ and solves the linear least-squares problem for the step:

It uses $J_F^\top J_F$ in place of the Hessian of $J$, which drops the second derivatives of $F$ (optimization 10). Adding $\alpha^2 I$ to $J_F^\top J_F$ gives the Levenberg–Marquardt method, a regularized step.

Earthquake location (Geiger's method, 1912). Stations at known positions $s_i$ record the arrival time $t_i$ of the first wave. In a medium with constant speed $v$,

where $x$ is the hypocentre and $t_0$ the origin time ($m$ has 4 entries in 3D; the figure uses 2D and 3 entries). Row $i$ of the Jacobian is $\big((x - s_i)^\top/(v\|x - s_i\|),\ 1\big)$: a unit vector from the station to the source, divided by $v$, and a 1 for $t_0$.

Station geometry. If all stations lie on one side, the unit vectors are almost parallel, the rows of $J_F$ are almost dependent, and $J_F^\top J_F$ is ill-conditioned. The linearized posterior covariance $\sigma_t^2(J_F^\top J_F)^{-1}$ (for timing noise $\sigma_t$) is then a long ellipse: the source can move along one direction (often away from the stations, with a later or earlier $t_0$) almost without changing the arrival times. Depth and origin time trade off in the same way when all stations are far away.

Interactive figure: loading…

The adjoint-state method #

Many forward maps are given by a large equation for a state $u$ (a wavefield, a temperature):

for example $j(u) = \tfrac12\|Pu - d\|^2$, where $P$ samples the state at the receivers. We need $\nabla_m J$ for millions of parameters.

Theorem 3 (adjoint-state gradient). Let $A(m)$ be invertible and differentiable in $m$. Let the adjoint state $\lambda$ solve

Then

for every $k$.

Proof (implicit differentiation, as in calculus 09).

1. Differentiate $A(m)u(m) = s$ with respect to $m_k$: $\frac{\partial A}{\partial m_k}u + A\frac{\partial u}{\partial m_k} = 0$, so $\frac{\partial u}{\partial m_k} = -A^{-1}\frac{\partial A}{\partial m_k}u$.
2. Chain rule: $\frac{\partial J}{\partial m_k} = \nabla_u j^\top\frac{\partial u}{\partial m_k} = -\nabla_u j^\top A^{-1}\frac{\partial A}{\partial m_k}u$.
3. The row vector $-\nabla_u j^\top A^{-1}$ is $\lambda^\top$, because $\lambda = -A^{-\top}\nabla_u j$. So $\frac{\partial J}{\partial m_k} = \lambda^\top\frac{\partial A}{\partial m_k}u$. $\blacksquare$

Cost. One forward solve for $u$ and one adjoint solve for $\lambda$, then a cheap product for each $k$. The number of solves does not depend on the number of parameters. With several sources, the cost is one forward and one adjoint solve per source.

The same result with a Lagrangian

1. Define $\mathcal L(u, \lambda, m) = j(u) + \lambda^\top(A(m)u - s)$. When $u = u(m)$ solves the state equation, $\mathcal L = J(m)$ for every $\lambda$.
2. $\frac{dJ}{dm_k} = \frac{\partial\mathcal L}{\partial u}\frac{\partial u}{\partial m_k} + \frac{\partial\mathcal L}{\partial m_k}$, with $\frac{\partial\mathcal L}{\partial u} = \nabla_u j^\top + \lambda^\top A$.
3. Choose $\lambda$ so that $\frac{\partial\mathcal L}{\partial u} = 0$: this is $A^\top\lambda = -\nabla_u j$. The unknown $\partial u/\partial m_k$ then drops out.
4. What remains is $\frac{\partial\mathcal L}{\partial m_k} = \lambda^\top\frac{\partial A}{\partial m_k}u$. $\blacksquare$

Lesson 03 is the special case where $A(m)u = s$ is the whole time-stepping scheme: the matrix is block lower triangular in time, so $A^\top$ is block upper triangular, and solving with it runs backward in time.

Full-waveform inversion #

In 1D, the acoustic wave equation with squared slowness $m(x) = 1/c(x)^2$ and a source at $x_s$ is

A receiver at $x_r$ records $u(x_r, t)$. FWI minimizes the waveform misfit

Theorem 4 (FWI gradient). Let $\lambda$ solve the adjoint wave equation with the residual as a source at the receiver and final conditions:

Then

Why this is true. This is Theorem 3 with $A(m)u = m\,u_{tt} - u_{xx}$. The derivative of $A(m)u$ with respect to $m(x)$ is $u_{tt}(x, \cdot)$, and the pairing $\lambda^\top(\dots)$ is an integral over time. The wave operator is symmetric in space, and its transpose in time turns the initial conditions into final conditions: the adjoint equation runs backward. Because the wave equation is time-reversible, "the residual played backward from the receiver" is a wave that travels back into the model.

Derivation with the Lagrangian

1. $\mathcal L = J + \int_0^T\!\!\int \lambda\,(m\,u_{tt} - u_{xx} - w\,\delta_{x_s})\,dx\,dt$.
2. Perturb $u$ by $\delta u$ with $\delta u = \delta u_t = 0$ at $t = 0$. Integrate by parts twice in $t$ and twice in $x$ (the boundary terms in $x$ vanish for a large domain or for matching boundary conditions): $\int\!\!\int \lambda\,(m\,\delta u_{tt} - \delta u_{xx}) = \int\!\!\int (m\,\lambda_{tt} - \lambda_{xx})\,\delta u + \big[\int m(\lambda\,\delta u_t - \lambda_t\,\delta u)\,dx\big]_{t=T}$.
3. The term from $J$ is $\int (u(x_r, t) - d(t))\,\delta u(x_r, t)\,dt = \int\!\!\int (u - d)\,\delta(x - x_r)\,\delta u$.
4. Setting all terms with $\delta u$ to zero gives the adjoint equation and the final conditions $\lambda = \lambda_t = 0$ at $t = T$.
5. What remains is $\partial\mathcal L/\partial m(x) = \int_0^T\lambda\,u_{tt}\,dt$. $\blacksquare$

We checked the discrete version of this formula (leapfrog time stepping, the discrete adjoint run backward) against finite differences of the misfit; they agree to 8 digits. The notebook repeats the check.

What the gradient looks like. At each point $x$ the gradient is the zero-lag cross-correlation of the forward field's acceleration with the adjoint field. In a reflection setting (source and receiver at the surface) the residual contains the reflected waves that the starting model misses. Sent back into the model, they meet the forward field exactly at the depth of the reflector. So the first gradient is an image of the reflectors: this is the link between FWI and reverse-time migration. Also, $10^7$ unknowns need only two simulations per source (forward and adjoint).

Interactive figure: loading…

Cycle skipping #

FWI minimizes $J$ with gradient methods, so it finds a local minimum. Consider one arrival: the predicted wavelet is the observed wavelet shifted in time by $\Delta$. For a wavelet of period $P$, the misfit $J(\Delta)$ has local minima near $\Delta = \pm P, \pm 2P, \dots$, where the predicted wavelet matches the wrong cycle of the observed one. If the starting model predicts the arrival wrong by more than about half a period, gradient descent moves toward the wrong cycle: this is cycle skipping. For the Ricker wavelet with peak frequency $f$ we use here, the basin of the correct minimum ends at $|\Delta| \approx 0.43/f$, a little less than half of the period $1/f$.

Remedies.

• Low frequencies first (Bunks et al. 1995). The period of a low-frequency wavelet is long, so the same time error is a small fraction of a period. Invert low-pass filtered data first, then add higher frequencies step by step.
• A good starting model, for example from travel-time tomography.
• Other misfits that compare shapes instead of point values, for example optimal-transport distances (Engquist and Froese 2014) or envelopes.

Interactive figure: loading…

Machine learning and inverse problems

• Learned regularizers. Replace $\alpha^2\|m\|^2$ by a penalty learned from examples of good models, for example an adversarially trained critic (Lunz, Öktem and Schönlieb 2018).
• Unrolled optimization. Take $K$ iterations of a gradient or proximal method, give each iteration learnable parameters, and train the whole chain end to end (LISTA, Gregor and LeCun 2010; learned primal-dual, Adler and Öktem 2018). The forward map and its adjoint stay inside the network.
• Generative priors. Plug-and-play methods use a denoiser as the prior step (Venkatakrishnan et al. 2013); diffusion posterior sampling uses a pretrained diffusion model as the prior and the gradient of the data misfit for guidance (Chung et al. 2023; lessons 12–13 for the diffusion side).
• Differentiable simulators. Write the forward model in an automatic-differentiation framework, and the adjoint-state gradient comes for free (the discrete adjoint of lesson 03). This also lets a neural network inside the simulator be trained through the physics.
• Amortized and simulation-based inference. Train a network on simulated pairs $(m, d)$ to output the posterior $p(m \mid d)$ for any new $d$ in one pass (Cranmer, Brehmer and Louppe 2020).
• In seismology, deep networks already pick the arrival times that earthquake location uses (for example PhaseNet, Zhu and Beroza 2019).

Answers to section 1. (1) $10^3$: the noise component $10^{-3}$ is divided by $\sigma_i = 10^{-6}$. (2) One half: $f_i = \alpha^2/(\alpha^2 + \alpha^2) = \tfrac12$. (3) Two: one forward and one adjoint wave simulation, whatever the number of unknowns (Theorem 3).

4. Worked example #

(a) Regularization in two dimensions. Let $v_1 = (1, 1)/\sqrt2$ and $v_2 = (1, -1)/\sqrt2$, and

$G$ is symmetric with eigenvalues $1$ and $0.01$, so its SVD has $u_i = v_i$, $\sigma_1 = 1$, $\sigma_2 = 0.01$, and $\kappa = 100$. Take $m_\text{true} = (1, 1)$ and noise $e = (0.01, -0.01)$.

1. Data: $Gm_\text{true} = (1, 1)$ (it lies along $v_1$), so $d = (1.01, 0.99)$.
2. Coefficients: $u_1^\top d = (1.01 + 0.99)/\sqrt2 = \sqrt2$ and $u_2^\top d = (1.01 - 0.99)/\sqrt2 = 0.02/\sqrt2$.
3. Naive solution: $\frac{\sqrt2}{1}v_1 + \frac{0.02/\sqrt2}{0.01}v_2 = (1, 1) + (1, -1) = (2, 0)$. A noise of size $0.014$ caused an error of size $\sqrt2 \approx 1.41$: amplification $1/\sigma_2 = 100$.
4. Tikhonov with $\alpha = 0.1$: $f_1 = \frac{1}{1 + 0.01} \approx 0.990$ and $f_2 = \frac{0.0001}{0.0001 + 0.01} = \frac{1}{101} \approx 0.0099$.
5. $m_\alpha = 0.990\,(1, 1) + 0.0099\,(1, -1) = (1.000, 0.980)$. The error is $0.0198$ instead of $1.41$.

The regularized solution lost almost nothing along $v_1$ (bias $1\%$) and removed almost all of the noise along $v_2$.

(b) The adjoint-state method by hand. Let

at $m = (0, 0)$.

1. Forward solve: $A = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}$, $A^{-1} = \frac13\begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$, so $u = A^{-1}s = (\tfrac23, \tfrac13)$.
2. Residual: $\nabla_u j = u - d = (\tfrac16, -\tfrac16)$.
3. Adjoint solve ($A$ is symmetric): $\lambda = -A^{-1}(\tfrac16, -\tfrac16) = -\tfrac13(\tfrac16, -\tfrac16) = (-\tfrac{1}{18}, \tfrac{1}{18})$.
4. Derivatives of the operator: $\partial A/\partial m_1 = e_1e_1^\top$, so $\frac{\partial A}{\partial m_1}u = (u_1, 0)$; likewise $\frac{\partial A}{\partial m_2}u = (0, u_2)$.
5. Gradient: $\partial J/\partial m_1 = \lambda_1u_1 = -\tfrac{1}{18}\cdot\tfrac23 = -\tfrac{1}{27}$ and $\partial J/\partial m_2 = \lambda_2u_2 = \tfrac{1}{18}\cdot\tfrac13 = \tfrac{1}{54}$.

Check: central finite differences of $J$ give $-0.0370370$ and $0.0185185$ ✓. Two linear solves gave both derivatives; with $10^7$ parameters it would still be two solves.

5. In code #

NumPy checks the deblurring numbers, the discrepancy principle and the adjoint-state example. The notebook has these cells and the exercises.

import numpy as np

# 1. A 1D blur (deconvolution) problem: d = G m + noise
n = 64
x = (np.arange(n) + 0.5) / n
w = 0.03                                            # width of the Gaussian blur
G = np.exp(-(x[:, None] - x[None, :]) ** 2 / (2 * w**2)) / (np.sqrt(2 * np.pi) * w * n)
m_true = ((x > 0.15) & (x < 0.35)) + 0.8 * np.exp(-((x - 0.68) / 0.07) ** 2)
sigma = 1e-3                                        # noise level
d = G @ m_true + sigma * np.random.default_rng(0).normal(size=n)

U, s, Vt = np.linalg.svd(G)
print(f"{s[0]:.3f} {s[40]:.1e} {s[0] / s[-1]:.1e}")  # 0.996 1.0e-03 3.3e+07
m_naive = np.linalg.solve(G, d)
print(f"{np.linalg.norm(m_naive - m_true):.1e}")    # 5.1e+04: the naive inverse is useless

# 2. Tikhonov through the filter factors f_i = s_i^2 / (s_i^2 + alpha^2)
beta = U.T @ d                                      # coefficients u_i^T d
def tikhonov(alpha):
    f = s**2 / (s**2 + alpha**2)
    return Vt.T @ (f * beta / s)

for alpha in [1e-4, 1e-3, 1e-2, 1e-1]:
    m = tikhonov(alpha)
    print(f"{alpha:.0e} residual {np.linalg.norm(G @ m - d):.2e} error {np.linalg.norm(m - m_true):.2f}")
# 1e-04 residual 3.09e-03 error 8.60
# 1e-03 residual 3.64e-03 error 1.17
# 1e-02 residual 4.91e-03 error 0.54
# 1e-01 residual 4.99e-02 error 0.69

# 3. Discrepancy principle: choose alpha with ||G m_alpha - d|| = sqrt(n) sigma (bisection in log alpha)
lo, hi = -6.0, 0.0
for _ in range(60):
    mid = (lo + hi) / 2
    r = np.linalg.norm(G @ tikhonov(10**mid) - d)
    lo, hi = (mid, hi) if r < np.sqrt(n) * sigma else (lo, mid)
a_dp = 10**lo
print(f"alpha {a_dp:.3f} error {np.linalg.norm(tikhonov(a_dp) - m_true):.2f}")   # alpha 0.025 error 0.57

# 4. Adjoint-state gradient for A(m) u = s, J = |u - d|^2 / 2, checked with finite differences
def A(m): return np.array([[2 + m[0], -1.0], [-1.0, 2 + m[1]]])
src, dat = np.array([1.0, 0.0]), np.array([0.5, 0.5])
def J(m): u = np.linalg.solve(A(m), src); return 0.5 * np.sum((u - dat) ** 2)
m0 = np.zeros(2)
u = np.linalg.solve(A(m0), src)                     # forward solve
lam = np.linalg.solve(A(m0).T, -(u - dat))          # adjoint solve: A^T lam = -grad_u J
grad = np.array([lam[0] * u[0], lam[1] * u[1]])     # lam^T (dA/dm_k) u, with dA/dm_k = e_k e_k^T
fd = [(J(m0 + 1e-6 * e) - J(m0 - 1e-6 * e)) / 2e-6 for e in np.eye(2)]
print(grad, np.round(fd, 8))                        # [-0.03703704  0.01851852] [-0.03703704  0.01851852]

The residual stays near $3\cdot10^{-3}$ for small $\alpha$ because the model can fit almost all of the data, noise included; the target of the discrepancy principle is $\sqrt{64}\cdot10^{-3} = 8\cdot10^{-3}$. The smallest error in this list is at $\alpha = 10^{-2}$; the discrepancy principle chooses a slightly larger $\alpha$, as expected.

6. Where it shows up #

• Global and regional seismology. Adjoint tomography fits full seismograms of many earthquakes with the adjoint-state gradient (Tromp, Tape and Liu 2005). Global models of the mantle have been built this way from hundreds of earthquakes (for example GLAD-M15, Bozdağ et al. 2016).
• Exploration geophysics. FWI is a standard step in building velocity models for reservoirs (review: Virieux and Operto 2009). The first gradient is closely related to reverse-time migration.
• Earthquake catalogs. Every earthquake catalog comes from solving a location problem for each event. Relative methods such as double-difference location (Waldhauser and Ellsworth 2000) reduce the effect of an imperfect velocity model.
• Weather and climate. 3D-Var and 4D-Var are regularized inverse problems with the forecast as the prior (lesson 06); 4D-Var uses the adjoint of the forecast model, as FWI uses the adjoint of the wave equation.
• Machine learning. Ridge regression is Tikhonov regularization with the data matrix as $G$. Weight decay is the same penalty on network weights. Learned reconstruction methods for CT and MRI combine the forward operator with trained networks.

7. Common mistakes #

• Inverting an ill-conditioned matrix directly. np.linalg.solve(G, d) returns an answer, but it is dominated by amplified noise. Look at the singular values and the Picard plot first.
• Choosing $\alpha$ by the best fit to the data. The residual always falls as $\alpha \to 0$; the best fit is the noisiest model. Use the discrepancy principle, the L-curve, GCV or a held-out set.
• Mixing up $\alpha$ and $\alpha^2$. With the objective $\|Gm - d\|^2 + \alpha^2\|m\|^2$ the filter factor switches at $\sigma_i = \alpha$ and the prior variance is $\sigma^2/\alpha^2$. Libraries often write $\lambda\|m\|^2$ with $\lambda = \alpha^2$.
• A wrong sign or a missing time reversal in the adjoint. The adjoint wave equation has final conditions and runs backward. Check every adjoint gradient against finite differences on a small problem before using it.
• Starting FWI with high frequencies or a poor model. The misfit is not convex; cycle skipping gives a wrong model with a small gradient. Start from low frequencies and a smooth, travel-time-based model.
• Reading the regularized solution as the truth. Components with small $\sigma_i$ come from the prior, not from the data. Report the posterior uncertainty or the resolution.

8. Check yourself #

1. State Hadamard's three conditions. Which one fails for deblurring, and why does a discrete blur matrix still "work" in exact arithmetic?
   Answer

   Existence, uniqueness, continuous dependence on the data. Stability fails: tiny changes in the data cause huge changes in the solution. A discrete blur matrix is usually invertible, so the solution depends continuously on $d$, but the constant (the condition number, here about $3\cdot10^7$) is so large that rounding and noise dominate.
2. Write the Tikhonov solution in the SVD basis. What is the largest possible amplification of a data error, and for which singular value does it happen?
   Answer

   $m_\alpha = \sum_i \frac{\sigma_i}{\sigma_i^2 + \alpha^2}(u_i^\top d)\,v_i$. The gain $\sigma/(\sigma^2 + \alpha^2)$ is largest at $\sigma = \alpha$ (set its derivative $(\alpha^2 - \sigma^2)/(\sigma^2 + \alpha^2)^2$ to zero), where it equals $1/(2\alpha)$.
3. The noise level is $\sigma = 0.01$ and there are $n_d = 400$ data. What residual norm does the discrepancy principle aim for?
   Answer

   $\|Gm_\alpha - d\| \approx \sqrt{n_d}\,\sigma = 20 \cdot 0.01 = 0.2$, the expected size of the noise vector.
4. Give the posterior covariance of the Tikhonov problem with noise $\mathcal N(0, \sigma^2 I)$, and its variance along $v_i$. What happens for $\sigma_i \to 0$?
   Answer

   $\Sigma_\text{post} = \sigma^2(G^\top G + \alpha^2 I)^{-1}$; along $v_i$ the variance is $\sigma^2/(\sigma_i^2 + \alpha^2)$. For $\sigma_i \to 0$ it tends to $\sigma^2/\alpha^2$, the prior variance: the data do not constrain that direction.
5. Four stations lie on a line, and an earthquake happens off the line. What can go wrong in the location, and why?
   Answer

   The mirror image of the source across the line has exactly the same distances to all stations, so it predicts the same arrival times: the problem has two solutions (uniqueness fails). Gauss–Newton converges to the one on the side of the starting guess. Near the line, the misfit is also very flat across it.
6. For $A(m) = L + \operatorname{diag}(m)$ with a fixed matrix $L$, give the adjoint-state gradient of $J = \tfrac12\|Pu - d\|^2$.
   Answer

   Solve $Au = s$ and $A^\top\lambda = -P^\top(Pu - d)$. Since $\partial A/\partial m_k = e_ke_k^\top$, $\partial J/\partial m_k = \lambda_k u_k$: the gradient is the elementwise product $\lambda \odot u$.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. State Hadamard's three conditions for a well-posed problem.

   Show answer

2. Write the naive SVD solution of $d = Gm + e$ and say why it fails.

   Show answer

3. What are the Tikhonov filter factors?

   Show answer

4. What does the discrete Picard condition say?

   Show answer

5. For which prior is Tikhonov regularization the MAP estimate, with noise $\mathcal N(0, \sigma^2 I)$?

   Show answer

6. What residual norm does the discrepancy principle aim for?

   Show answer

7. Write one Gauss–Newton step for $\min \tfrac12\|F(m) - d\|^2$.

   Show answer

8. Give the adjoint-state gradient for $A(m)u = s$ and $J = j(u)$, and its cost.

   Show answer

9. Give the FWI gradient for $m\,u_{tt} - u_{xx} = s$.

   Show answer

10. When does FWI cycle-skip?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why not just invert the matrix?" — Small singular values amplify the noise by $1/\sigma_i$. Show the singular values and the Picard plot, then regularize.
• "Explain ridge regression as a Bayesian method." — It is MAP with a Gaussian prior $\mathcal N(0, (\sigma^2/\alpha^2)I)$ on the weights; the posterior covariance is $\sigma^2(X^\top X + \alpha^2 I)^{-1}$.
• "How do you choose the regularization weight?" — Discrepancy principle if the noise level is known; L-curve or GCV otherwise; cross-validation on held-out data in ML.
• "How do you get the gradient for millions of parameters in a PDE-constrained problem?" — The adjoint-state method: one forward and one adjoint solve; derive it with a Lagrangian; check it against finite differences.
• "What is cycle skipping and how do you avoid it?" — A local minimum when the arrival is off by more than half a period; use low frequencies first, a good starting model, or a misfit that is less oscillatory.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Build the blur matrix, see the naive inverse fail, and fix it with truncated SVD and Tikhonov.
2. Choose $\alpha$ with the discrepancy principle and with the corner of the L-curve; compare with the best $\alpha$.
3. Locate an earthquake with Gauss–Newton from 5 stations, and compare the uncertainty ellipses of two station geometries.
4. The adjoint-state gradient for a 1D problem $(L + \operatorname{diag}(m))u = s$ with receivers, checked with finite differences.
5. Misfit curves of a 1D wave simulation for a low and a high frequency: count the local minima (cycle skipping).

11. Watch and read #

Regularization Methods - Part 1: Introduction to Inverse ProblemsFranks Mathematics · 27 min

Ill-posedness and the SVD view of a linear inverse problem.

Regularization Methods - Part 2: Tikhonov RegularizationFranks Mathematics · 32 min

Tikhonov regularization and its filter factors.

Chap 2: Hadamard & Picard Conditions, Singular Value Expansion, Naive Reconstruction - 1DTUdk · 21 min

The singular value expansion, the Picard condition and why the naive reconstruction fails.

Lectures on FWI 1: overviewMrCrazyPhysicist · 23 min

What full-waveform inversion is and why cycle skipping matters.

Lecture 4 Part 2: Nonlinear Root Finding, Optimization, and Adjoint Gradient MethodsMIT OpenCourseWare · 18.S096 Matrix Calculus (Edelman and Johnson) · 44 min

The adjoint method for $A(p)x = b$, as in Theorem 3.

• A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation (SIAM 2005): the Bayesian view of inverse problems, with many geophysical examples.
• S. G. Johnson, Notes on adjoint methods for 18.335: the adjoint of $A(p)x = b$ in two pages.
• R.-E. Plessix, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophysical Journal International 167 (2006).
• J. Virieux and S. Operto, An overview of full-waveform inversion in exploration geophysics, Geophysics 74 (2009).
• A. M. Stuart, Inverse problems: a Bayesian perspective, Acta Numerica 19 (2010).
• J. Tromp, C. Tape and Q. Liu, Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, Geophysical Journal International 160 (2005).