Matched filter — derivation, covariance flavours, and the link to 3D-Var

Matched filter for hyperspectral methane retrieval¶

The matched filter (MF) is the maximum-SNR linear detector for a known target signature in additive Gaussian noise. In hyperspectral methane retrieval it is the workhorse baseline — fast, closed-form, and surprisingly competitive with full variational retrievals when the signal is weak and the scene is statistically homogeneous. This note derives the MF from scratch, lays out the four covariance flavours implemented in plumax.matched_filter, and shows the exact sense in which the MF is the one-step limit of the dual / PSAS 3D-Var solver from the companion 00_variational_derivation.md.

The derivation borrows notation from Funk et al. ^[1], Theiler & Foy ^[2], and Thompson et al. ^[3]; the operator-based implementation mirrors the structured-linear-algebra view of Golub & Van Loan ^[4] and the Woodbury-friendly design in gaussx.

1. The signal model¶

We observe a per-pixel radiance spectrum $x \in \mathbb{R}^B$ (one vector per pixel of the scene, with $B$ spectral bands). The model is

where

• $\mu \in \mathbb{R}^B$ is the background mean spectrum (what you would observe at a pixel free of plume enhancement),

• $t \in \mathbb{R}^B$ is the target signature (the spectral shape that a unit VMR enhancement would leave in $x$),

• $\alpha \ge 0$ is the plume amplitude (the unknown — pixels with $\alpha = 0$ are background, $\alpha > 0$ is a plume),

• $\Sigma \in \mathbb{R}^{B \times B}$ is the background clutter covariance (noise + surface heterogeneity + calibration residual), assumed Gaussian.

Under this model, the log-likelihood of a pixel is, up to additive constants,

Setting $\partial \ell / \partial \alpha = 0$ gives the maximum-likelihood estimator of the amplitude:

This is the matched filter. It is unbiased ($\mathbb{E}[\hat\alpha] = \alpha$), has the minimum variance among all unbiased linear estimators (Gauss–Markov), and maximises SNR in the Neyman–Pearson sense for detecting the alternative $\alpha > 0$ against the null $\alpha = 0$.

1.1 The SNR and detection threshold¶

Under the null ($\alpha = 0$) the score $\hat\alpha(x)$ is $\mathcal{N}(0, 1/(t^\top \Sigma^{-1} t))$. Under an alternative with amplitude $\alpha$ it is $\mathcal{N}(\alpha, 1/(t^\top \Sigma^{-1} t))$. So the detection SNR is

and the Gaussian detection threshold for a false-alarm rate (FAR) is

Both are implemented in matched_filter.core as matched_filter_snr and detection_threshold.

2. Practical decomposition: precompute once, apply per pixel¶

Equation (5) is deceptively cheap. Per pixel, it’s a dot product of length $B$ — except $\Sigma^{-1}$ appears. Done naively, you invert a $B \times B$ matrix once, but that discards the chance to exploit any structure in $\Sigma$.

The right decomposition is

Then per pixel

— one scalar solve and one scalar, followed by one dot product per pixel (a BLAS call when vectorised over the image). This is exactly what apply_image in core.py does:

def apply_image(cube, mean, cov_op, target):
    w = gx.solve(cov_op, target)      # ← Σ⁻¹ t, one solve
    z = jnp.dot(target, w)             # ← scalar
    residual = cube - mean             # ← broadcast
    scores = einops.einsum(residual, w, "H W B, B -> H W")
    return scores / z

The key move is gx.solve(cov_op, target). Because cov_op is an abstract lineax.AbstractLinearOperator, the solve dispatches on structure: a dense operator triggers a Cholesky, a diagonal operator triggers a per-element divide, and a low-rank update λI + V D V^\top triggers Woodbury — with no change to the MF kernel above.

3. Four flavours of $\Sigma$¶

