01 · Feldman–Hájek: when are two Gaussian priors "the same kind of thing"?¶
Companion to notes.tex, chapter 1. Every figure in the notes is made here (saved to figs/).
The question. We have a reference prior $P$ (our stationary SPDE) and an alternative $Q$ (wells, leakage, a lens...). In finite dimensions they always have a density ratio $dQ/dP$. On a function space, i.e. as the mesh is refined, that ratio may stop existing. Then $P$ and $Q$ are mutually singular: they put their mass on disjoint sets, and one perfectly resolved sample tells you which one generated it.
The one-sentence intuition we will build: two Gaussian measures are equivalent iff you cannot tell them apart from a single infinitely-resolved sample.
Plan
- Finite dimensions: equivalence is the default.
- Infinitely many coordinates: Kakutani and the $\sum t_j^2$ condition.
- Feldman–Hájek in general (in the notes; here we just use it).
- Random functions: why changing the range $\kappa$ is harmless but changing the amplitude $\tau$ is not.
- Back to the aquifer in 1D: wells, leakage, lens.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.sparse import diags
from scipy.sparse.linalg import splu
from pathlib import Path
FIG = Path("figs"); FIG.mkdir(exist_ok=True)
rng = np.random.default_rng(2026)
P_COL, Q_COL, C3, C4 = "#2a78d6", "#eb6834", "#1baf7a", "#eda100" # P = reference (blue), Q = alternative (orange)
plt.rcParams.update({"figure.dpi": 100, "savefig.dpi": 160, "savefig.bbox": "tight",
"axes.spines.top": False, "axes.spines.right": False,
"axes.grid": True, "grid.color": "#e4e3df", "font.size": 10,
"axes.titlesize": 10.5, "legend.frameon": False})
def save(fig, name):
fig.savefig(FIG / f"{name}.pdf"); fig.savefig(FIG / f"{name}.png")
def bayes_error(lP, lQ):
# error of the best one-sample test 'say Q iff log dQ/dP > 0', equal prior odds
return 0.5 * ((lP > 0).mean() + (lQ < 0).mean())
1 · Finite dimensions: equivalence is the default¶
$P=\mathcal N(0,1)$, $Q=\mathcal N(0,1.5^2)$. The density ratio $\frac{dQ}{dP}(x)=\frac{1}{1.5}\exp\!\big(\tfrac{x^2}{2}(1-1.5^{-2})\big)$ is finite and positive for every $x$: every set that is impossible under one is impossible under the other. Equivalent.
The only way to be singular in finite dimensions is for one measure to be degenerate: $Q=\mathcal N(0,\mathrm{diag}(1,0))$ lives on the line $\{x_2=0\}$, which has $P$-probability zero.
x = np.linspace(-5, 5, 801)
p, q = stats.norm.pdf(x, 0, 1), stats.norm.pdf(x, 0, 1.5)
fig, ax = plt.subplots(1, 3, figsize=(12.5, 3.4))
ax[0].plot(x, p, color=P_COL, label=r"$P=\mathcal{N}(0,1)$")
ax[0].plot(x, q, color=Q_COL, label=r"$Q=\mathcal{N}(0,1.5^2)$")
ax[0].set_title("(a) two Gaussians on the line"); ax[0].legend(fontsize=9)
ax[1].plot(x, q / p, color="k"); ax[1].set_yscale("log")
ax[1].set_title(r"(b) $dQ/dP$: finite and $>0$ everywhere $\Rightarrow$ equivalent")
zP = rng.standard_normal((400, 2)); zQ = np.c_[rng.standard_normal(400), np.zeros(400)]
ax[2].scatter(*zP.T, s=6, color=P_COL, alpha=.5, label=r"$P=\mathcal{N}(0,I_2)$")
ax[2].scatter(*zQ.T, s=6, color=Q_COL, label=r"$Q=\mathcal{N}(0,\mathrm{diag}(1,0))$")
ax[2].set_aspect("equal"); ax[2].set_xlim(-4, 4); ax[2].set_ylim(-3.5, 3.5)
ax[2].set_title("(c) singular: Q lives on a line P never hits"); ax[2].legend(fontsize=8, loc="upper left")
save(fig, "fh01_finite_dim"); plt.show()
2 · Infinitely many coordinates: equivalence can break¶
Take $x=(x_1,\dots,x_n)$ with $P=\mathcal N(0,I)$ and $Q=\mathcal N(0,\mathrm{diag}(1+t_j))$. The log likelihood ratio of one sample is a sum of independent terms
$$\ell_n(x)=\log\frac{dQ}{dP}(x)=\sum_{j=1}^n\Big[-\tfrac12\log(1+t_j)+\tfrac12\,\tfrac{t_j}{1+t_j}\,x_j^2\Big].$$
For small $t_j$: the gap between its mean under $Q$ and under $P$ is $\approx\tfrac12\sum t_j^2$ and its standard deviation is $\approx\sqrt{\tfrac12\sum t_j^2}$. So the two histograms separate like $\sqrt{\sum_j t_j^2/2}$: if $\sum t_j^2=\infty$ one sample decides between P and Q with certainty.
First the i.i.d. case $t_j=0.2$ for all $j$ (every coordinate has 20% more variance).
def llr_samples(t, reps, under, rng, chunk=500):
# log dQ/dP for P = N(0, I), Q = N(0, diag(1+t)), on `reps` draws from P or Q
t = np.asarray(t, float); a = t / (1 + t)
const = -0.5 * np.log1p(t).sum()
sd = np.ones_like(t) if under == "P" else np.sqrt(1 + t)
out = np.empty(reps)
for i in range(0, reps, chunk):
k = min(chunk, reps - i)
out[i:i + k] = const + 0.5 * (rng.standard_normal((k, t.size)) * sd) ** 2 @ a
return out
ns = [1, 10, 100, 1000, 10000]
fig, ax = plt.subplots(1, 5, figsize=(15, 2.9))
for a_, n in zip(ax, ns):
t = np.full(n, 0.2)
lP, lQ = llr_samples(t, 4000, "P", rng), llr_samples(t, 4000, "Q", rng)
bins = np.linspace(min(lP.min(), lQ.min()), max(lP.max(), lQ.max()), 60)
a_.hist(lP, bins, color=P_COL, alpha=.6, density=True, label="sample from P")
a_.hist(lQ, bins, color=Q_COL, alpha=.6, density=True, label="sample from Q")
a_.axvline(0, color="k", lw=.8)
a_.set_title(f"n = {n}: test error {bayes_error(lP, lQ):.3f}"); a_.set_yticks([])
a_.set_xlabel(r"$\log\, dQ/dP$ of one sample")
ax[0].legend(fontsize=8, loc="upper left")
save(fig, "fh02_iid_separation"); plt.show()
With $n=1$ the best test is barely better than a coin. With $n=10^4$ coordinates it never errs: from one sample you just compute $\frac1n\sum x_j^2\to 1$ or $1.2$ (law of large numbers). In the limit $n\to\infty$, $P$ and $Q$ live on the disjoint sets $\{\frac1n\sum x_j^2\to1\}$ and $\{\frac1n\sum x_j^2\to1.2\}$: mutually singular, even though each coordinate is a harmless 20% change.
Now let the perturbation decay with $j$. Three choices:
- $t_j=0.2$ (i.i.d.): $\sum t_j^2=\infty$, fast;
- $t_j=j^{-1/4}$: $\sum t_j^2=\sum j^{-1/2}=\infty$, slowly;
- $t_j=1/j$: $\sum t_j^2=\pi^2/6<\infty$.
Kakutani's theorem (1948) says a product of equivalent measures is equivalent iff the product of the per-coordinate Hellinger affinities $\rho_j=\int\sqrt{p_jq_j}=\big(2\sqrt{1+t_j}/(2+t_j)\big)^{1/2}$ stays positive; since $-\log\rho_j\approx t_j^2/16$, that is again exactly $\sum t_j^2<\infty$.
nmax = 10_000; jj = np.arange(1, nmax + 1)
seqs = {r"$t_j=0.2$ (iid)": np.full(nmax, 0.2),
r"$t_j=j^{-1/4}$": jj ** -0.25,
r"$t_j=1/j$": 1.0 / jj}
cols = [Q_COL, C4, C3]
nchk = [10, 100, 1000, 10000]
err = {k: [] for k in seqs}; hists = {}
for k, t in seqs.items():
for n in nchk:
lP, lQ = llr_samples(t[:n], 2000, "P", rng), llr_samples(t[:n], 2000, "Q", rng)
err[k].append(bayes_error(lP, lQ))
hists[k] = (lP, lQ) # keep the n = 10^4 histograms
fig = plt.figure(figsize=(14, 6.4)); gs = fig.add_gridspec(2, 3, hspace=.45)
ax = [fig.add_subplot(gs[0, i]) for i in range(3)]
for (k, t), c in zip(seqs.items(), cols):
ax[0].loglog(jj, np.cumsum(t ** 2), color=c, label=k)
rho = np.exp(np.cumsum(0.5 * np.log(2 * np.sqrt(1 + t) / (2 + t))))
ax[1].semilogx(jj, rho, color=c, label=k)
ax[2].semilogx(nchk, err[k], "o-", color=c, label=k)
ax[0].set_title(r"(a) $\sum_{j\leq n} t_j^2$ (Feldman–Hájek quantity)"); ax[0].set_xlabel("n"); ax[0].legend(fontsize=8)
ax[1].set_title(r"(b) Hellinger affinity $\prod_{j\leq n}\rho_j$ (overlap of P and Q)"); ax[1].set_xlabel("n")
ax[2].set_title("(c) error of best one-sample test"); ax[2].set_xlabel("n"); ax[2].set_ylim(-.02, .52)
ax[2].axhline(.5, color="#8a8984", lw=.8, ls=":")
for i, ((k, (lP, lQ)), c) in enumerate(zip(hists.items(), cols)):
a_ = fig.add_subplot(gs[1, i])
bins = np.linspace(np.percentile(np.r_[lP, lQ], .2), np.percentile(np.r_[lP, lQ], 99.8), 60)
a_.hist(lP, bins, color=P_COL, alpha=.6, density=True, label="from P")
a_.hist(lQ, bins, color=Q_COL, alpha=.6, density=True, label="from Q")
a_.axvline(0, color="k", lw=.8); a_.set_yticks([])
a_.set_title(f"({'def'[i]}) {k}, n = 10⁴: error {bayes_error(lP, lQ):.3f}")
a_.set_xlabel(r"$\log\, dQ/dP$ of one sample")
if i == 0: a_.legend(fontsize=8)
save(fig, "fh03_decay_kakutani"); plt.show()
Read it: for $t_j=1/j$ the overlap (b) and the test error (c) level off — adding more coordinates never makes the two measures distinguishable. That is equivalence. For $t_j=j^{-1/4}$ they keep falling (slowly) to zero: singular.
2b · Means: the Cameron–Martin condition¶
Same idea for a mean shift $Q=\mathcal N(m,I)$: $\ell_n(x)=\sum m_jx_j-\tfrac12\|m\|^2$ is exactly Gaussian, $\mathcal N(\mp\tfrac12\|m\|^2,\|m\|^2)$, so the best test has error $\Phi(-\|m\|/2)$. Equivalent iff $\sum m_j^2<\infty$ — the shift must lie in the Cameron–Martin space of $P$.
fig, ax = plt.subplots(1, 2, figsize=(10, 3.2))
for m, c, lab in [(1 / jj, C3, r"$m_j=1/j$"), (jj ** -0.5, Q_COL, r"$m_j=1/\sqrt{j}$")]:
M = np.cumsum(m ** 2)
ax[0].loglog(jj, M, color=c, label=lab)
ax[1].semilogx(jj, stats.norm.cdf(-np.sqrt(M) / 2), color=c, label=lab)
ax[0].set_title(r"(a) $\|m\|^2$ over the first n coordinates"); ax[0].legend(); ax[0].set_xlabel("n")
ax[1].set_title(r"(b) one-sample test error $\Phi(-\|m\|/2)$"); ax[1].set_xlabel("n"); ax[1].set_ylim(0, .52)
save(fig, "fh04_mean_shift"); plt.show()
4 · Random functions: range $\kappa$ versus amplitude $\tau$¶
Now the coordinates are the sine coefficients of a random function on $[0,1]$ with $u(0)=u(1)=0$: $$(\kappa^2-\partial_x^2)\,u=\tau\,\xi \quad\Longrightarrow\quad u(x)=\sum_j c_j\sqrt2\sin(j\pi x),\qquad c_j\sim\mathcal N\!\Big(0,\ \frac{\tau^2}{(\kappa^2+\lambda_j)^2}\Big),\quad \lambda_j=(j\pi)^2 .$$ This is the 1D version of our reference prior. Both $P$ and $Q$ are diagonal in the sine basis, so everything from §2 applies with $$1+t_j=\frac{\text{var}_Q(c_j)}{\text{var}_P(c_j)} .$$
- Range change $\kappa_0=6\to\kappa_1=2$: $\;t_j=\big(\tfrac{36+\lambda_j}{4+\lambda_j}\big)^2-1\approx\frac{64}{\pi^2j^2}\to0$, $\sum t_j^2<\infty$ → equivalent.
- Amplitude change $\tau_0=1\to\tau_1=1.3$: $\;t_j=0.69$ for every $j$ → singular (it is §2's i.i.d. case).
Draw all three with the same white noise and look.
J = 4096
jJ = np.arange(1, J + 1); lam = (np.pi * jJ) ** 2
svar = lambda tau, kap: tau ** 2 / (kap ** 2 + lam) ** 2
xg = np.linspace(0, 1, 1001)
Phi = np.sqrt(2) * np.sin(np.pi * np.outer(xg, jJ[:1024]))
xi = rng.standard_normal(J)
cases = [("reference τ=1, κ=6", 1.0, 6.0, P_COL),
("range change τ=1, κ=2", 1.0, 2.0, C3),
("amplitude change τ=1.3, κ=6", 1.3, 6.0, Q_COL)]
fig, ax = plt.subplots(1, 4, figsize=(16, 3.4))
for lab, tau, kap, c in cases:
ax[0].plot(xg, Phi @ (xi[:1024] * np.sqrt(svar(tau, kap)[:1024])), color=c, label=lab)
ax[1].loglog(jJ, svar(tau, kap), color=c, label=lab)
ax[0].set_title("(a) one sample each, same white noise"); ax[0].legend(fontsize=7.5, loc="lower left")
ax[1].set_title(r"(b) coefficient variances $\tau^2/(\kappa^2+\lambda_j)^2$"); ax[1].set_xlabel("mode j")
t_kap = svar(1, 2) / svar(1, 6) - 1; t_tau = svar(1.3, 6) / svar(1, 6) - 1
ax[2].semilogx(jJ, t_kap, color=C3, label=r"range: $t_j\to 0$"); ax[2].semilogx(jJ, t_tau, color=Q_COL, label=r"amplitude: $t_j=0.69$")
ax[2].set_yscale("log"); ax[2].set_title(r"(c) whitened ratio $t_j$"); ax[2].set_xlabel("mode j"); ax[2].legend(fontsize=8)
ax[3].loglog(jJ, np.cumsum(t_kap ** 2), color=C3, label="range change")
ax[3].loglog(jJ, np.cumsum(t_tau ** 2), color=Q_COL, label="amplitude change")
ax[3].set_title(r"(d) $\sum_{j\leq n}t_j^2$: converges vs grows $\propto n$"); ax[3].set_xlabel("n modes"); ax[3].legend(fontsize=8)
save(fig, "fh05_range_vs_amplitude"); plt.show()
print(f"sum t_j^2 (range change, all {J} modes): {np.sum(t_kap**2):.2f}")
sum t_j^2 (range change, all 4096 modes): 103.81
The surprise. In (a) the range change ($\kappa=2$, green) looks by far the most different — much bigger, smoother excursions. The amplitude change (orange) looks like the reference with a slightly larger scale. Yet the green one is equivalent and the orange one is singular.
Equivalence is decided by the fine scales (high modes), because there are infinitely many of them. Large-scale differences, however dramatic to the eye, live in finitely many modes and can never be confirmed with certainty from a single sample.
4b · What you can learn from one sample (microergodicity)¶
The flip side: from one sample resolved to $N$ modes, fit $(\tau,\kappa)$ by maximum likelihood. $\tau$ is pinned down by the fine scales and becomes exact as $N\to\infty$. $\kappa$ is not — its Fisher information stays bounded. 300 replicate samples from the reference.
Ns, R = [16, 64, 256, 1024, 4096], 300
kgrid = np.linspace(0, 30, 601)
fits = {}
for N in Ns:
lamN = lam[:N]
c = rng.standard_normal((R, N)) * np.sqrt(svar(1, 6)[:N])
W = (kgrid[:, None] ** 2 + lamN[None, :]) ** 2 # (kappa grid) x (modes)
tau2 = (c ** 2) @ W.T / N # profiled tau^2 for each kappa
prof = -0.5 * N * np.log(tau2) + np.log(kgrid[:, None] ** 2 + lamN).sum(1) # profile log-lik
i = prof.argmax(1)
fits[N] = (np.sqrt(tau2[np.arange(R), i]), kgrid[i])
fig, ax = plt.subplots(1, 3, figsize=(14, 3.5))
ax[0].boxplot([fits[N][0] for N in Ns], positions=range(len(Ns)), widths=.5, showfliers=False)
ax[0].axhline(1, color=P_COL, ls="--", lw=1); ax[0].set_xticks(range(len(Ns)), Ns)
ax[0].set_title(r"(a) $\hat\tau$ from ONE sample: shrinks to the truth"); ax[0].set_xlabel("modes resolved N")
ax[1].boxplot([fits[N][1] for N in Ns], positions=range(len(Ns)), widths=.5, showfliers=False)
ax[1].axhline(6, color=P_COL, ls="--", lw=1); ax[1].set_xticks(range(len(Ns)), Ns)
ax[1].set_title(r"(b) $\hat\kappa$ from ONE sample: never settles"); ax[1].set_xlabel("modes resolved N")
Itau = 2 * jJ / 1.0 ** 2
Ikap = np.cumsum(0.5 * (4 * 6 / (36 + lam)) ** 2)
ax[2].loglog(jJ, Itau, color=Q_COL, label=r"$I_\tau(N)=2N/\tau^2$")
ax[2].loglog(jJ, Ikap, color=C3, label=r"$I_\kappa(N)\to$ const")
ax[2].set_title("(c) Fisher information in one sample"); ax[2].set_xlabel("modes resolved N"); ax[2].legend(fontsize=8)
save(fig, "fh06_microergodic"); plt.show()
for N in Ns:
th, kh = fits[N]
print(f"N={N:5d} sd(tau_hat)={th.std():.4f} sd(kappa_hat)={kh.std():.2f} P(kappa_hat=0)={np.mean(kh == 0):.2f}")
N= 16 sd(tau_hat)=0.2649 sd(kappa_hat)=3.82 P(kappa_hat=0)=0.00 N= 64 sd(tau_hat)=0.0967 sd(kappa_hat)=2.43 P(kappa_hat=0)=0.00 N= 256 sd(tau_hat)=0.0434 sd(kappa_hat)=2.42 P(kappa_hat=0)=0.00 N= 1024 sd(tau_hat)=0.0228 sd(kappa_hat)=2.55 P(kappa_hat=0)=0.00 N= 4096 sd(tau_hat)=0.0105 sd(kappa_hat)=2.52 P(kappa_hat=0)=0.00
The spread of $\hat\tau$ falls like $1/\sqrt N$ (0.26 → 0.01), while the spread of $\hat\kappa$ stalls at about 2.5 no matter how fine the sample: the fine scales say nothing about $\kappa$. That is why $\kappa$ is the "cheap" direction the stationary fit M1 slid down in the main project (to $\hat\kappa\approx0$) when the wells added extra large-scale variance: moving $\kappa$ costs almost no likelihood.
5 · Back to the aquifer, in 1D¶
Finite differences on $N$ interior nodes, $h=1/(N+1)$, for $$-\big(K u'\big)' + \kappa^2(x)\,u = q,\qquad u(0)=u(1)=0 ,$$ written as $A u = q$ with $q=\xi$ white noise ($\xi_i\sim\mathcal N(0,1/h)$). Then $\Sigma=h^{-1}A^{-1}A^{-\top}$ and precision $Q=h\,A^\top A$. Reference: $K\equiv1$, $\kappa^2\equiv36$.
| truth | change | where in the equation |
|---|---|---|
| well | $q=\xi+a\,g$, $a\sim\mathcal N(0,c)$, $c=50$ | source |
| leakage | $\kappa^2(x)=36-30\,e^{-((x-.75)/.06)^2}$ | zeroth order |
| lens | $K(x)=e^{\,b(x)}$, $b\approx1$ on $[0.4,0.6]$ | principal part |
The whitened operator is $T=\Sigma_0^{-1/2}\Sigma_1\Sigma_0^{-\top/2}-I$. With $\Sigma_0^{1/2}=h^{-1/2}A_0^{-1}$: $I+T = BB^\top$, $B=A_0A_1^{-1}$, so the eigenvalues of $T$ are the squared singular values of $A_0A_1^{-1}$ minus one.
def grid1d(N):
h = 1.0 / (N + 1)
return h, h * np.arange(1, N + 1), h * (np.arange(N + 1) + 0.5) # spacing, nodes, faces
def operator(N, K=lambda x: np.ones_like(x), k2=lambda x: np.full_like(x, 36.0)):
# tridiagonal FD matrix for -(K u')' + k2 u with u(0)=u(1)=0; K evaluated at cell faces
h, x, xf = grid1d(N)
Kf = K(xf)
return diags([-Kf[1:-1] / h**2, (Kf[:-1] + Kf[1:]) / h**2 + k2(x), -Kf[1:-1] / h**2], [-1, 0, 1], format="csc")
lens_b = lambda x: 0.5 * (np.tanh((x - 0.4) / 0.01) - np.tanh((x - 0.6) / 0.01))
K_lens = lambda x: np.exp(lens_b(x))
k2_leak = lambda x: 36 - 30 * np.exp(-((x - 0.75) / 0.06) ** 2)
sig_g, c_well = 0.03, 50.0
g_well = lambda x: np.exp(-((x - 0.25) / sig_g) ** 2) / np.sqrt(sig_g * np.sqrt(np.pi / 2)) # ∫ g² = 1
TRUTH = {"well (source)": dict(color=Q_COL), "leakage (zeroth order)": dict(color=C3), "lens (principal part)": dict(color=C4)}
# --- picture of the three truths ---------------------------------------------------------
N = 1024; h, x, _ = grid1d(N)
A0, Al, Ak = operator(N), operator(N, k2=k2_leak), operator(N, K=K_lens)
xi = rng.standard_normal(N) / np.sqrt(h)
u_ref = splu(A0).solve(xi)
u_true = {"well (source)": splu(A0).solve(xi + np.sqrt(c_well) * g_well(x)),
"leakage (zeroth order)": splu(Al).solve(xi),
"lens (principal part)": splu(Ak).solve(xi)}
coef = {"well (source)": (g_well(x), r"well shape $g(x)$"),
"leakage (zeroth order)": (k2_leak(x), r"$\kappa^2(x)$ (reference 36)"),
"lens (principal part)": (K_lens(x), r"$K(x)$ (reference 1)")}
fig, ax = plt.subplots(2, 3, figsize=(14, 5.6), sharex=True)
for i, name in enumerate(TRUTH):
f, lab = coef[name]
ax[0, i].plot(x, f, color=TRUTH[name]["color"]); ax[0, i].set_title(f"{name}: {lab}")
ax[1, i].plot(x, u_ref, color="#8a8984", lw=1.2, label="reference")
ax[1, i].plot(x, u_true[name], color=TRUTH[name]["color"], lw=1.4, label=name.split()[0])
ax[1, i].legend(fontsize=8, loc="lower left"); ax[1, i].set_xlabel("x")
ax[1, 0].set_ylabel("head sample u(x), same noise")
save(fig, "fh07_aquifer_1d_setup"); plt.show()
def T_eigs(N, which):
h, x, _ = grid1d(N)
A0 = operator(N).toarray()
if which.startswith("well"):
t = np.zeros(N); t[0] = c_well * h * np.sum(g_well(x) ** 2) # rank one: c * h g'g ≈ c ∫ g²
return t
A1 = (operator(N, k2=k2_leak) if which.startswith("leak") else operator(N, K=K_lens)).toarray()
s = np.linalg.svd(A0 @ np.linalg.inv(A1), compute_uv=False)
return np.sort(s ** 2 - 1)
Ngrid = [32, 64, 128, 256, 512, 1024, 2048]
res = {name: {"fro": [], "kl": []} for name in TRUTH}
spec = {}
for N in Ngrid:
for name in TRUTH:
t = T_eigs(N, name)
res[name]["fro"].append(np.sqrt(np.sum(t ** 2)))
res[name]["kl"].append(0.5 * np.sum(t - np.log1p(t)))
if N == 1024: spec[name] = t
print(f"{'N':>6}" + "".join(f"{n.split()[0]:>12}" for n in TRUTH) + " <- ||T||_F")
for i, N in enumerate(Ngrid):
print(f"{N:6d}" + "".join(f"{res[n]['fro'][i]:12.3f}" for n in TRUTH))
N well leakage lens <- ||T||_F
32 50.000 0.638 2.340
64 50.000 0.636 3.180
128 50.000 0.635 4.461
256 50.000 0.635 6.252
512 50.000 0.635 8.796
1024 50.000 0.635 12.407
2048 50.000 0.635 17.522
fig, ax = plt.subplots(1, 3, figsize=(14.5, 3.7))
for name, st in TRUTH.items():
ax[0].loglog(Ngrid, res[name]["fro"], "o-", color=st["color"], label=name)
ax[1].loglog(Ngrid, res[name]["kl"], "o-", color=st["color"], label=name)
ref = np.array(Ngrid, float)
ax[0].loglog(ref, res["lens (principal part)"]["fro"][0] * np.sqrt(ref / ref[0]), ":", color="#8a8984", label=r"$\propto\sqrt{N}$")
ax[1].loglog(ref, res["lens (principal part)"]["kl"][0] * ref / ref[0], ":", color="#8a8984", label=r"$\propto N$")
ax[0].set_title(r"(a) $\|T\|_F$ under refinement"); ax[0].set_xlabel("nodes N"); ax[0].legend(fontsize=8)
ax[1].set_title(r"(b) $\mathrm{KL}(Q\|P)=\frac{1}{2}\sum(t-\log(1+t))$"); ax[1].set_xlabel("nodes N")
for name in ["leakage (zeroth order)", "lens (principal part)"]:
t = spec[name]; ax[2].plot(np.arange(t.size) / t.size, t, color=TRUTH[name]["color"], label=name)
ax[2].axhline(np.exp(-2) - 1, color="#8a8984", ls=":", lw=1, label=r"$K^{-2}-1=e^{-2}-1$")
ax[2].set_title("(c) all eigenvalues of T at N = 1024 (sorted)"); ax[2].set_xlabel("fraction of modes"); ax[2].legend(fontsize=8)
save(fig, "fh08_aquifer_refinement"); plt.show()
The 1D version of the main project's F1/F2: the well and the leakage bump give flat $\|T\|_F$ (equivalent); the lens grows like $\sqrt N$ (in 2D, $N=n^2$, so $\propto n$). Panel (c) shows why: about 20% of the eigenvalues of $T$ (the share of the domain the lens occupies) sit at the symbol value $K^{-2}-1=e^{-2}-1\approx-0.86$ — and that share is a fixed fraction of all modes, so their number grows with $N$.
5b · The one-sample test, for the aquifer¶
Exactly the experiment of §2, but now with the real operators: $\log\frac{d\mu_1}{d\mu_0}(u)=\log|A_1|-\log|A_0|-\tfrac h2\big(\|A_1u\|^2-\|A_0u\|^2\big)$ (and a Woodbury form for the well).
def logdet(A):
return np.log(np.abs(splu(A).U.diagonal())).sum()
def llr_aquifer(N, name, R, rng):
h, x, _ = grid1d(N)
A0 = operator(N); lu0 = splu(A0)
xi0 = rng.standard_normal((N, R)) / np.sqrt(h); xi1 = rng.standard_normal((N, R)) / np.sqrt(h)
if name.startswith("well"):
g = g_well(x); s = h * g @ g
u0 = lu0.solve(xi0)
u1 = lu0.solve(xi1 + np.sqrt(c_well) * g[:, None] * rng.standard_normal(R))
ll = lambda u: -0.5 * np.log1p(c_well * s) + 0.5 * c_well * (h * g @ (A0 @ u)) ** 2 / (1 + c_well * s)
return ll(u0), ll(u1)
A1 = operator(N, k2=k2_leak) if name.startswith("leak") else operator(N, K=K_lens)
u0, u1 = lu0.solve(xi0), splu(A1).solve(xi1)
dld = logdet(A1) - logdet(A0)
ll = lambda u: dld - 0.5 * h * (np.sum((A1 @ u) ** 2, 0) - np.sum((A0 @ u) ** 2, 0))
return ll(u0), ll(u1)
Nshow = [64, 512, 4096]
fig, ax = plt.subplots(3, 3, figsize=(13, 7.6))
for r, name in enumerate(TRUTH):
for cidx, N in enumerate(Nshow):
lP, lQ = llr_aquifer(N, name, 1500, rng)
a_ = ax[r, cidx]
hi = 97 if name.startswith("well") else 99.7 # the well's LLR has a long right tail
bins = np.linspace(np.percentile(np.r_[lP, lQ], .3), np.percentile(np.r_[lP, lQ], hi), 60)
a_.hist(lP, bins, color=P_COL, alpha=.6, density=True, label="from reference")
a_.hist(lQ, bins, color=TRUTH[name]["color"], alpha=.6, density=True, label="from truth")
a_.axvline(0, color="k", lw=.8); a_.set_yticks([])
a_.set_title(f"{name}, N = {N}: error {bayes_error(lP, lQ):.3f}", fontsize=9.5)
ax[r, 0].legend(fontsize=8, loc="upper left")
for a_ in ax[-1]: a_.set_xlabel(r"$\log\, d\mu_1/d\mu_0$ of one head sample")
fig.tight_layout()
save(fig, "fh09_aquifer_one_sample_test"); plt.show()
Same story as §2: for the well and the leakage bump the test error stays put as the mesh is refined (equivalent: a finer mesh gives no new evidence). For the lens it drops toward zero: a finely resolved head field carries the lens's fingerprint in every small-scale wiggle (singular).
Take-aways¶
- Equivalence = "same null sets" = a density $dQ/dP$ exists = one perfect sample cannot tell them apart.
- Diagonal case: equivalent iff $\sum_j t_j^2<\infty$ (variances) and $\sum_j m_j^2<\infty$ (means). Feldman–Hájek is this, written basis-free: $T$ Hilbert–Schmidt + mean in Cameron–Martin.
- It is decided at fine scales. Big, eye-catching large-scale changes (range $\kappa$, wells, leakage) are harmless; uniform small-scale changes (amplitude $\tau$, conductivity $K$) are fatal.
- Practical meaning: what a fine sample determines exactly ($\tau$, $K$) you cannot "correct" with finitely many directions; what it never determines ($\kappa$) is a cheap knob a fit can slide along.
Check yourself¶
- Why is $t_j=1/\sqrt{j}$ a borderline case? (Compute $\sum t_j^2$.)
- In 2D the range change gives $t\sim 1/\lambda$ and there are $\sim\Lambda$ modes with $\lambda\le\Lambda$. Is $\sum t^2$ finite? What about 4D?
- A constant conductivity $K\equiv 2$: equivalent or singular? (Hint: factor it out of the operator.)