Gradient verification
Every sensitivity a recipe exposes is verified three ways, and the three catch different things.
The Taylor test checks that the gradient is the derivative of the objective the forward actually computes. With \(\Phi(m) = \tfrac12\|F(m) - d\|^2\),
\[ \Phi(m_0 + h\,\delta m) - \Phi(m_0) - h\,\langle \nabla\Phi(m_0), \delta m\rangle \;=\; \mathcal{O}(h^2), \]
so the second-order error falls with a fitted log-log slope of 2. The assertion is one-sided (slope > 1.7): converging faster than \(h^2\) means the second-order term nearly vanishes for that perturbation, which is not an error, so an upper bound would fail correct gradients.
The adjoint dot product checks the propagator on its own, independently of any sensitivity: \(\langle Fs, r\rangle = \langle s, F^{\!\top} r\rangle\). A recipe whose adjoint is only the forward stencils run backward in time — a pseudo-adjoint — fails this even when its gradients look plausible, and the Taylor slopes alone cannot tell the two apart. It is the test to fix first when both are failing: no sensitivity can recover from an adjoint that is not the transpose.
The transpose test checks one kernel against one operator, with no wave equation and no adjoint in between: it takes the operator the propagation applies, differentiates it numerically with respect to the parameter at a few points, and compares. It is the instrument that separates a transposition mistake from a wavefield one, and it is sharp where the other two are not — every ISIC kernel defect below was 145, 80 and 58 per cent out pointwise here while its Taylor slope passed.
Regenerate every figure on this page with
python scripts/plot_gradients.pywhich runs the same checks as tests/test_gradient.py on a reduced model (N=101, nt=1000) that reproduces every slope to about 0.03.
Adjoint dot product
The mismatch is measured against \(\|Fs\|\,\|r\|\), the scale of the bilinear form, rather than against \(\langle Fs, r\rangle\) itself. With a random \(r\) that inner product is a largely cancelling sum whose value swings by an order of magnitude with the number of time steps, so dividing by it reports the conditioning of one draw rather than the accuracy of the transpose: the same exact adjoint reads 1.9e-5, 1.2e-4 and 1.1e-5 at 200, 400 and 800 steps.
Every recipe is exact without a free surface. With one, only the acoustic clears the bar — see Known limitations.
Acoustic and anisotropic-acoustic
Imaging-condition parameters
Z and mZ are the reflection-FWI pair, and neither is a stored field, so neither is perturbed directly. \(Z = \rho c\) at fixed m is a pure density move, which is why its direction is b. mZ is squared slowness at fixed impedance, so it moves both fields at once – \(\delta\ln c = -\tfrac12\delta\ln m\) and \(\delta\ln\rho = +\tfrac12\delta\ln m\) leave \(Z\) alone.
Both instruments, because neither subsumes the other. The panels above show the Taylor errors; the ratio \(\mathrm{d}\Phi/\mathrm{d}h \,/\, \langle g,\delta m\rangle\) by central difference is printed alongside. For Z the slope is the sharper of the two: scaling a correct gradient by 1.1 takes it from 2.045 to 0.492 and flipping its sign takes it to 1.020. For mZ it is not. mZ is a transmission parameter, weakly seen by a surface acquisition, so \(h\langle g,\delta m\rangle\) stays small against the curvature and the same three cases read 1.910 exact, 1.995 ten per cent out, and 1.395 sign-flipped — a wrong gradient scoring better than the right one. Widening the background offset does not rescue it (1.575 against 1.704 at an offset of 0.10). The ratio separates all three, so both are asserted.
The fit window matters as much as the bar: \(h\in[0.015,\,0.25]\). Above it the quadratic term dominates \(\Phi(h)-\Phi(0)\) and below it error2 is at the float32 floor of a misfit summed over the whole record, and a fit spanning either reports neither — over eight decades of \(h\) the same exact Z gradient fits 2.43 and the same exact mZ gradient fits 1.78.
Both are built from the same two correlations the buoyancy gradient is built from, a mass one \(K_m\) and a stiffness one \(K_s\):
\[ g_m = -\mathrm{d}t\, b\, K_m, \quad g_b = -\mathrm{d}t\,(m K_m + K_s), \quad g_{mZ} = -\mathrm{d}t\,(m K_m - K_s), \]
so Z is the buoyancy gradient and mZ is \(2(m/b)g_m - g_b\) pointwise. That identity is the sharpest check available on the condition, and test_isic_slowness_composition holds it to 1e-6 – worth having, because mZ differences two kernels that can cancel to one part in 16, which amplifies a finite difference’s own error by the same factor.
What was wrong
Every entry in the first column was a kernel that passed its Taylor test.
| before | after | |
|---|---|---|
iso-acoustic Z |
0.906 | 1.0000 |
iso-acoustic mZ |
1.081 | 0.9996 |
tti-zhang Z |
0.980 | 1.0002 |
tti-zhang mZ |
-0.920 | 0.9992 |
tti-fletcher Z |
1.027 | 0.9995 |
tti-fletcher mZ |
0.80 / 1.26 | 0.9924 |
The right-hand column no longer reproduces on this branch: the same test now reads 0.9357, 0.9363 and 0.9354 for the three Z entries and 0.9247 (tti-zhang) and 1.1421 (tti-fletcher) for mZ, so five of the six fail their 5 per cent bar. It is not the kernels – each still passes the transpose test below – and it is not the imaging condition either: the same 0.935 turns up on vp, and on the Born/adjoint dot product of every recipe (see Born). Something scales the accumulated gradient by about 1.07, and it arrived after the column was measured.
Three separate defects, each verified in isolation by a transpose test that holds no wave equation – just the operator the propagation applies, and a central difference of it with respect to b:
- The stiffness kernel was taken at the node.
grad(u).grad(ua)is \(\mathrm{d}/\mathrm{d}b\) of \(\nabla\!\cdot\!(b\nabla u)\) in the continuum but not for the discrete objective:bis sampled as a two-point average onto the half cells, so the transpose is the staggered derivative product averaged back. The node-centred form is 145 per cent out pointwise. - Zhang’s rotated part dropped the cross terms (80 per cent). \(G_{zz}\) pairs each dimension’s derivative of
uwith the whole rotated derivative ofua, not just the same dimension’s. Zeroing the tilt took the assembled gradient from 0.993 to 0.9995, which is what identified it. - Fletcher expanded the transpose by the product rule (58 per cent). \(d_a[c(d_c q)w] - c(d_ad_cq)w = c(d_cq)(d_aw) + (d_ac)(d_cq)w\) holds in the continuum, but \(d_ad_cq\) is the compact second-derivative kernel and not \(d_a\) applied to \(d_cq\). Written as it stands the kernel is exact.
And one that is not a kernel at all: the saved field did not carry the points its own transpose reads. A stencil of the form \(d_a(d_cq\cdot w)\) reaches space_order points into the saved wavefield where the propagation reached space_order//2, and Fletcher allocates the narrow halo. devito neither raises nor falls back – the buoyancy kernel came back at half its value with the forward run bit-identical. test_saved_field_halo_covers_ sensitivity compares what the kernel reads against what the field allocates, symbolically, rather than waiting for a gradient test to notice.
The background has to be smooth
The mZ sweep starts from a background offset from the truth – from the truth itself the residual is identically zero and F0 is 0.0 to the last bit. That offset was cell-to-cell random, and on a background that oscillates every cell the Fletcher m gradient reads 0.75 against a central difference where it reads 1.000 on a smooth one. m on Fletcher is the pseudo-adjoint’s gradient and not the strict Fréchet derivative, which its own docstring says; a background that rough is not a model, and measuring an imaging condition on one says nothing about it. Smoothing the offset moved mZ from 0.80 to 0.99 with no change to any kernel.
Born
Both are linearized too, on all three recipes. The Born sources are transposed from the accumulators rather than derived again – the identities above pin them:
\[ J_Z = J_b, \qquad J_{mZ}\,q = J_m\!\left(2\tfrac{m}{b}q\right) - J_b(q), \]
and since \(J_m(x) = x\,L(u)/m\) the second collapses to the same two terms as the first with the stiffness one added instead of subtracted. Against the adjoint, \(\langle J q, d\rangle / \langle q, J^{\!\top} d\rangle\) reads 1.000005 (Z) and 1.000028 (mZ); as a Fréchet derivative the second-order error falls at 4.075 and 3.824.
mZ needs a larger perturbation than the other parameters. Its two legs cancel to about a half, so the linearized signal sits closer to the float32 floor of differencing two forwards: at 5 per cent error2 stops falling after two steps (3.81, 1.01, 0.41, 0.33) and the fitted rate reads 2.5 on a Born source that is exact.
What bounds it from above is not the amplitude of q but the exponent h*q/b, since mZ scales the velocity by exp(-h*q/b): scaled like Z that reaches 1.36, a factor 3.9 on the velocity, which the acoustic forward absorbs and Fletcher’s does not – it goes unstable and the sweep comes back NaN. The sweep bounds it at 0.4 instead and reads the best ratio rather than the mean, test_jacobian_forward’s convention and for its reason.
On the anisotropic recipes
The same two identities, with the stiffness written per physics. Both TTI recipes carry the buoyancy inside their stiffness operators – Gxx/Gzz for Zhang, the rotated H1/H2 for Fletcher – so the mass term is only half of d/db. Writing \(H^{(q)}\) for the stiffness evaluated with q in place of b, which is its derivative in that direction since every one of them is linear in the buoyancy,
\[ J_b\,q = (q/b)\,H - H^{(q)}, \]
the first term being Zhang’s mass term rewritten with the PDE identity and Fletcher’s vp**2/b prefactors, which its stencil carries instead of a mass term. q is passed to the operators as their b argument – they take one, and each is linear in it – rather than substituted afterwards, which would match nothing: the buoyancy enters those stencils already sampled at the half cells.
The Zhang buoyancy source used to be the mass term alone, dp*H0/b, with no counterpart to the stiffness term – the buoyancy Born source of a different PDE, one whose stiffness holds no b. Nothing caught it because no Born test covered a density parameter outside iso-acoustic; test_jacobian_forward now does. Dropping the stiffness leg is what that catches: the second-order error then falls as h and no ratio in the six imaging-condition sweeps reaches 2.2, against 3.5 and up with it.
Measured: as Fréchet derivatives, b converges at 4.10 (Zhang) and 4.22 (Fletcher), Z at 4.10 to 4.22 and mZ at 3.46 to 3.91 across the three recipes. Against the adjoint, \(\langle Jq, d\rangle/\langle q,
J^{\!\top}d\rangle\) reads 0.930 to 0.938 – for every recipe and every parameter, vp included.
That last offset is on the gradient side of the pair, not the linearization: the impedance gradient reads the same 0.935 against a plain central difference of the objective, with no Born anywhere in it, on all three recipes. Both measurements say the accumulated gradient is currently about 7 per cent too large; the Born sources are what converge at 4.
That last check has to run on a smooth background. Fletcher’s gradient is the pseudo-adjoint’s rather than the strict Fréchet derivative, so it is only the transpose where the two agree: over three layers its vp ratio falls to 0.55 and its b ratio to -0.84 – a cancelling inner product against a gradient that is 45 per cent out pointwise can land anywhere, sign included – where the same Born sources read 0.93 on one layer. It is the recipe’s gradient and not the linearization that the layers break, and vp, whose pair predates the density ones, measures the same both ways.
Which way the scale separation goes
Z sums the mass and stiffness kernels and mZ differences them, and the two annihilate opposite scattering angles: on a homogeneous model with one source-receiver pair, where the only event is the transmitted wave, Z retains 11 per cent of the cross-correlation kernel’s energy away from the two singular points and mZ retains 80. So Z is the reflectivity condition and mZ the transmission one, which is what the parameterization says they are and what JUDI’s ic="isic" / ic="fwi" pair does with the same signs. examples/imaging_conditions.ipynb has the picture.
Worth recording how nearly this was called the other way round. The test that covered it compared the spectral centroid of two whole images, and with the kernels corrected it read 1.08 where it wanted 1.15 — a fail on a condition that removes four fifths of the artifact it is supposed to remove. Smoothing the images first, to get past Fletcher’s grid-scale ringing, sent it to 0.86 and made it look like the ordering had reversed. It had not: a centroid over a whole image is not sensitive to the thing being claimed, and the fix was to measure the artifact directly rather than to move the bar.
Elastic
vp and vs do not appear: an elastic model stores the Lamé parameters, the stiffness is assembled from them, and the propagation never reads a velocity, so perturbing model.vp leaves it untouched and there is nothing to Taylor test against. They reach the stiffness through \(\lambda = v_p^2/\rho - 2v_s^2\) and \(\mu = v_s^2/\rho\), and test_elastic_velocity_chain_rule checks them against the lam/mu gradients, which are tested here.
Known limitations
| Case | Symptom |
|---|---|
Any recipe but iso-acoustic, fs=True |
the propagator is not the transpose — 6e-5 (tti-selfadjoint) to 2e-3 (iso-elastic) |
Any elastic gradient, fs=True |
slope 1.17–1.30, downstream of the above |
tti-fletcher epsilon/delta, fs=True |
slope 1.50 / 1.48 |
tti-elastic epsilon/delta |
slope 1.50 / 1.42 — exact on vti-elastic, so it is the rotation |
tti-fletcher b/Z images |
grid-scale ringing from the expanded operator form — the gradient is exact (0.9995), the image is not smooth |
Z/mZ Born on the TTI recipes |
tti-fletcher has no b Born source and tti-zhang’s carries only the mass term |
These are xfailed in tests/test_gradient.py rather than silently loosened; the bar stays at 2 for everything else.
Two things that were not the cause, both ruled out by measurement rather than by reading the code, and both worth not re-deriving:
- The free-surface mirror is not the problem. It is part of the modeling and is self-adjoint: the image folds the stencil into one that is symmetric on the interior, so applying the same mirror in both directions is the transpose.
iso-acousticdemonstrates it, passing the dot product atfs=True. Replacing the mirror with an accumulating form makes it worse. - The acquisition depth is a probe, not a verdict. Moving the source and receivers 30 cells down takes every recipe back to exact, but at 1.25 cells the interpolation stencil never reaches the halo — depth moves the acquisition out of the band near the surface where the adjoint is inaccurate, rather than switching a component on or off.