The four reference files in jej_vc_snippets/methane_retrieval/ originally wrote out four separate implementations of equation (5), one per covariance model. In the operator-based view they collapse into four estimation routines that all return the same type — a lineax.AbstractLinearOperator.

The Ledoit–Wolf and OAS estimators come from sklearn.covariance; the randomised low-rank SVD comes from sklearn.decomposition.TruncatedSVD. Because MF estimation is off the JAX gradient path — the background is treated as a fixed input — we lean on scikit-learn freely. Only the per-pixel application stays in JAX, vectorised with jax.vmap / einsum for GPU-friendly batching.

3.1 Why low-rank + Tikhonov matters¶

For $n_{\text{samples}} \gg B$ and a homogeneous scene, the spectral covariance is approximately low-rank: a handful of surface-reflectance PCs explain most of the variance. Writing $\Sigma = \tau I + V D V^\top$ with $V \in \mathbb{R}^{B \times k}$ and $D \in \mathbb{R}^{k \times k}$ diagonal lets us apply Woodbury:

The inner inverse is $k \times k$, so $w = \Sigma^{-1} t$ costs $\mathcal{O}(B k + k^3)$ instead of $\mathcal{O}(B^3)$. For a typical hyperspectral methane setup ($B = 400$, $k = 16$), this is a ~1000× speedup on the MF precompute — visible to the user only in that estimate_cov_lowrank returns a gaussx.LowRankUpdate instead of a dense operator.

4. Relationship to 3D-Var¶

3D-Var minimises (see 00_variational_derivation.md, eq. 2):

where $d = y - H(x_b)$ is the innovation and $H'$ is the tangent-linear of the forward operator. The PSAS / dual path solves in observation space:

The MF is exactly the PSAS analysis with a flat (infinite-variance) prior and identity state-space coupling.

Take $B = \sigma_B^2 I$ with $\sigma_B \to \infty$, $R = \Sigma$, and $H' \delta x = \alpha t$ (the tangent-linear sees a unit plume amplitude as the target signature). Then equation (13) becomes

As $\sigma_B \to \infty$ the second term saturates, and a short calculation gives

So running a single CG iteration of run_dual_psas with the flat prior and linearised $H'$ gives the MF score as a by-product. The MF is therefore not a heuristic: it is the exact, closed-form, flat-prior dual 3D-Var solution, and the SNR / threshold formulas follow without further work.

The corollary is that when $B$ is not flat — i.e. when we have real prior information about the VMR field (spatial smoothness, a climatological mean, etc.) — the full 3D-Var retrieval strictly dominates the MF. See notebook 01_mf_vs_3dvar.ipynb for a numerical comparison on the same scene.

5. Where the target signature $t$ comes from¶

For methane, $t(\nu)$ is the spectral signature of a unit VMR enhancement. From the Beer–Lambert forward model $L(\nu) = L_0(\nu)\,\exp(-a(\nu)\,\Delta\mathrm{VMR})$, the Taylor expansion at the background state gives

Rather than hard-coding this, we take the JVP of the existing observation operator:

t = linear_target_from_obs(obs, vmr_background, pattern="uniform")
#   internally: jax.jvp(obs.forward, (vmr_background,), (ones,))

This guarantees the target is exactly consistent with whatever forward model (PSF, GSD, SRF) the retrieval uses — no drift between “target generator” and “forward model” implementations. For strong plumes where Beer–Lambert saturation matters, nonlinear_target_from_obs(obs, x_b, amplitude=α) computes the exact finite-amplitude response $L(x_b + \alpha\,\delta x) - L(x_b)$ instead.

6. Summary¶

The pieces of plumax.matched_filter:

The two companion notebooks:

• 01_mf_vs_3dvar.ipynb — MF versus full 3D-Var on the same HAPI methane scene used in the assimilation tutorial, so you can see where the flat-prior MF gives up information.

• 02_covariance_flavours.ipynb — dense vs shrunk vs low-rank vs GMM-clustered $\Sigma$ on a heterogeneous scene, comparing detection power and runtime.