Paper deep dive
Two-Sample Testing via Path-based Inference
Eshant English, Wei-Cheng Lai, Yanfeng Yang, Kenji Fukumizu, Taiji Suzuki, Christoph Lippert
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Modern deep generative models are primarily studied for their ability to generate realistic samples, yet the generative dynamics they learn can also serve as objects of statistical inference. We develop this idea for two-sample testing, the problem of deciding whether the same distribution generated two finite datasets. Using stochastic interpolants, we connect both distributions to a shared Gaussian bottleneck, so that each half of the resulting path is a Gaussian channel acting on a single population. We prove that the null hypothesis holds if and only if the population denoiser, or equivalently, the velocity fields of the two halves, coincide at any single noise level, which amounts to a reflection symmetry of the path about the bottleneck. Deviations from this symmetry yield a continuum of two-sample witnesses, which we estimate via held-out regression risks on learned denoisers and velocities and aggregate along the path; under an information-theoretic weighting, the aggregated discrepancy equals the Jeffreys divergence between the noise-smoothed distributions. Calibrating the resulting statistics by permutation yields tests that are valid in finite samples for any trained networks and consistent when the fields are learned accurately. On a synthetic benchmark and three image benchmarks, the proposed tests improve power over the strongest baseline by up to 33 percentage points at an equal total sample budget, with the best choice of regression representation and path weighting depending on the data modality. These results show that generative paths provide a principled representation for statistical testing, extending stochastic-interpolant models beyond generation.
Tags
Links
- Source: https://arxiv.org/abs/2610.05684v1
- Canonical: https://arxiv.org/abs/2610.05684v1
Trouble viewing inline? Open PDF directly →
Full Text
162,479 characters extracted from source content.
Expand or collapse full text
Two-Sample Testing via Path-based Inference Eshant English †thanks: Equal contribution, sorted in alphabetical order. Contact email: eshant.english@outlook.co.nz, wei-cheng.lai@hpi.de Affiliation: The University of Tokyo, Tokyo, Japan Affiliation: Hasso Plattner Institute for Digital Engineering, Germany Wei-Cheng Lai Affiliation: Hasso Plattner Institute for Digital Engineering, Germany Yanfeng Yang Affiliation: Institute of Statistical Mathematics, Japan Affiliation: The Graduate University for Advanced Studies, SOKENDAI, Japan Kenji Fukumizu Affiliation: Institute of Statistical Mathematics, Japan Taiji Suzuki Affiliation: The University of Tokyo, Tokyo, Japan Affiliation: RIKEN Centre for Advanced Intelligence Project (AIP), Tokyo, Japan Christoph Lippert Affiliation: Hasso Plattner Institute for Digital Engineering, Germany Affiliation: Hasso Plattner Institute for Digital Health at the Icahn School of Medicine at Mount Sinai, USA Abstract Modern deep generative models are primarily studied for their ability to generate realistic samples, yet the generative dynamics they learn can also serve as objects of statistical inference. We develop this idea for two-sample testing, the problem of deciding whether the same distribution generated two finite datasets. Using stochastic interpolants, we connect both distributions to a shared Gaussian bottleneck, so that each half of the resulting path is a Gaussian channel acting on a single population. We prove that the null hypothesis holds if and only if the population denoiser, or equivalently, the velocity fields of the two halves, coincide at any single noise level, which amounts to a reflection symmetry of the path about the bottleneck. Deviations from this symmetry yield a continuum of two-sample witnesses, which we estimate via held-out regression risks on learned denoisers and velocities and aggregate along the path; under an information-theoretic weighting, the aggregated discrepancy equals the Jeffreys divergence between the noise-smoothed distributions. Calibrating the resulting statistics by permutation yields tests that are valid in finite samples for any trained networks and consistent when the fields are learned accurately. On a synthetic benchmark and three image benchmarks, the proposed tests improve power over the strongest baseline by up to 33 percentage points at an equal total sample budget, with the best choice of regression representation and path weighting depending on the data modality. These results show that generative paths provide a principled representation for statistical testing, extending stochastic-interpolant models beyond generation. †shortheadings: Two-Sample Testing via Path-based Inference / Two-Sample Testing via Path-based Inference†firstpage: 1 keywords Two-Sample Testing, Generative Modelling, Dynamic Transport Contents 1 Introduction 2 Preliminaries 2.1 Stochastic Interpolants 2.2 Two-Sample Testing 3 Path-Based Two-Sample Testing 3.1 A Gaussian-Bottleneck Interpolant 3.2 Denoiser and Velocity Fields 3.3 Path Invariance Characterises the Null 3.4 Held-Out Risk Witnesses 3.5 Path Aggregation 3.6 Held-Out Permutation Calibration 3.7 Consistency 3.8 Relationship between the Four Variants 3.9 Monte Carlo Approximation over the Path 4 Experiments 4.1 High-Dimensional Gaussian Mixtures 4.2 MNIST 4.3 OrganAMNIST 4.4 Fashion-MNIST 4.5 Type-I Error 4.6 Samples from the Fitted Models 4.7 Ablation on Model Architecture 5 Observation-Level Asymmetry under Corruption 6 Related Work 7 Discussion and Limitations A Experimental Details and Notations B Per-Benchmark Power and Type-I Error Curves References 1 Introduction Dynamical-transport generative models [20, 10], including diffusion models [54], flow matching [37] and stochastic interpolants [1, 2], are the state of the art for data synthesis. The underlying idea is simple: given two distributions ρ0 _0 and ρ1 _1, one constructs a continuous dynamical process (xt)t∈[0,1](x_t)_t∈[0,1] such that xt∼ρ0x_t _0 at t=0t=0 and xt∼ρ1x_t _1 at t=1t=1. The choice of t∈[0,1]t∈[0,1] is purely a unifying convention; diffusion models, for instance, instead run a fixed noising process over an extended or infinite time horizon to reach a Gaussian. Conventions differ across implementations, but in each case, a deep network is trained to approximate the field that generates this process: a velocity field b(t,x)b(t,x) in the deterministic case, or a drift and score in the stochastic case. Combining this simple modelling principle with the expressivity of deep networks has yielded strong empirical results and enabled applications in high-fidelity computer vision [50], audio and speech synthesis [30], and computational biology, such as protein design [59]. Despite this abundance of generative applications, the direct use of these models for statistical inference remains largely unexplored. The continuous transformation carrying a sample from ρ0 _0 to ρ1 _1 traces out a trajectory, which we call a path, describing how a sample evolves from ρ0 _0 into a sample from ρ1 _1. This path may follow stochastic or deterministic dynamics; we restrict our attention to deterministic dynamics, under which each sample follows a unique, identifiable path. We argue that the path itself, not merely its endpoints, is a rich statistical object in its own right. At the population level, suppose we want to test whether the distribution of collected data has changed over time. We can construct a dynamical transport model that carries samples from ρ0 _0 to ρ1 _1. If the dynamics realised optimal transport, then under ρ0=ρ1 _0= _1 the optimal map would be the identity so that every sample would remain at rest, and any movement along the path would be evidence of a distributional change. Learned dynamics are rarely optimal, however, so this picture serves only as intuition. In Section 3, we show that an analogous property holds exactly for an interpolant that routes both distributions through a shared Gaussian bottleneck [1]: the null hypothesis holds if and only if the path is reflection symmetric about the bottleneck (Propositions 4 and 1). We refer to this broader line of reasoning as path-based inference. Paths are particularly well-suited to this role because of the stability properties of the models that generate them. Kadkhodaie et al. [25], Mlodozeniec et al. [44] showed that these models are stable in the sense that sampling from the same seed and the same latent vector yields visually similar outputs across different trained models, reflecting a stability of the learned score, or velocity, field whenever the target distribution is held fixed. Briq et al. [9] showed an analogous stability of the path dynamics itself under ODE-based sampling. This might suggest that it suffices to test only the initial and final points of the path, but doing so discards information. English and Suzuki [16] showed that paths can differ substantially even when their final mappings are visually near-identical, and, more importantly, such endpoint similarity can itself be an artefact of modelling bias rather than genuine distributional closeness. Additionally, these models are known to lock onto low-dimensional data manifolds in high-dimensional settings [46, 8], making them well-suited to high-dimensional inference problems. Intuitively, path-based inference lets us base statistical conclusions on the entire generative process connecting two datasets, rather than on their endpoints alone. Figure 1: Path-based two-sample testing with a Gaussian-bottleneck interpolant. Samples from ρ0 _0 are transported by the probability flow (solid curves) to the shared bottleneck (0,12d)N(0, 12I_d) at t=12t= 12 (open circles) and onwards to ρ1 _1, while dashed curves show the first half of each path mirrored under t↦1−t 1-t. Left: under ℋ0H_0, the path is reflection symmetric about the bottleneck, so matched slices such as ρ1/4 _1/4 and ρ3/4 _3/4 coincide, each trajectory retraces its first half, and the mirrored curves are hidden beneath the flow. Right: under ℋ1H_1, the second half of each path departs from the mirror image of the first, and this asymmetry, measured through the denoiser and velocity fields at matched noise levels, is the two-sample signal. We develop this idea for the classical problem of two-sample testing [35]: given two finite datasets, determine whether they were drawn from the same distribution. Two-sample testing underlies many statistical and machine learning tasks, including detecting distribution shift [47], checking whether generated data match the true data distribution [51], and testing whether a distribution differs across demographic groups [4]. Although classical two-sample testing is well studied [21, 3], existing methods still struggle with structured, high-dimensional data, especially in the low-sample regimes common in resource-limited settings such as biomedicine [48]. A further desideratum is interpretability: beyond a binary reject or fail-to-reject decision, we would like to identify which observations drive a rejection, for instance, to flag and remove corrupted data from a newly collected batch while preserving the rest of the distribution. Path-based inference is naturally suited to this, since our test statistics are differences of averages of observation-level scores, so the observations that contribute most to a rejection can be read off directly (Figures 7 and 5). For a comprehensive positioning with respect to related works, see Section 6. We build on the stochastic interpolant framework, which unifies the major dynamical-transport generative models under a single formalism, and make the following contributions: (i) We show that an interpolant routing both distributions through a shared Gaussian bottleneck [1] turns two-sample testing into a test of path symmetry, and we prove that the null hypothesis holds if and only if the denoiser, or equivalently velocity, fields of its two halves coincide at any single noise level (Section 3.3); (i) we turn this characterisation into denoising-risk and velocity-risk statistics that are evaluated on held-out data and aggregated along the path, and we show that an information-theoretic weighting recovers the Jeffreys divergence between the noise-smoothed distributions (Section 3.5); (i) we prove that permutation calibration yields tests that are valid in finite samples for any trained networks and consistent when the fields are learned accurately, with an explicit account of how estimation errors affect power (Sections 3.6 and 3.7); and (iv) on a synthetic benchmark and three image benchmarks, the proposed tests improve power over the strongest baseline by up to 33 percentage points at an equal total sample budget (Section 4). 2 Preliminaries 2.1 Stochastic Interpolants We work within the stochastic interpolant framework of Albergo et al. [1], which unifies flow-based and diffusion-based generative models under a single construction. Let ρ0 _0 and ρ1 _1 be probability densities on ℝdR^d, and let (x0,x1)∼ν(x_0,x_1) ν be a coupling with marginals ρ0 _0 and ρ1 _1; the simplest choice is the independent coupling ν=ρ0⊗ρ1ν= _0 _1. Let z∼(0,d)z (0,I_d) be a latent variable independent of (x0,x1)(x_0,x_1). The stochastic interpolant is the process xt=αtx0+βtx1+γtz,t∈[0,1],x_t= _t\,x_0+ _t\,x_1+ _t\,z, t∈[0,1], (1) where the coefficients αt,βt,γt _t, _t, _t are sufficiently smooth and satisfy the boundary conditions α0=β1=1,α1=β0=0,γ0=γ1=0,γt>0for t∈(0,1). _0= _1=1, _1= _0=0, _0= _1=0, _t>0\;\;for t∈(0,1). (2) These conditions guarantee that xt=0=x0x_t=0=x_0 and xt=1=x1x_t=1=x_1, so the time-dependent density ρt _t of xtx_t interpolates between ρ0 _0 and ρ1 _1. A standard choice is αt=1−t _t=1-t, βt=t _t=t and γt=2t(1−t) _t= 2t(1-t), and the same latent coefficient γt _t is used in our construction in Section 3. The latent term γtz _tz smooths ρt _t at intermediate times; setting γt≡0 _t≡ 0 recovers the deterministic interpolants underlying flow matching [37] and rectified flow [40]. Equation 1 specifies the path at the level of individual samples: conditional on (x0,x1,z)(x_0,x_1,z), the trajectory t↦xtt x_t is fully determined, and its law defines the modelling path (ρt)t∈[0,1]( _t)_t∈[0,1]. The conditional velocity along this path is available in closed form, x˙t=α˙tx0+β˙tx1+γ˙tz, x_t= α_t\,x_0+ β_t\,x_1+ γ_t\,z, (3) where dots denote time derivatives. However, x˙t x_t depends on the endpoints and the latent, which are not accessible from xtx_t alone. The generative dynamics are therefore obtained by averaging over them. Velocity field: The marginal velocity field is the conditional expectation of equation 3 given the current position, b(t,x)=[x˙t|xt=x].b(t,x)=E\! [ x_t\, |\,x_t=x ]. (4) Albergo and Vanden-Eijnden [2] show that ρt _t satisfies the transport equation ∂tρt+∇⋅(bρt)=0,ρt=0=ρ0. _t _t+∇· (b\, _t )=0, _t=0= _0. (5) Consequently, solutions of the probability flow ODE X˙t=b(t,Xt),X0∼ρ0, X_t=b(t,X_t), X_0 _0, (6) satisfy Xt∼ρtX_t _t for all t∈[0,1]t∈[0,1], and in particular X1∼ρ1X_1 _1. Integrating equation 6 backward in time from X1∼ρ1X_1 _1 likewise transports ρ1 _1 to ρ0 _0. Although the interpolant equation 1 involves the latent z, the transport equation 6 is deterministic: stochasticity enters only through the construction of ρt _t, and each initial condition X0X_0 is mapped to a unique trajectory. Denoisers and score: Equivalently, the dynamics can be parameterised through the denoisers η0(t,x)=[x0|xt=x],η1(t,x)=[x1|xt=x],ηz(t,x)=[z|xt=x]. _0(t,x)=E\! [x_0\, |\,x_t=x ], _1(t,x)=E\! [x_1\, |\,x_t=x ], _z(t,x)=E\! [z\, |\,x_t=x ]. (7) By linearity of conditional expectation, b(t,x)=α˙tη0(t,x)+β˙tη1(t,x)+γ˙tηz(t,x),x=αtη0(t,x)+βtη1(t,x)+γtηz(t,x),b(t,x)= α_t\, _0(t,x)+ β_t\, _1(t,x)+ γ_t\, _z(t,x), x= _t\, _0(t,x)+ _t\, _1(t,x)+ _t\, _z(t,x), (8) so that only two of the three denoisers need to be learned. Moreover, for Gaussian z, the latent denoiser determines the score of ρt _t, ∇logρt(x)=−γt−1ηz(t,x),t∈(0,1),∇ _t(x)=- _t^-1\, _z(t,x), t∈(0,1), (9) which connects the framework to score-based diffusion models. For any ϵ≥0ε≥ 0, the SDE dXt=[b(t,Xt)+ϵ∇logρt(Xt)]dt+2ϵdWtdX_t= [b(t,X_t)+ε\,∇ _t(X_t) ]\,dt+ 2ε\,dW_t shares the marginals ρt _t; throughout this work, we use only the deterministic case ϵ=0ε=0. Learning: Each conditional expectation above is the unique minimiser of a quadratic regression objective. The velocity field minimises ℒb[b^]=∫01|b^(t,xt)−x˙t|2t,L_b[ b]= _0^1E\, | b(t,x_t)- x_t |^2\,dt, (10) and the latent denoiser minimises ℒηz[η^z]=∫01|η^z(t,xt)−z|2t,L_ _z[ η_z]= _0^1E\, | η_z(t,x_t)-z |^2\,dt, (11) with analogous objectives for η0 _0 and η1 _1 obtained by replacing z with x0x_0 or x1x_1. Here, the expectation is taken over (x0,x1)∼ν(x_0,x_1) ν and z∼(0,d)z (0,I_d), with xtx_t and x˙t x_t given by equation 1 and equation 3. These objectives coincide, up to additive constants independent of b b and η^z η_z, with those of Albergo and Vanden-Eijnden [2]. In practice, b b and η^z η_z are parameterised by neural networks, and the objectives are estimated by Monte Carlo sampling of t∼[0,1]t [0,1], the endpoints and the latent. 2.2 Two-Sample Testing Consider two independent datasets 0=x0(i)i=1n0D_0=\x_0^(i)\_i=1^n_0 and 1=x1(j)j=1n1D_1=\x_1^(j)\_j=1^n_1 in ℝdR^d, whose members are drawn i.i.d. from densities ρ0 _0 and ρ1 _1, respectively. The densities themselves are unknown, and we only have access to the two finite datasets. The two-sample problem asks whether these datasets are consistent with a single underlying density, which we phrase as a test of the null hypothesis against its complement, ℋ0:ρ0=ρ1versusℋ1:ρ0≠ρ1.H_0 _0= _1 _1 _0≠ _1. (12) Our validity guarantee is non-asymptotic, holding at finite n0n_0 and n1n_1, whereas power is analysed as the sample sizes grow. A test is a binary decision rule Δ(0,1)∈0,1 (D_0,D_1)∈\0,1\, where Δ=1 =1 corresponds to rejecting ℋ0H_0. We allow Δ to be randomised, that is, to depend on internal randomness ξ drawn independently of the data, such as the sample splitting, network training, Monte Carlo draws, and permutations used in Section 3. We write ℙρ0,ρ1P_ _0, _1 for the joint law of (0,1,ξ)(D_0,D_1,ξ) when the datasets are drawn from ρ0 _0 and ρ1 _1; the law of ξ is fixed by the test and is therefore left implicit. More generally, subscripts on ℙP and E indicate the distributions of the random variables involved, as in x∼ρ[f(x)]E_x ρ[f(x)], and unsubscripted probabilities and expectations are taken over all random variables in the expression. A false rejection, or type I error, occurs when Δ=1 =1 although ρ0=ρ1 _0= _1. We require its probability to be controlled uniformly over the set (ℝd)P(R^d) of probability densities on ℝdR^d, supρ0∈(ℝd)ℙρ0,ρ0[Δ(0,1)=1]≤κ, _ _0 (R^d)P_ _0, _0 [ (D_0,D_1)=1 ]≤κ, (13) for a prescribed κ∈(0,1)κ∈(0,1), in which case Δ is said to have level κ. Conversely, a type I error occurs when Δ=0 =0 although ρ0≠ρ1 _0≠ _1. For a fixed alternative, if ℙρ0,ρ1[Δ(0,1)=0]≤δP_ _0, _1 [ (D_0,D_1)=0 ]≤δ (14) for some δ∈(0,1)δ∈(0,1), the test is said to have power 1−δ1-δ against that alternative. The goal is to construct tests that attain level κ while maximising power against a broad class of alternatives. 3 Path-Based Two-Sample Testing We work in the setting of Section 2.2. Our construction connects each of ρ0 _0 and ρ1 _1 to a common Gaussian bottleneck that depends on neither population. Rather than comparing ρ0 _0 and ρ1 _1 only at their endpoints, we compare the regression fields induced along the two halves of this path. Under ℋ0H_0, the two fields coincide at every matched noise level. In contrast, under ℋ1H_1, they differ at every matched noise level, so their discrepancy supplies a continuum of two-sample witnesses indexed by position along the path. Throughout, population fields are distinguished from their neural estimates, which we denote with hats. 3.1 A Gaussian-Bottleneck Interpolant Following Albergo et al. [1], we specialise the interpolant equation 1 to the schedule αt=cos2(πt) 1t≤1/2,βt=cos2(πt) 1t≥1/2,γt=2t(1−t). _t= ^2(π t)\,1_\t≤ 1/2\, _t= ^2(π t)\,1_\t≥ 1/2\, _t= 2t(1-t). (15) Its structure is most transparent in the radial coordinate s=|t−12|∈[0,12]s=|t- 12|∈[0, 12], which measures distance from the midpoint. Each radius s corresponds to the two mirrored times ts0=12−s,ts1=12+s,t_s^0= 12-s, t_s^1= 12+s, (16) and we define the radial schedule as=sin2(πs),σs=12−2s2.a_s= ^2(π s), _s= 12-2s^2. (17) For s∈(0,12)s∈(0, 12), this schedule satisfies as>0a_s>0 and σs>0 _s>0, and its derivatives as′a_s and σs′ _s are finite; we use these facts without further mention. The following lemma records the properties of the schedule that we need. Lemma 1 (Gaussian-bottleneck interpolant). The schedule equation 15 satisfies the boundary conditions equation 2, with αt,βt∈C1([0,1]) _t, _t∈ C^1([0,1]) and γt2∈C∞([0,1]) _t^2∈ C^∞([0,1]). Moreover, for every s∈[0,12]s∈[0, 12] and ℓ∈0,1 ∈\0,1\, xtsℓ=asxℓ+σsz.x_t_s =a_s\,x_ + _s\,z. (18) In particular, the law of xtx_t does not depend on the coupling ν, and the midpoint density is the pivot ρ1/2=(0,12d) _1/2=N (0, 12I_d ). Proof. The proof consists of a direct evaluation of the schedule. For the boundary values, α0=cos20=1 _0= ^20=1, α1=0 _1=0, β0=0 _0=0, β1=cos2π=1 _1= ^2π=1, and γ0=γ1=0 _0= _1=0, while γt2=2t(1−t)>0 _t^2=2t(1-t)>0 for t∈(0,1)t∈(0,1). The function γt2 _t^2 is a polynomial and hence lies in C∞([0,1])C^∞([0,1]). On [0,12][0, 12], αt=cos2(πt) _t= ^2(π t) has derivative −πsin(2πt)-π (2π t), and both αt _t and its derivative vanish at t=12t= 12, matching the identically zero branch on (12,1]( 12,1]. Hence αt∈C1([0,1]) _t∈ C^1([0,1]), and the same argument applies to βt _t. For the channel representation, let s∈[0,12]s∈[0, 12]. At t=ts0=12−st=t_s^0= 12-s, we have αt=cos2(π2−πs)=sin2(πs)=as _t= ^2( π2-π s)= ^2(π s)=a_s and βt=0 _t=0, since either t<12t< 12 or t=12t= 12, where cos2(π2)=0 ^2( π2)=0. Moreover, γt2=2(12−s)(12+s)=12−2s2=σs2 _t^2=2( 12-s)( 12+s)= 12-2s^2= _s^2. Substituting into equation 1 gives xts0=asx0+σszx_t_s^0=a_sx_0+ _sz. At t=ts1t=t_s^1, the identity cos2(π2+πs)=sin2(πs) ^2( π2+π s)= ^2(π s) gives xts1=asx1+σszx_t_s^1=a_sx_1+ _sz in the same way. Since z is independent of (x0,x1)(x_0,x_1), the law of xtsℓx_t_s depends only on the marginal ρℓ _ and not on the coupling ν. At s=0s=0, we have a0=0a_0=0 and σ02=12 _0^2= 12, so x1/2=σ0z∼(0,12d)x_1/2= _0z (0, 12I_d). ∎ Lemma 1 shows that each half of the path is a Gaussian channel acting on a single population. We write ζsℓ=ρtsℓ,ℓ∈0,1, _s = _t_s , ∈\0,1\, (19) for the density of channel ℓ at radius s, so that ζs0 _s^0 and ζs1 _s^1 are the interpolant densities at the two mirrored times. As s↑12s 12, we have as→1a_s→ 1 and σs→0 _s→ 0, so the channels recover the data densities, ζ1/2ℓ=ρℓ _1/2 = _ , while at s=0s=0 both reduce to the pivot, ζ00=ζ01=(0,12d) _0^0= _0^1=N (0, 12I_d ). In the original time coordinate, the construction therefore realises the transport ρ0⟶(0,12d)⟵ρ1, _0 (0, 12I_d ) _1, (20) while in the radial coordinate, it consists of two outward paths emanating from the same pivot. The coefficients αt _t and βt _t are only once differentiable at the bottleneck. This suffices for all results below, which we establish directly for the schedule equation 15 rather than through the general theory of Albergo et al. [1]. Throughout, we impose the following condition on the two populations. Assumption 1 (Moments). The densities ρ0 _0 and ρ1 _1 have finite second moments. Assumption 1 guarantees that the fields of Section 3.2 are well defined and that the squared-error risks of Section 3.4 are finite. The Jeffreys identity of Proposition 7 additionally relies on a technical regularity condition, Assumption 2, which we state in Section 3.5. 3.2 Denoiser and Velocity Fields For ℓ∈0,1 ∈\0,1\ and s∈(0,12)s∈(0, 12), we define the posterior-mean denoiser and the outward velocity of channel ℓ as ηsℓ(x)=[xℓ|xtsℓ=x],usℓ(x)=[∂sxtsℓ|xtsℓ=x], _s (x)=E\! [x_ \, |\,x_t_s =x ], u_s (x)=E\! [ _sx_t_s \, |\,x_t_s =x ], (21) where ∂sxtsℓ=as′xℓ+σs′z _sx_t_s =a_s x_ + _s z and primes denote derivatives in s. These are the fields of Section 2.1 expressed in the radial coordinate: ηs0=η0(ts0,⋅) _s^0= _0(t_s^0,·) and ηs1=η1(ts1,⋅) _s^1= _1(t_s^1,·), while us0=−b(ts0,⋅)u_s^0=-b(t_s^0,·) and us1=b(ts1,⋅)u_s^1=b(t_s^1,·), the sign reflecting that s runs against t on the first half. Under Assumption 1, both conditional expectations exist and admit continuous versions on ℝdR^d, constructed for the denoiser in Lemma 2 and for the velocity in the proof of Lemma 3. We use these versions throughout, so that identities between fields hold at every point. Write φσ(x)=(2πσ2)−d/2exp(−|x|2/(2σ2)) _σ(x)=(2πσ^2)^-d/2 \! (-|x|^2/(2σ^2) ) for the density of (0,σ2d)N(0,σ^2I_d). The following lemma collects the properties of the channel densities used throughout. Lemma 2 (Gaussian smoothing). Let s∈[0,12)s∈[0, 12) and ℓ∈0,1 ∈\0,1\. The channel density is ζsℓ(x)=∫ℝdφσs(x−asy)ρℓ(y)y, _s (x)= _R^d _ _s(x-a_sy)\, _ (y)\,dy, (22) which is infinitely differentiable and strictly positive on ℝdR^d. If, in addition, s>0s>0 and Assumption 1 holds, then the posterior density of xℓx_ given xtsℓ=x_t_s =x is φσs(x−asy)ρℓ(y)/ζsℓ(x) _ _s(x-a_sy)\, _ (y)/ _s (x), and ηsℓ(x)=1ζsℓ(x)∫ℝdyφσs(x−asy)ρℓ(y)y _s (x)= 1 _s (x) _R^dy\, _ _s(x-a_sy)\, _ (y)\,dy (23) is a continuous version of the conditional expectation in equation 21. Proof. By Lemma 1, xtsℓ=asxℓ+σszx_t_s =a_sx_ + _sz with z∼(0,d)z (0,I_d) independent of xℓx_ . Conditionally on xℓ=yx_ =y, the variable xtsℓx_t_s therefore has density φσs(⋅−asy) _ _s(·-a_sy), and integrating over y∼ρℓy _ gives equation 22. Every partial derivative of φσs _ _s, of any order, is bounded on ℝdR^d. Hence, the corresponding derivative of the integrand in equation 22 is dominated by a constant multiple of ρℓ(y) _ (y), uniformly in x, and dominated convergence allows differentiation under the integral sign to all orders. This shows ζsℓ∈C∞(ℝd) _s ∈ C^∞(R^d). Strict positivity follows because φσs>0 _ _s>0 everywhere and ρℓ _ integrates to one. The posterior density follows from Bayes’ rule. Finally, the integrand in the numerator of equation 23 is continuous in x and bounded by (2πσs2)−d/2|y|ρℓ(y)(2π _s^2)^-d/2\,|y|\, _ (y), which is integrable under Assumption 1. The numerator is therefore continuous by dominated convergence, and so is the ratio, since the denominator is continuous and strictly positive. ∎ Lemma 3 (Denoiser, velocity and score). Under Assumption 1, for every s∈(0,12)s∈(0, 12) and ℓ∈0,1 ∈\0,1\: (i) the velocity satisfies usℓ(x)=λsηsℓ(x)+σs′σsxu_s (x)= _s\, _s (x)+ _s _s\,x with λs=as′−asσs′σs=πsin(2πs)+2sasσs2 _s=a_s -a_s _s _s=π (2π s)+ 2s\,a_s _s^2; (i) the denoiser satisfies Tweedie’s identity, ηsℓ(x)=1as(x+σs2∇logζsℓ(x)) _s (x)= 1a_s (x+ _s^2\,∇ _s (x) ); and (i) the contrasts satisfy us1−us0=λs(ηs1−ηs0)=χs(∇logζs1−∇logζs0)u_s^1-u_s^0= _s( _s^1- _s^0)= _s(∇ _s^1-∇ _s^0) with χs=λsσs2as>0 _s= _s _s^2a_s>0. Proof. Parts (i) and (i) are computed directly, and part (i) follows by subtraction. Fix s∈(0,12)s∈(0, 12) and ℓ∈0,1 ∈\0,1\, and abbreviate xt=xtsℓx_t=x_t_s . Both ∂sxt=as′xℓ+σs′z _sx_t=a_s x_ + _s z and xℓx_ are integrable under Assumption 1, so all conditional expectations below exist. For (i), linearity of conditional expectation gives usℓ(x)=as′ηsℓ(x)+σs′[z∣xt=x]u_s (x)=a_s \, _s (x)+ _s \,E[z x_t=x]. Since z=(xt−asxℓ)/σsz=(x_t-a_sx_ )/ _s, we have [z∣xt=x]=(x−asηsℓ(x))/σsE[z x_t=x]=(x-a_s _s (x))/ _s. Substituting and collecting terms yields (i). The right-hand side of (i) is continuous by Lemma 2, and it provides the continuous version of usℓu_s used throughout. The explicit form of λs _s follows from as′=2πsin(πs)cos(πs)=πsin(2πs)a_s =2π (π s) (π s)=π (2π s) and from differentiating σs2=12−2s2 _s^2= 12-2s^2, which gives σs′/σs=−2s/σs2 _s / _s=-2s/ _s^2. For (i), differentiating equation 22 under the integral sign, as justified in Lemma 2, gives ∇ζsℓ(x)=−1σs2∫ℝd(x−asy)φσs(x−asy)ρℓ(y)dy=−xσs2ζsℓ(x)+asσs2ζsℓ(x)ηsℓ(x),∇ _s (x)=- 1 _s^2 _R^d(x-a_sy)\, _ _s(x-a_sy)\, _ (y)\,dy=- x _s^2\, _s (x)+ a_s _s^2\, _s (x)\, _s (x), (24) where the second equality uses equation 23. Dividing by ζsℓ(x)>0 _s (x)>0 and rearranging yields (i). For (i), subtracting (i) for ℓ=0 =0 from (i) for ℓ=1 =1 cancels the common affine term and gives the first equality in (i). Subtracting (i) for the two channels gives ηs1−ηs0=(σs2/as)(∇logζs1−∇logζs0) _s^1- _s^0=( _s^2/a_s)(∇ _s^1-∇ _s^0), and multiplying by λs _s gives the second equality. Finally, on (0,12)(0, 12) we have sin(2πs)>0 (2π s)>0, as>0a_s>0 and σs2>0 _s^2>0, so λs>0 _s>0 and hence χs>0 _s>0. ∎ The affine term in Lemma 3 (i) is known and common to both channels, so it cancels in every contrast. Denoiser, velocity and score contrasts therefore encode the same population difference, up to known positive scalar functions of s. 3.3 Path Invariance Characterises the Null Proposition 4 (Path invariance). Suppose Assumption 1 holds, and fix any s∈(0,12)s∈(0, 12). Then the following are equivalent: (i) ℋ0H_0 holds (ρ0=ρ1 _0= _1); (i) ζs0=ζs1 _s^0= _s^1; (i) ηs0=ηs1 _s^0= _s^1; (iv) us0=us1u_s^0=u_s^1. Proof. We prove the cycle (i) ⇒ (i) ⇒ (i) ⇒ (i), and separately the equivalence (i) ⇔ (iv). The key step is (i) ⇒ (i), which inverts the Gaussian channel. The implication (i) ⇒ (i) holds because, by Lemma 1, the law of xtsℓx_t_s is determined by ρℓ _ alone, so ρ0=ρ1 _0= _1 implies ζs0=ζs1 _s^0= _s^1. For (i) ⇒ (i), the right-hand side of Lemma 3 (i) depends on channel ℓ only through ζsℓ _s , so equal densities give equal denoisers. For (i) ⇒ (i), Lemma 3 (i) shows that equal denoisers give ∇logζs0=∇logζs1∇ _s^0=∇ _s^1 on ℝdR^d. Both densities are smooth and strictly positive by Lemma 2, so log(ζs0/ζs1) ( _s^0/ _s^1) is constant on the connected set ℝdR^d. Since both densities integrate to one, the constant is zero, and hence ζs0=ζs1 _s^0= _s^1. Writing f^ℓ(k)=x∼ρℓ[eik⋅x] f_ (k)=E_x _ [e^ik· x ] for the characteristic function of ρℓ _ , the representation equation 18 gives the characteristic function of ζsℓ _s as f^ℓ(ask)e−σs2|k|2/2 f_ (a_sk)\,e^- _s^2|k|^2/2. The Gaussian factor never vanishes, so f^0(ask)=f^1(ask) f_0(a_sk)= f_1(a_sk) for all k∈ℝdk ^d. Since as>0a_s>0, this gives f^0=f^1 f_0= f_1, and Lévy’s uniqueness theorem yields ρ0=ρ1 _0= _1. Finally, (i) ⇔ (iv) follows from Lemma 3 (i): since us1−us0=λs(ηs1−ηs0)u_s^1-u_s^0= _s( _s^1- _s^0) with λs>0 _s>0, one difference vanishes if and only if the other does. ∎ In the original time coordinate, condition (iv) is a reflection antisymmetry of the velocity field about the bottleneck, b(12−s,x)+b(12+s,x)=0,b( 12-s,x)+b( 12+s,x)=0, (25) and condition (i), which reads ρ1/2−s=ρ1/2+s _1/2-s= _1/2+s, is the corresponding reflection symmetry of the interpolant density. Consequently, whenever the probability flow ODE equation 6 is well posed, every trajectory satisfies X1−t=XtX_1-t=X_t under ℋ0H_0, since t↦X1−t X_1-t solves the same ODE by equation 25 and agrees with XtX_t at t=12t= 12. The transport from ρ0 _0 to ρ1 _1 is then the identity map, which makes the optimal-transport intuition of Section 1 exact for this construction. Proposition 4 has two consequences. Under ℋ0H_0, the path fields coincide at every radius. Under ℋ1H_1, they differ at every radius, so any single matched noise level already carries two-sample evidence. This motivates testing ℋ0H_0 through discrepancies between the fields of the two channels. 3.4 Held-Out Risk Witnesses Suppose that η^s0,u^s0 η_s^0, u_s^0 and η^s1,u^s1 η_s^1, u_s^1 are estimated on a training split that is independent of the observations used for testing. For a held-out observation ω∈ℝdω ^d, we draw an independent latent z∼(0,d)z (0,I_d) and form the corrupted point and its exact outward derivative, ωs=asω+σsz,Us(ω,z)=as′ω+σs′z. _s=a_s\,ω+ _s\,z, U_s(ω,z)=a_s \,ω+ _s \,z. (26) We define the denoising-risk and velocity-risk contrasts ψsden(ω,z) _s^den(ω,z) =|ω−η^s1(ωs)|2−|ω−η^s0(ωs)|2, = |ω- η_s^1( _s) |^2- |ω- η_s^0( _s) |^2, (27) ψsvel(ω,z) _s^vel(ω,z) =|Us(ω,z)−u^s1(ωs)|2−|Us(ω,z)−u^s0(ωs)|2. = |U_s(ω,z)- u_s^1( _s) |^2- |U_s(ω,z)- u_s^0( _s) |^2. (28) Each contrast is a single scalar function evaluated identically on held-out observations from both datasets. The group label of ω is never supplied to the witness; the labels 00 and 11 in equation 27 and equation 28 refer only to the two learned fields. Throughout, the latent z∼(0,d)z (0,I_d) is drawn independently of all other variables. We write, for instance, ω∼ρ0,zE_ω _0,\,z for the expectation over a held-out observation ω∼ρ0ω _0 and the latent. Lemma 5 (Population risk contrasts). Suppose Assumption 1 holds and that the learned fields equal their population counterparts. Then for every s∈(0,12)s∈(0, 12), ω∼ρ0,z[ψsden(ω,z)] _ω _0,\,z [ _s^den(ω,z) ] =x∼ζs0[|ηs0(x)−ηs1(x)|2], =E_x _s^0 [| _s^0(x)- _s^1(x)|^2 ], ω∼ρ1,z[ψsden(ω,z)] _ω _1,\,z [ _s^den(ω,z) ] =−x∼ζs1[|ηs0(x)−ηs1(x)|2], =-E_x _s^1 [| _s^0(x)- _s^1(x)|^2 ], (29) ω∼ρ0,z[ψsvel(ω,z)] _ω _0,\,z [ _s^vel(ω,z) ] =x∼ζs0[|us0(x)−us1(x)|2], =E_x _s^0 [|u_s^0(x)-u_s^1(x)|^2 ], ω∼ρ1,z[ψsvel(ω,z)] _ω _1,\,z [ _s^vel(ω,z) ] =−x∼ζs1[|us0(x)−us1(x)|2], =-E_x _s^1 [|u_s^0(x)-u_s^1(x)|^2 ], (30) where both sides may be infinite. Consequently, the population differences of held-out means, Dsrep=ω∼ρ0,z[ψsrep(ω,z)]−ω∼ρ1,z[ψsrep(ω,z)],rep∈den,vel,D_s^rep=E_ω _0,\,z [ _s^rep(ω,z) ]-E_ω _1,\,z [ _s^rep(ω,z) ], ∈\den,vel\, (31) take values in [0,∞][0,∞] and satisfy Dsvel=λs2DsdenD_s^vel= _s^2D_s^den. Moreover, Dsrep=0D_s^rep=0 if and only if ℋ0H_0 holds. Proof. The identities follow from the orthogonality of regression residuals, and the final claim follows from Proposition 4. Fix s∈(0,12)s∈(0, 12), and consider first the denoising contrast for ω∼ρ0ω _0; in this paragraph, E denotes the expectation over ω∼ρ0ω _0 and z. By Lemma 1, the pair (ω,ωs)(ω, _s) has the same joint law as (x0,xts0)(x_0,x_t_s^0), so ωs∼ζs0 _s _s^0 and ηs0(ωs)=[ω∣ωs] _s^0( _s)=E[ω _s]. Write e=ω−ηs0(ωs)e=ω- _s^0( _s) and g=ηs0(ωs)−ηs1(ωs)g= _s^0( _s)- _s^1( _s), so that ω−ηs1(ωs)=e+gω- _s^1( _s)=e+g and ψsden(ω,z)=|e+g|2−|e|2 _s^den(ω,z)=|e+g|^2-|e|^2. Since e is the residual of a conditional expectation, |e|2=|ω|2−|ηs0(ωs)|2≤ω∼ρ0|ω|2<∞E|e|^2=E|ω|^2-E| _s^0( _s)|^2 _ω _0|ω|^2<∞ under Assumption 1. The negative part of ψsden _s^den is therefore integrable, and its expectation is well defined in (−∞,∞](-∞,∞]. If |g|2<∞E|g|^2<∞, then the cross term [e⋅g]E[e· g] is integrable by the Cauchy–Schwarz inequality and vanishes by the tower property, because g is a function of ωs _s and [e∣ωs]=0E[e _s]=0. Expanding the square then gives [ψsden(ω,z)]=|g|2=x∼ζs0[|ηs0(x)−ηs1(x)|2],E [ _s^den(ω,z) ]=E|g|^2=E_x _s^0 [| _s^0(x)- _s^1(x)|^2 ], (32) which is the first identity in equation 29. If instead |g|2=∞E|g|^2=∞, then e+ge+g is not square integrable, because e is, and both sides equal +∞+∞. For ω∼ρ1ω _1, the pair (ω,ωs)(ω, _s) has the same law as (x1,xts1)(x_1,x_t_s^1), so ωs∼ζs1 _s _s^1. The same argument with the roles of the two denoisers exchanged gives ω∼ρ1,z[ψsden(ω,z)]=−x∼ζs1[|ηs0(x)−ηs1(x)|2]∈[−∞,0]E_ω _1,\,z[ _s^den(ω,z)]=-E_x _s^1[| _s^0(x)- _s^1(x)|^2]∈[-∞,0], which is the second identity in equation 29. The velocity identities equation 30 follow identically, replacing ω by the target Us(ω,z)U_s(ω,z). This target is square integrable under Assumption 1, and usℓ(ωs)=[Us(ω,z)∣ωs]u_s ( _s)=E[U_s(ω,z) _s] when ω∼ρℓω _ . Since the first identity in each row is non-negative and the second is non-positive, their difference DsrepD_s^rep is well defined in [0,∞][0,∞]. By Lemma 3 (i), |us0−us1|2=λs2|ηs0−ηs1|2|u_s^0-u_s^1|^2= _s^2| _s^0- _s^1|^2 pointwise, so Dsvel=λs2DsdenD_s^vel= _s^2D_s^den. If ℋ0H_0 holds, then ηs0=ηs1 _s^0= _s^1 by Proposition 4, and Dsden=0D_s^den=0. Conversely, if Dsden=0D_s^den=0, then ηs0=ηs1 _s^0= _s^1 holds ζs0 _s^0-almost everywhere. This equality then holds Lebesgue-almost everywhere because ζs0>0 _s^0>0, and everywhere because both fields are continuous. Proposition 4 then gives ℋ0H_0. The same conclusion holds for DsvelD_s^vel, since λs>0 _s>0. ∎ 3.5 Path Aggregation We aggregate the local contrasts over the truncated interval ℐτ=[τ,12−τ]I_τ=[τ, 12-τ], for some τ∈(0,14)τ∈(0, 14), using a positive weight w(s)w(s), Φwrep(ω)=z[∫ℐτw(s)ψsrep(ω,z)s],rep∈den,vel. _w^rep(ω)=E_z\! [ _I_τw(s)\, _s^rep(ω,z)\,ds ], ∈\den,vel\. (33) We consider two weighting schemes. The uniform scheme sets wU(s)=1w_U(s)=1 and assigns equal path measure to all retained radii. The Jeffreys scheme uses the information-theoretically canonical weights wJden(s)=λs2χs=λsasσs2,wJvel(s)=1χs,w_J^den(s)= _s^2 _s= _s\,a_s _s^2, w_J^vel(s)= 1 _s, (34) whose origin is made precise in Proposition 7. Given held-out subsets of 0D_0 and 1D_1 of sizes n0ten_0^te and n1ten_1^te, the empirical statistic is the difference of held-out group means, T^wrep=1n0te∑i=1n0teΦ^wrep(x0(i))−1n1te∑j=1n1teΦ^wrep(x1(j)), T_w^rep= 1n_0^te _i=1^n_0^te _w^rep (x_0^(i) )- 1n_1^te _j=1^n_1^te _w^rep (x_1^(j) ), (35) where Φ^wrep _w^rep denotes the Monte Carlo estimate of Section 3.9. Normalising w by a constant is immaterial, since a deterministic rescaling does not affect the permutation test of Section 3.6. Corollary 6 (Population target). Suppose the conditions of Lemma 5 hold, that the sampling density μ of Section 3.9 is strictly positive on ℐτI_τ, and that ∫ℐτw(s)ω∼ρℓ,z[|ψsrep(ω,z)|]s<∞ _I_τw(s)\,E_ω _ ,\,z [| _s^rep(ω,z)| ]\,ds<∞ for ℓ∈0,1 ∈\0,1\. Then the population mean of T^wrep T_w^rep is Twrep=∫ℐτw(s)Dsreps≥0,T_w^rep= _I_τw(s)\,D_s^rep\,ds≥ 0, (36) with equality if and only if ℋ0H_0 holds. Proof. We compute the mean of a single Monte Carlo term, then apply linearity. Let ω∼ρℓω _ , s∼μs μ and z∼(0,d)z (0,I_d) be independent. By Tonelli’s theorem, the integrability condition in the statement gives ω∼ρℓ,s∼μ,z[w(s)μ(s)|ψsrep(ω,z)|]=∫ℐτw(s)ω∼ρℓ,z[|ψsrep(ω,z)|]s<∞,E_ω _ ,\,s μ,\,z\! [ w(s)μ(s)\, | _s^rep(ω,z) | ]= _I_τw(s)\,E_ω _ ,\,z [| _s^rep(ω,z)| ]\,ds<∞, (37) where we used μ>0μ>0 on ℐτI_τ. Fubini’s theorem then gives ω∼ρℓ,s∼μ,z[w(s)μ(s)ψsrep(ω,z)]=∫ℐτw(s)ω∼ρℓ,z[ψsrep(ω,z)]s,E_ω _ ,\,s μ,\,z\! [ w(s)μ(s)\, _s^rep(ω,z) ]= _I_τw(s)\,E_ω _ ,\,z [ _s^rep(ω,z) ]\,ds, (38) and the inner expectations are finite for almost every s∈ℐτs _τ. Each term of equation 70 has this mean for an observation from ρℓ _ , so by linearity the mean of T^wrep T_w^rep is ∫ℐτw(s)(ω∼ρ0,z[ψsrep(ω,z)]−ω∼ρ1,z[ψsrep(ω,z)])s=∫ℐτw(s)Dsreps. _I_τw(s) (E_ω _0,\,z [ _s^rep(ω,z) ]-E_ω _1,\,z [ _s^rep(ω,z) ] )\,ds= _I_τw(s)\,D_s^rep\,ds. (39) The integrand is non-negative by Lemma 5, and w>0w>0, so the integral vanishes if and only if Dsrep=0D_s^rep=0 for almost every s∈ℐτs _τ. If ℋ0H_0 holds, then Dsrep=0D_s^rep=0 for every s. Conversely, since τ<14τ< 14, the interval ℐτI_τ has positive length, so if Dsrep=0D_s^rep=0 for almost every s∈ℐτs _τ, it vanishes for at least one s, and Lemma 5 gives ℋ0H_0. ∎ We now make precise the origin of the Jeffreys weights equation 34. The identity below relies on the following technical condition, which is used nowhere else in the paper. Assumption 2 (Regularity of the divergences). For every S∈(0,12)S∈(0, 12), the maps s↦KL(ζs0∥ζs1)s ( _s^0\,\|\, _s^1) and s↦KL(ζs1∥ζs0)s ( _s^1\,\|\, _s^0) are absolutely continuous on every compact subinterval of (0,S](0,S]. Their derivatives may be computed by differentiating under the integral sign and integrating by parts with vanishing boundary terms. For a density ζ on ℝdR^d and a measurable map f:ℝd→ℝdf ^d ^d, we write ‖f‖L2(ζ)=(x∼ζ[|f(x)|2])1/2.\|f\|_L^2(ζ)= (E_x ζ [|f(x)|^2 ] )^1/2. (40) Proposition 7 (Jeffreys identity). Under Assumptions 1 and 2, for every S∈(0,12)S∈(0, 12), S:=∫0SDsvelχsds=∫0SwJden(s)Dsdends=KL(ζS0∥ζS1)+KL(ζS1∥ζS0). J_S:= _0^S D_s^vel _s\,ds= _0^Sw_J^den(s)\,D_s^den\,ds=KL ( _S^0\, \|\, _S^1 )+KL ( _S^1\, \|\, _S^0 ). (41) Proof. The middle equality in equation 41 follows from Lemma 5, since wJden(s)Dsden=(λs2/χs)Dsden=Dsvel/χsw_J^den(s)\,D_s^den=( _s^2/ _s)\,D_s^den=D_s^vel/ _s. For the remaining equality, we derive continuity equations for the two channels, differentiate the resulting divergences, and integrate from the pivot. Within this proof, we write Fs=KL(ζs0∥ζs1)F_s=KL( _s^0\| _s^1) and Gs=KL(ζs1∥ζs0)G_s=KL( _s^1\| _s^0). We begin with the continuity equations. Let f∈Cc∞(ℝd)f∈ C_c^∞(R^d), and let s range over a compact subinterval of (0,12)(0, 12). The map s↦xtsℓs x_t_s is differentiable with derivative ∂sxtsℓ=as′xℓ+σs′z _sx_t_s =a_s x_ + _s z. This derivative is dominated by an integrable variable, uniformly in s, because as′a_s and σs′ _s are bounded on the subinterval and xℓx_ and z are integrable. Since ∇f∇ f is bounded, we may differentiate under the expectation, and conditioning on xtsℓx_t_s gives dsx∼ζsℓ[f(x)]=[∇f(xtsℓ)⋅∂sxtsℓ]=x∼ζsℓ[∇f(x)⋅usℓ(x)]. dds\,E_x _s [f(x) ]=E [∇ f(x_t_s )· _sx_t_s ]=E_x _s [∇ f(x)· u_s (x) ]. (42) Hence ∂sζs0+∇⋅(us0ζs0)=0 _s _s^0+∇·(u_s^0\, _s^0)=0 and ∂sζs1+∇⋅(us1ζs1)=0 _s _s^1+∇·(u_s^1\, _s^1)=0 in the weak sense. Next, we differentiate the divergences. Fix s∈(0,12)s∈(0, 12). Assumption 2 permits differentiating Fs=∫ζs0log(ζs0/ζs1)xF_s= _s^0 ( _s^0/ _s^1)\,dx under the integral sign and testing the continuity equations above against the functions log(ζs0/ζs1) ( _s^0/ _s^1) and ζs0/ζs1 _s^0/ _s^1, which are not compactly supported. Differentiating gives ∂sFs=∫∂sζs0logζs0ζs1x+∫∂sζs0x−∫ζs0ζs1∂sζs1x. _sF_s= _s _s^0\, _s^0 _s^1\,dx+ _s _s^0\,dx- _s^0 _s^1\, _s _s^1\,dx. (43) The middle term vanishes because ζs0 _s^0 integrates to one for all s. Substituting the continuity equations and integrating by parts, the first term becomes x∼ζs0[us0(x)⋅∇log(ζs0/ζs1)(x)]E_x _s^0[u_s^0(x)·∇ ( _s^0/ _s^1)(x)]. Similarly, the third term becomes −∫ζs1us1⋅∇(ζs0/ζs1)dx=−x∼ζs0[us1(x)⋅∇log(ζs0/ζs1)(x)]- _s^1\,u_s^1·∇( _s^0/ _s^1)\,dx=-E_x _s^0[u_s^1(x)·∇ ( _s^0/ _s^1)(x)], using ζs1∇(ζs0/ζs1)=ζs0∇log(ζs0/ζs1) _s^1∇( _s^0/ _s^1)= _s^0∇ ( _s^0/ _s^1). Hence ∂sFs=x∼ζs0[(us0(x)−us1(x))⋅∇logζs0ζs1(x)]. _sF_s=E_x _s^0 [ (u_s^0(x)-u_s^1(x) )·∇ _s^0 _s^1(x) ]. (44) By Lemma 3 (i), ∇log(ζs0/ζs1)=(us0−us1)/χs∇ ( _s^0/ _s^1)=(u_s^0-u_s^1)/ _s, so ∂sFs=‖us0−us1‖L2(ζs0)2/χs _sF_s=\|u_s^0-u_s^1\|_L^2( _s^0)^2/ _s. Exchanging the roles of the two channels gives ∂sGs=‖us0−us1‖L2(ζs1)2/χs _sG_s=\|u_s^0-u_s^1\|_L^2( _s^1)^2/ _s, and adding the two yields ∂s(Fs+Gs)=Dsvelχs≥0. _s(F_s+G_s)= D_s^vel _s≥ 0. (45) It remains to control the divergences near the pivot. Write ζs0 _s^0 and ζs1 _s^1 as mixtures of (asy,σs2d)N(a_sy, _s^2I_d) and (asy′,σs2d)N(a_sy , _s^2I_d), respectively, over the same mixing law (y,y′)∼ρ0⊗ρ1(y,y ) _0 _1. By joint convexity of the Kullback–Leibler divergence, Fs≤(y,y′)∼ρ0⊗ρ1[KL((asy,σs2d)∥(asy′,σs2d))]=as22σs2(y,y′)∼ρ0⊗ρ1[|y−y′|2],F_s _(y,y ) _0 _1 [KL (N(a_sy, _s^2I_d)\, \|\,N(a_sy , _s^2I_d) ) ]= a_s^22 _s^2\,E_(y,y ) _0 _1 [|y-y |^2 ], (46) and the same bound holds for GsG_s. The expectation on the right is finite under Assumption 1. Since as→0a_s→ 0 and σs2→12 _s^2→ 12 as s↓0s 0, both FsF_s and GsG_s tend to zero, and both are finite for every s<12s< 12. To conclude, let 0<ϵ<S<120<ε<S< 12. By Assumption 2, Fs+GsF_s+G_s is absolutely continuous on [ϵ,S][ε,S], so integrating equation 45 gives FS+GS−Fϵ−Gϵ=∫ϵSDsvelχss.F_S+G_S-F_ε-G_ε= _ε^S D_s^vel _s\,ds. (47) As ϵ↓0ε 0, the left-hand side converges to FS+GSF_S+G_S by equation 46. The right-hand side converges to SJ_S by monotone convergence, since the integrand is non-negative. This proves equation 41. ∎ The Jeffreys-weighted statistics thus accumulate the local path discrepancy into a global information divergence between the two channels. By Corollary 6, their population target on ℐτI_τ is 1/2−τ−τJ_1/2-τ-J_τ. Remark 8. Under the assumptions of Proposition 7, the divergence SJ_S is non-decreasing in S, and as S↑12S 12 it converges in [0,∞][0,∞] to the Jeffreys divergence KL(ρ0∥ρ1)+KL(ρ1∥ρ0)KL( _0\| _1)+KL( _1\| _0). Proof. Monotonicity is immediate, since SJ_S is the integral of a non-negative function over (0,S](0,S]. For the limit, it suffices by Proposition 7 to show that KL(ζS0∥ζS1)→KL(ρ0∥ρ1)KL( _S^0\| _S^1) ( _0\| _1), since the reverse divergence is treated identically. The densities ζS0 _S^0 and ζS1 _S^1 are the images of ρ0 _0 and ρ1 _1 under the same Markov kernel y↦(aSy,σS2d)y (a_Sy, _S^2I_d). The data-processing inequality therefore gives KL(ζS0∥ζS1)≤KL(ρ0∥ρ1)KL( _S^0\| _S^1) ( _0\| _1) for every S. As S↑12S 12, we have aS→1a_S→ 1 and σS→0 _S→ 0, so xtSℓ→xℓx_t_S → x_ in L2L^2 and hence in distribution. Since the Kullback–Leibler divergence is jointly lower semicontinuous with respect to weak convergence, lim infS↑1/2KL(ζS0∥ζS1)≥KL(ρ0∥ρ1) _S 1/2KL( _S^0\| _S^1) ( _0\| _1). Combining the two bounds gives convergence in [0,∞][0,∞]. ∎ 3.6 Held-Out Permutation Calibration We split each dataset into a training part and a held-out part, ℓ=ℓtr∪ℓteD_ =D_ ^tr _ ^te for ℓ∈0,1 ∈\0,1\. The fields η^ℓ η and u^ℓ u are fitted on ℓtrD_ ^tr alone, and all statistics are evaluated on te=(0te,1te)D^te=(D_0^te,D_1^te). Conditionally on the training data and on the randomness of training, the learned fields are fixed, and every held-out observation is mapped through the same scalar witness Φ . Writing Φ^i0=Φ^(x0(i)) _i^0= (x_0^(i)) and Φ^j1=Φ^(x1(j)) _j^1= (x_1^(j)) for the held-out scores, we use the generic statistic T^=1n0te∑i=1n0teΦ^i0−1n1te∑j=1n1teΦ^j1. T= 1n_0^te _i=1^n_0^te _i^0- 1n_1^te _j=1^n_1^te _j^1. (48) For each permutation r=1,…,Rr=1,…,R, the pooled scores are randomly relabelled while preserving the group sizes, yielding T^(r) T^(r). The two-sided permutation p-value is p^=1+∑r=1R|T^(r)|≥|T^|R+1, p= 1+ _r=1^R1 \| T^(r)|≥| T| \R+1, (49) and the test is Δ(0,1)=p^≤κ (D_0,D_1)=1\ p≤κ\. Proposition 9 (Finite-sample validity). Let G denote the training data together with all randomness used in training. Suppose that the held-out data are independent of G, that the Monte Carlo draws are i.i.d. across observations and independent of the data and of the group labels, and that the permutations are drawn independently of everything else. Then, under ℋ0H_0, for any fitted fields, ℙ[p^≤κ|]≤κ.P [ p≤κ\, |\,G ]≤κ. (50) In particular, Δ has level κ in the sense of equation 13. Proof. We first show that the held-out scores are conditionally exchangeable given G, then that the observed statistic is exchangeable with its permuted copies, and we conclude with a counting argument. Write nte=n0te+n1ten^te=n_0^te+n_1^te, and let W=(W1,…,Wnte)W=(W_1,…,W_n^te) be the pooled vector of held-out scores, with the first n0ten_0^te entries from 0teD_0^te. Each score is obtained by applying a G-measurable function to one held-out observation and its own Monte Carlo draws. Under ℋ0H_0, the held-out observations are i.i.d. from ρ0 _0 and independent of G, and the Monte Carlo draws are i.i.d. across observations and independent of the group labels. Hence, conditionally on G, the entries of W are i.i.d., and in particular W is exchangeable. No moment condition is needed here, since the learned fields are defined for any data. To pass from the scores to the statistics, write Wg=(Wg(1),…,Wg(nte))W_g=(W_g(1),…,W_g(n^te)) for a permutation g of 1,…,nte\1,…,n^te\, so that (Wg)g′=Wg∘g′(W_g)_g =W_g g . Let h(W)h(W) denote the difference of group means equation 48. Let g0g_0 be the identity, and let g1,…,gRg_1,…,g_R be the i.i.d. uniform permutations used by the test, so that T^(r)=h(Wgr) T^(r)=h(W_g_r). Draw an auxiliary uniform permutation g∗g_ , independent of everything else, and set V=Wg∗V=W_g_ . By exchangeability, V has the same conditional law as W and is independent of (g∗,g1,…,gR)(g_ ,g_1,…,g_R). Moreover, Wgr=Vg~rW_g_r=V_ g_r with g~0=g∗−1 g_0=g_ ^-1 and g~r=g∗−1∘gr g_r=g_ ^-1 g_r for r≥1r≥ 1. The map (g∗,g1,…,gR)↦(g~0,…,g~R)(g_ ,g_1,…,g_R) ( g_0,…, g_R) is a bijection, so g~0,…,g~R g_0,…, g_R are i.i.d. uniform and independent of V. Consequently, the vector (A0,…,AR)(A_0,…,A_R) with Ar=|h(Vg~r)|=|T^(r)|A_r=|h(V_ g_r)|=| T^(r)| is exchangeable conditionally on G, where T^(0)=T T^(0)= T. Finally, for j∈0,…,Rj∈\0,…,R\, define p^j=(R+1)−1#r:Ar≥Aj p_j=(R+1)^-1\#\r:A_r≥ A_j\, so that p^=p^0 p= p_0. Let J=j:p^j≤κJ=\j: p_j≤κ\. If J is non-empty, choose j∗∈Jj ∈ J minimising Aj∗A_j . Every j∈Jj∈ J then satisfies Aj≥Aj∗A_j≥ A_j , so |J|≤#r:Ar≥Aj∗≤κ(R+1)|J|≤\#\r:A_r≥ A_j \≤κ(R+1). By exchangeability, each p^j p_j has the same conditional law, and therefore ℙ[p^≤κ|]=1R+1∑j=0Rℙ[p^j≤κ|]=[|J||]R+1≤κ.P [ p≤κ\, |\,G ]= 1R+1 _j=0^RP [ p_j≤κ\, |\,G ]= E [|J|\, |\,G ]R+1≤κ. (51) This is equation 50. Taking expectations over G shows that the unconditional rejection probability is at most κ for every ρ0∈(ℝd) _0 (R^d), which is the level condition equation 13. ∎ Proposition 9 requires neither moment assumptions nor any accuracy of the learned fields. Consistency of the learned denoisers or velocities is required for power, through Corollary 6, but not for validity. 3.7 Consistency Proposition 9 holds for arbitrary fitted fields, so the quality of the learned fields enters only through power. To quantify this, we measure the error of each learned field under both channel densities, εℓ|k(s)=(x∼ζsk[|η^sℓ(x)−ηsℓ(x)|2])1/2,ℓ,k∈0,1, _ |k(s)= (E_x _s^k [| η_s (x)- _s (x)|^2 ] )^1/2, ,k∈\0,1\, (52) and we denote by D^srep D_s^rep the contrast equation 31 computed with the learned fields in place of the population ones. The diagonal errors ε0|0 _0|0 and ε1|1 _1|1 measure how well each field is learned on its own channel. The cross errors ε1|0 _1|0 and ε0|1 _0|1 measure how well each field extrapolates to the other channel, on which the witness also evaluates it. Lemma 10 (Risk contrasts with learned fields). Suppose Assumption 1 holds, fix s∈(0,12)s∈(0, 12), and let the learned denoisers be fixed and square integrable under ζs0 _s^0 and ζs1 _s^1. Then D^sden=(x∼ζs0[|ηs0(x)−η^s1(x)|2]−ε0|0(s)2)+(x∼ζs1[|ηs1(x)−η^s0(x)|2]−ε1|1(s)2). D_s^den= (E_x _s^0 [| _s^0(x)- η_s^1(x)|^2 ]- _0|0(s)^2 )+ (E_x _s^1 [| _s^1(x)- η_s^0(x)|^2 ]- _1|1(s)^2 ). (53) In particular, D^sden=0 D_s^den=0 under ℋ0H_0 for any learned fields, and |D^sden−Dsden|≤2Dsden(ε1|0(s)+ε0|1(s))+∑ℓ,k∈0,1εℓ|k(s)2. | D_s^den-D_s^den |≤ 2 D_s^den\, ( _1|0(s)+ _0|1(s) )+ _ ,k∈\0,1\ _ |k(s)^2. (54) The same statements hold for the velocity contrast, with u in place of η throughout. Each bracket in equation 53 is positive when, on the corresponding channel, the matching learned field is closer to the population field than the mismatched one. The test, therefore, retains power whenever each model fits its own population better than the other, which is considerably weaker than the consistency of either model. Proof of Lemma 10. The identity follows from the orthogonality argument of Lemma 5, now applied to the learned fields, and the bound follows by expanding the squares around the population fields. Fix s∈(0,12)s∈(0, 12), and abbreviate ∥⋅∥k=∥⋅∥L2(ζsk)\|·\|_k=\|·\|_L^2( _s^k) for k∈0,1k∈\0,1\. Let ω∼ρ0ω _0 and z be independent, let E denote the expectation over both, and write e=ω−ηs0(ωs)e=ω- _s^0( _s). As in the proof of Lemma 5, ωs∼ζs0 _s _s^0 and [e∣ωs]=0E[e _s]=0, and both e and ηs0(ωs) _s^0( _s) are square integrable, with second moments at most ω∼ρ0|ω|2E_ω _0|ω|^2. For any measurable f with ‖f‖0<∞\|f\|_0<∞, the variable (ηs0−f)(ωs)( _s^0-f)( _s) is therefore square integrable and a function of ωs _s. The tower property then gives |ω−f(ωs)|2=|e|2+‖ηs0−f‖02.E |ω-f( _s) |^2=E|e|^2+\| _s^0-f\|_0^2. (55) Applying this with f=η^s1f= η_s^1 and f=η^s0f= η_s^0, which are square integrable under ζs0 _s^0 by assumption, and subtracting, we obtain ω∼ρ0,z[ψsden(ω,z)]=‖ηs0−η^s1‖02−ε0|0(s)2.E_ω _0,\,z [ _s^den(ω,z) ]=\| _s^0- η_s^1\|_0^2- _0|0(s)^2. (56) The same argument for ω∼ρ1ω _1, with residual ω−ηs1(ωs)ω- _s^1( _s), gives ω∼ρ1,z[ψsden(ω,z)]=ε1|1(s)2−‖ηs1−η^s0‖12.E_ω _1,\,z [ _s^den(ω,z) ]= _1|1(s)^2-\| _s^1- η_s^0\|_1^2. (57) All terms are finite, and subtracting the second identity from the first yields equation 53. Under ℋ0H_0, Proposition 4 gives ζs0=ζs1 _s^0= _s^1 and ηs0=ηs1=:ηs _s^0= _s^1=: _s, so the two norms coincide and we write ∥⋅∥\|·\| for both. The identity equation 53 then reads D^sden=(‖ηs−η^s1‖2−‖ηs−η^s0‖2)+(‖ηs−η^s0‖2−‖ηs−η^s1‖2)=0. D_s^den= (\| _s- η_s^1\|^2-\| _s- η_s^0\|^2 )+ (\| _s- η_s^0\|^2-\| _s- η_s^1\|^2 )=0. (58) If any of the four errors is infinite, the right-hand side of equation 54 is infinite and there is nothing to prove, so we assume all four are finite. Write A=ηs0−ηs1A= _s^0- _s^1 and Bℓ=η^sℓ−ηsℓB_ = η_s - _s , so that ‖B1‖0=ε1|0(s)\|B_1\|_0= _1|0(s) and ‖B0‖1=ε0|1(s)\|B_0\|_1= _0|1(s). Then ηs0−η^s1=A−B1 _s^0- η_s^1=A-B_1 and ηs1−η^s0=−(A+B0) _s^1- η_s^0=-(A+B_0). The norms of A are finite, since ‖A‖0≤‖ηs0−η^s1‖0+ε1|0(s)\|A\|_0≤\| _s^0- η_s^1\|_0+ _1|0(s) and similarly for ‖A‖1\|A\|_1, and Lemma 5 gives Dsden=‖A‖02+‖A‖12D_s^den=\|A\|_0^2+\|A\|_1^2. Expanding the squares, ‖A−B1‖02−‖A‖02 \|A-B_1\|_0^2-\|A\|_0^2 =−2x∼ζs0[A(x)⋅B1(x)]+ε1|0(s)2, =-2\,E_x _s^0 [A(x)· B_1(x) ]+ _1|0(s)^2, (59) ‖A+B0‖12−‖A‖12 \|A+B_0\|_1^2-\|A\|_1^2 =2x∼ζs1[A(x)⋅B0(x)]+ε0|1(s)2. =2\,E_x _s^1 [A(x)· B_0(x) ]+ _0|1(s)^2. Subtracting DsdenD_s^den from equation 53 therefore gives D^sden−Dsden=−2x∼ζs0[A(x)⋅B1(x)]+2x∼ζs1[A(x)⋅B0(x)]+ε1|0(s)2+ε0|1(s)2−ε0|0(s)2−ε1|1(s)2. split D_s^den-D_s^den&=-2\,E_x _s^0 [A(x)· B_1(x) ]+2\,E_x _s^1 [A(x)· B_0(x) ]\\ & + _1|0(s)^2+ _0|1(s)^2- _0|0(s)^2- _1|1(s)^2. split (60) By the Cauchy–Schwarz inequality, the two cross terms are bounded in absolute value by 2‖A‖0ε1|0(s)2\|A\|_0\, _1|0(s) and 2‖A‖1ε0|1(s)2\|A\|_1\, _0|1(s). Since ‖A‖k≤Dsden\|A\|_k≤ D_s^den for both k, the triangle inequality yields equation 54. Velocity contrast: The same three steps apply with the target Us(ω,z)U_s(ω,z) in place of ω and the velocities in place of the denoisers. This uses that UsU_s is square integrable, that usℓ(ωs)=[Us(ω,z)∣ωs]u_s ( _s)=E[U_s(ω,z) _s] for ω∼ρℓω _ , and that us0=us1u_s^0=u_s^1 under ℋ0H_0 by Proposition 4(iv). ∎ Proposition 11 (Consistency). Suppose the sampling conditions of Proposition 9 hold, that Assumption 1 holds, that μ is strictly positive on ℐτI_τ, and that the integrability condition of Corollary 6 holds for the population fields. Suppose further that ℋ1H_1 holds with Twrep<∞T_w^rep<∞, and that the number of permutations satisfies R≥1/κ−1R≥ 1/κ-1. For ℓ∈0,1 ∈\0,1\, let Mℓ=[Φ^(ω)2|]M_ =E [ (ω)^2\, |\,G ] denote the conditional second moment of the score of a held-out observation ω∼ρℓω _ . Consider a sequence of experiments along which n0te,n1te→∞n_0^te,n_1^te→∞, and suppose that along this sequence (a) the errors equation 52 of the learned fields of representation reprep satisfy ∫ℐτw(s)εℓ|k(s)2s→0 _I_τw(s)\, _ |k(s)^2\,ds→ 0 in probability for all ℓ,k∈0,1 ,k∈\0,1\; (b) for every s∈ℐτs _τ, the learned fields are square integrable under ζs0 _s^0 and ζs1 _s^1, and M0M_0 and M1M_1 are bounded in probability. Then ℙ[Δ(0,1)=1]→1P [ (D_0,D_1)=1 ]→ 1. Condition (b) holds, for instance, if ρ0 _0 and ρ1 _1 have finite fourth moments, μ is bounded away from zero on ℐτI_τ, w is bounded on ℐτI_τ as are the weights equation 34, and the learned fields grow at most linearly, |η^sℓ(x)|≤C(1+|x|)| η_s (x)|≤ C(1+|x|), with C bounded in probability. The same holds for the velocity representation, with u^sℓ u_s in place of η^sℓ η_s . Condition (a) makes precise the role of the cross errors: near the pivot, the two channels overlap heavily, and extrapolation is mild, whereas near s=12s= 12 the channels separate and the cross errors are hardest to control. This motivates both the truncation τ and the uniform weighting, which, unlike the Jeffreys weights of the denoiser representation, does not up-weight the low-noise region. Proof of Section 3.7. We show that the observed statistic converges in probability to the positive population target, while every permuted statistic converges to zero. Within this proof, we drop the superscript on the held-out sizes, write n=n0+n1n=n_0+n_1, T=TwrepT=T_w^rep and M=max(M0,M1)M= (M_0,M_1), and let T~=[T^∣] T=E[ T ] denote the conditional mean of the statistic. Condition on G, so that the learned fields are fixed. For an observation ω∼ρℓω _ with Monte Carlo draws (sk,zk)k=1K(s_k,z_k)_k=1^K, the score is Φ^(ω)=K−1∑k=1KXk (ω)=K^-1 _k=1^KX_k with Xk=w(sk)ψskrep(ω,zk)/μ(sk)X_k=w(s_k)\, _s_k^rep(ω,z_k)/μ(s_k). Given ω and G, the XkX_k are i.i.d., so [Φ^(ω)2|ω,]=Var(X1∣ω,)K+[X1∣ω,]2≥[X12∣ω,]K,E [ (ω)^2\, |\,ω,G ]= Var(X_1 ω,G)K+E[X_1 ω,G]^2≥ E[X_1^2 ω,G]K, (61) and averaging over ω gives [X12∣]≤KMℓ<∞E[X_1^2 ]≤ KM_ <∞. In particular, X1X_1 is conditionally integrable. Since the learned fields are G-measurable and (ω,s1,z1)(ω,s_1,z_1) is independent of G, and since μ>0μ>0 on ℐτI_τ, Tonelli’s theorem gives [|X1||]=∫ℐτw(s)ω∼ρℓ,z[|ψsrep(ω,z)|]s,E [|X_1|\, |\,G ]= _I_τw(s)\,E_ω _ ,\,z [| _s^rep(ω,z)| ]\,ds, (62) with the learned fields inside ψsrep _s^rep held fixed. Fubini’s theorem then gives [Φ^(ω)∣]=∫ℐτw(s)ω∼ρℓ,z[ψsrep(ω,z)]sE[ (ω) ]= _I_τw(s)\,E_ω _ ,\,z[ _s^rep(ω,z)]\,ds, and by linearity, T~=∫ℐτw(s)D^sreps. T= _I_τw(s)\, D_s^rep\,ds. (63) Since T<∞T<∞ and w>0w>0, we have Dsrep<∞D_s^rep<∞ for almost every s∈ℐτs _τ. By condition (b), Lemma 10 applies at every s∈ℐτs _τ. Integrating equation 54 against w and applying the Cauchy–Schwarz inequality to the first term gives |T~−T|≤2T1/2(2∫ℐτw(s)(ε1|0(s)2+ε0|1(s)2)s)1/2+∑ℓ,k∈0,1∫ℐτw(s)εℓ|k(s)2s, | T-T |≤ 2\,T^1/2 (2 _I_τw(s) ( _1|0(s)^2+ _0|1(s)^2 )\,ds )^1/2+ _ ,k∈\0,1\ _I_τw(s)\, _ |k(s)^2\,ds, (64) where we used (ε1|0+ε0|1)2≤2(ε1|02+ε0|12)( _1|0+ _0|1)^2≤ 2( _1|0^2+ _0|1^2). By condition (a), the right-hand side converges to zero in probability, so T~→T T→ T in probability. Under ℋ1H_1, Corollary 6 gives T>0T>0. Conditionally on G, the held-out scores are independent, with the n0n_0 scores from 0teD_0^te identically distributed and likewise the n1n_1 scores from 1teD_1^te. Hence Var(T^∣)≤M0/n0+M1/n1≤M(1/n0+1/n1)Var( T )≤ M_0/n_0+M_1/n_1≤ M(1/n_0+1/n_1). For any δ>0δ>0 and C>0C>0, Chebyshev’s inequality gives ℙ[|T^−T~|>δ]≤ℙ[M>C]+Cδ2(1n0+1n1).P [| T- T|>δ ] [M>C]+ Cδ^2 ( 1n_0+ 1n_1 ). (65) The second term vanishes as n0,n1→∞n_0,n_1→∞, and the first can be made arbitrarily small by choosing C large, since M is bounded in probability. Hence T^−T~→0 T- T→ 0 in probability, and combined with Step 2, T^→T T→ T in probability. In particular, ℙ[|T^|≥T/2]→1P[| T|≥ T/2]→ 1. Let W be the pooled vector of held-out scores, with mean W¯ W and sample variance S2=(n−1)−1∑i=1n(Wi−W¯)2S^2=(n-1)^-1 _i=1^n(W_i- W)^2. For a uniform permutation g, independent of W, the permuted statistic is m0−m1m_0-m_1, where m0m_0 is the mean of a uniformly random subset of n0n_0 pooled scores and m1m_1 is the mean of its complement. Since n0m0+n1m1=nW¯n_0m_0+n_1m_1=n W, we have m0−m1=(n/n1)(m0−W¯)m_0-m_1=(n/n_1)(m_0- W). By the standard moments of sampling without replacement, [m0∣W]=W¯E[m_0 W]= W and Var(m0∣W)=n1S2/(nn0)Var(m_0 W)=n_1S^2/(n\,n_0). Hence, the permuted statistic has a conditional mean of zero and a conditional variance n2n12⋅n1S2nn0=(1n0+1n1)S2. n^2n_1^2· n_1S^2n\,n_0= ( 1n_0+ 1n_1 )S^2. (66) Moreover, S2≤n−1⋅1n∑i=1nWi2S^2≤ nn-1· 1n _i=1^nW_i^2, and the conditional mean of n−1∑iWi2n^-1 _iW_i^2 given G is at most M. By the conditional Markov inequality, S2S^2 is therefore bounded in probability. Chebyshev’s inequality, applied conditionally on W, gives for each r∈1,…,Rr∈\1,…,R\ ℙ[|T^(r)|>T/4|W]≤16T2(1n0+1n1)S2,P [| T^(r)|>T/4\, |\,W ]≤ 16T^2 ( 1n_0+ 1n_1 )S^2, (67) which converges to zero in probability. By bounded convergence, ℙ[|T^(r)|>T/4]→0P[| T^(r)|>T/4]→ 0, and a union bound over the fixed number R of permutations gives ℙ[max1≤r≤R|T^(r)|>T/4]→0P[ _1≤ r≤ R| T^(r)|>T/4]→ 0. Step 5: conclusion: By Steps 3 and 4, the event |T^|≥T/2∩maxr|T^(r)|≤T/4\| T|≥ T/2\∩\ _r| T^(r)|≤ T/4\ has probability tending to one. On this event, every indicator in equation 49 vanishes, so p^=1/(R+1)≤κ p=1/(R+1)≤κ by the choice R≥1/κ−1R≥ 1/κ-1, and the test rejects. ∎ Algorithm 1 Path-based two-sample test (PBI) 1: datasets 0,1D_0,D_1; representation rep∈den,velrep∈\den,vel\; path weight w (Table 1); truncation τ∈(0,14)τ∈(0, 14); sampling density μ>0μ>0 on ℐτI_τ; Monte Carlo draws K; permutations R≥1/κ−1R≥ 1/κ-1; level κ 2: decision Δ∈0,1 ∈\0,1\ and p-value p p 3: Split ℓ=ℓtr∪ℓteD_ =D_ ^tr _ ^te at random for ℓ∈0,1 ∈\0,1\ ⊳ independent of the data values 4: for ℓ∈0,1 ∈\0,1\ do 5: f^ℓ←FitField(ℓtr,rep) f ← FitField(D_ ^tr,rep) ⊳ Algorithm 2; training data only 6: end for 7: Pool the held-out data as (ω(i),ℓi)i=1N(ω^(i), _i)_i=1^N, with N=n0te+n1teN=n_0^te+n_1^te and group labels ℓi∈0,1 _i∈\0,1\ 8: for i=1,…,Ni=1,…,N do 9: Wi←Score(ω(i),f^0,f^1)W_i← Score(ω^(i); f^0, f^1) ⊳ Algorithm 3; label not used 10: end for 11: T^←MeanDiff(W,ℓ) T← MeanDiff(W, ) ⊳ statistic equation 48 12: for r=1,…,Rr=1,…,R do 13: ℓ(r)← ^(r)← uniformly random permutation of (ℓ1,…,ℓN)( _1,…, _N) 14: T^(r)←MeanDiff(W,ℓ(r)) T^(r)← MeanDiff(W, ^(r)) ⊳ reuses scores; no network evaluations 15: end for 16: p^←(1+∑r=1R|T^(r)|≥|T^|)/(R+1) p← (1+ _r=1^R1\| T^(r)|≥| T|\ )/(R+1) 17: return Δ=p^≤κ =1\ p≤κ\ and p p 18: 19: function MeanDiff(W,ℓW, ) 20: return 1|i:ℓi=0|∑i:ℓi=0Wi−1|i:ℓi=1|∑i:ℓi=1Wi 1|\i: _i=0\| _i: _i=0W_i\;-\; 1|\i: _i=1\| _i: _i=1W_i 21: end function Algorithm 2 FitField: fitting the field of one channel on its training split alone. No validation data or early stopping is used; the network after nitern_iter iterations is returned (Appendix A). 1: training data ℓtrD_ ^tr; representation reprep; network fθ(s,x)f_θ(s,x); training density μtr _tr on ℐτI_τ; iterations nitern_iter; batch size B 2: fitted field f^ℓ f , estimating ηsℓ _s if rep=denrep=den and usℓu_s if rep=velrep=vel 3: for iter=1,…,niteriter=1,…,n_iter do 4: Draw x(1),…,x(B)x^(1),…,x^(B) from ℓtrD_ ^tr, and independent sj∼μtrs_j _tr and zj∼(0,d)z_j (0,I_d) 5: for j=1,…,Bj=1,…,B do 6: xsj←asjx(j)+σsjzjx_s_j← a_s_jx^(j)+ _s_jz_j ⊳ Gaussian channel equation 18 7: yj←Target(x(j),zj,sj)y_j← Target(x^(j),z_j,s_j) 8: end for 9: Update θ by a gradient step on 1B∑j=1B|fθ(sj,xsj)−yj|2 1B _j=1^B |f_θ(s_j,x_s_j)-y_j |^2 10: end for 11: return f^ℓ=fθ f =f_θ 12: 13: function Target(x,z,sx,z,s) 14: if rep=denrep=den then 15: return x ⊳ denoiser target 16: else 17: return as′x+σs′za_s x+ _s z ⊳ outward velocity target 18: end if 19: end function Algorithm 3 Score: Monte Carlo estimate of the observation-level score. 1: held-out observation ω; fitted fields f^0,f^1 f^0, f^1; representation reprep; path weight w; sampling density μ on ℐτI_τ; Monte Carlo draws K 2: score Φ^wrep(ω) _w^rep(ω) 3: for k=1,…,Kk=1,…,K do 4: Draw sk∼μs_k μ and zk∼(0,d)z_k (0,I_d) ⊳ fresh draws for each observation 5: ωsk←askω+σskzk _s_k← a_s_kω+ _s_kz_k 6: yk←Target(ω,zk,sk)y_k← Target(ω,z_k,s_k) ⊳ Algorithm 2 7: ψk←|yk−f^1(sk,ωsk)|2−|yk−f^0(sk,ωsk)|2 _k← |y_k- f^1(s_k, _s_k) |^2- |y_k- f^0(s_k, _s_k) |^2 ⊳ local contrast 8: end for 9: return 1K∑k=1Kw(sk)μ(sk)ψk 1K _k=1^K w(s_k)μ(s_k)\, _k ⊳ importance weighting equation 70 Table 1: The four PBI variants, each combining a regression representation with a path weight and run with Algorithms 1, 2 and 3. The two Jeffreys weights give the same population target, the Jeffreys divergence of Proposition 7 up to path truncation. All weights are positive and bounded on IτI_τ. Variant Fitted fields Regression target y Path weight w(s)w(s) PBI-Denoiser_U η^s0,η^s1 η_s^0, η_s^1 ω 11 PBI-Denoiser_J η^s0,η^s1 η_s^0, η_s^1 ω λsas/σs2 _sa_s/ _s^2 PBI-Velocity_U u^s0,u^s1 u_s^0, u_s^1 as′ω+σs′za_s ω+ _s z 11 PBI-Velocity_J u^s0,u^s1 u_s^0, u_s^1 as′ω+σs′za_s ω+ _s z as/(λsσs2)a_s/( _s _s^2) 3.8 Relationship between the Four Variants The method yields four closely related test statistics, UniformJeffreys weightedDenoising riskT^UdenT^JdenVelocity riskT^UvelT^Jvel. array[]c|c&Uniform&Jeffreys weighted\\ risk& T_U^den& T_J^den\\[2.84526pt] Velocity risk& T_U^vel& T_J^vel. array (68) The rows differ in the regression representation used to probe the path, and the columns differ in how evidence is aggregated across radii. By Lemma 5, wJdenDsden=Dsvel/χsw_J^denD_s^den=D_s^vel/ _s at every radius, so the two Jeffreys variants share the population target identified in Proposition 7. Their finite-sample behaviour can nevertheless differ because the denoiser and the velocity regression have different noise and approximation properties. The uniform variants omit the information-theoretic reweighting. They provide a useful alternative when the Jeffreys weights amplify the difficult low-noise region near s=12s= 12, where σs→0 _s→ 0. 3.9 Monte Carlo Approximation over the Path In practice, neither the latent expectation nor the path integral in equation 33 is evaluated analytically. For each held-out observation ω(i)ω^(i), we draw independent sik∼μ,zik∼(0,d),k=1,…,K,s_ik μ, z_ik (0,I_d), k=1,…,K, (69) where μ is a sampling density that is strictly positive on ℐτI_τ, and use the importance-weighted estimate Φ^wrep(ω(i))=1K∑k=1Kw(sik)μ(sik)ψsikrep(ω(i),zik). _w^rep (ω^(i) )= 1K _k=1^K w(s_ik)μ(s_ik)\, _s_ik^rep (ω^(i),z_ik ). (70) Uniform sampling corresponds to constant μ, and the factor w/μw/μ permits nonuniform sampling in s without changing the population target equation 36. Increasing K reduces only the Monte Carlo component of the variance. The K corrupted copies of a single observation are therefore treated as one statistical unit and averaged before testing. 4 Experiments We evaluate the path-based inference (PBI) tests of Section 3 on one synthetic benchmark and three image benchmarks. Our study addresses two questions: (i) whether the proposed tests achieve higher power than existing two-sample tests as the sample size increases, and (i) how the regression representation (denoiser or velocity) and the path weighting (uniform or Jeffreys) affect performance across data modalities. We consider the four variants PBI-Denoiser_U, PBI-Denoiser_J, PBI-Velocity_U and PBI-Velocity_J (Table 1); the complete procedure is given in Algorithms 1, 2 and 3. The suffixes U and J denote the uniform and Jeffreys path weights of Section 3.5. Each variant is a valid test in its own right by Proposition 9, and we report all four to study the effects of these design choices. We compare against AutoTST [33], C2ST-L [11], C2ST-S [42], MMD-D and MMD-O [38], and MMDAgg [52]. Evaluation protocol: All methods are compared under the same total budget of 2n2n observations per population, irrespective of whether they require sample splitting. Methods that separate learning from testing (PBI, AutoTST, C2ST-L, C2ST-S, MMD-D and MMD-O) split each dataset into nℓtr=nℓte=n^tr_ =n^te_ =n observations; their learned components are fitted only on the training split and frozen before the test is evaluated on the independent held-out split. MMDAgg requires no training split and is instead evaluated on all 2n2n observations per population. Thus, n in all figures and tables indexes the split-based allocation, and every method receives 4n4n observations across the two populations, so the comparison accounts for both the data used to learn a statistic and the data used to evaluate it. For each of five seeds, the learned components are fitted once, and the test is then applied to 50 independently drawn held-out test sets of n observations per population. The power for a seed is the fraction of these 50 tests that reject ℋ0H_0 at level κ=0.05κ=0.05, and we report the mean over the five seeds (250 tests in total), quoting gains in percentage points (p) over the strongest baseline at each n. Implementation: Throughout this section, the two fields are fitted by separate networks, one per population, each trained on its own training split alone (the Separate architecture of Section 4.7); this applies to both the denoiser and the velocity representations. No validation set or early stopping is used, and we evaluate the final checkpoint after a fixed training budget. The training-time variable is restricted to s∈[0.01,0.49]s∈[0.01,0.49], and the time-embedding dimension is 64. Image inputs are normalised to [−1,1][-1,1], whereas HDGM samples are left unnormalised. Dataset-specific backbones and optimisation settings, together with the baseline implementations, are given in Appendix A. Figure 2: Test power at level κ=0.05κ=0.05 as a function of the per-population sample size n (ntr=nte=n_tr=n_te=n; every method receives 4n4n observations in total). MNIST, OrganAMNIST and Fashion-MNIST use class contamination at a fraction 0.10.1 (3→ 8, left→ lung, sneaker→ boot). Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds for the PBI variants. Highlighted baselines are those that are strongest at some sample size; the remaining baselines are shown in light grey. The HDGM axis is logarithmic. Reporting: Figure 2 summarises power on all four benchmarks. For each benchmark, a first table reports power at each sample size as the mean over five seeds ± the unbiased sample variance across seeds; the bands in Figure 2 instead show ±1± 1 standard deviation. A second table averages power over the sample sizes of the benchmark and compares each PBI variant with every baseline: each cell reports the absolute difference in percentage points and, in parentheses, the difference relative to that baseline’s power. Relative differences become large when a baseline’s power is close to the level κ, so they should be read alongside the absolute differences. 4.1 High-Dimensional Gaussian Mixtures Data: The high-dimensional Gaussian mixture (HDGM) benchmark of Liu et al. [38] consists of equal-weight, two-component Gaussian mixtures in d=10d=10 dimensions. Let c∼Bernoulli(12)c ( 12) denote the component indicator. Under both ρ0 _0 and ρ1 _1, the components have means 0 and 0.5 1d0.5\,1_d, for c=0c=0 and c=1c=1 respectively, and unit marginal variances. Under ρ0 _0, the coordinates are independent within each component. Under ρ1 _1, the first two coordinates are correlated, with correlation +0.5+0.5 in the first component and −0.5-0.5 in the second. Concretely, we draw v∼(,d)v (0,I_d), replace v2v_2 by ±12v1+32v2± 12v_1+ 32v_2 under ρ1 _1, with the sign ++ for c=0c=0 and −- for c=1c=1, and add the component mean 0.5c 1d0.5\,c\,1_d. The two densities therefore share their component weights, means and coordinate-wise marginal variances, and differ only in a localised, component-dependent covariance structure. We evaluate n∈500,1000,2000n∈\500,1000,2000\. Results: Figures 2, 2 and 3 show that the velocity representation is particularly effective on HDGM. At n=500n=500, PBI-Velocity_J reaches power 0.2400.240, exceeding the strongest baseline, C2ST-L, by 8.88.8 p. The advantage widens at n=1000n=1000, where PBI-Velocity_J attains 0.6840.684 against 0.4560.456 for C2ST-L, a gain of 22.822.8 p. At n=2000n=2000, PBI-Velocity_U performs best, reaching 0.9600.960 and exceeding the strongest baseline, AutoTST, by 13.213.2 p. The Jeffreys weighting markedly improves the denoiser, raising its power from 0.1120.112, 0.3080.308 and 0.7480.748 to 0.1800.180, 0.5000.500 and 0.9320.932 across the three sample sizes. As a result, PBI-Denoiser_J exceeds the strongest baseline at every sample size, whereas PBI-Denoiser_U remains below it. These results indicate that the velocity representation is well suited to the localised covariance shift in this benchmark, and that the Jeffreys weighting is beneficial in this low-dimensional setting. Table 2: Test power on HDGM (d=10d=10; localised covariance shift, Appendix A) at level κ=0.05κ=0.05 with n training and n test observations per population. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds (not the standard deviation). Bold marks the best method and underlining the best baseline in each column. Sample size n Method 500500 10001000 20002000 PBI-Denoiser_U 0.112±0.00650.112\, ± 0.0065 0.308±0.02010.308\, ± 0.0201 0.748±0.00650.748\, ± 0.0065 PBI-Denoiser_J 0.180±0.00980.180\, ± 0.0098 0.500±0.05520.500\, ± 0.0552 0.932±0.00230.932\, ± 0.0023 PBI-Velocity_U 0.200±0.01200.200\, ± 0.0120 0.656±0.02110.656\, ± 0.0211 0.960±0.00340.960\, ± 0.0034 PBI-Velocity_J 0.240±0.01100.240\, ± 0.0110 0.684±0.02310.684\, ± 0.0231 0.952±0.00310.952\, ± 0.0031 AutoTST 0.136±0.00450.136\, ± 0.0045 0.248±0.04330.248\, ± 0.0433 0.828¯±0.0117 0.828\, ± 0.0117 C2ST-L 0.152¯±0.0055 0.152\, ± 0.0055 0.456¯±0.0433 0.456\, ± 0.0433 0.552±0.01290.552\, ± 0.0129 C2ST-S 0.096±0.00350.096\, ± 0.0035 0.296±0.02870.296\, ± 0.0287 0.400±0.01580.400\, ± 0.0158 MMD-D 0.148±0.02090.148\, ± 0.0209 0.284±0.11750.284\, ± 0.1175 0.312±0.12630.312\, ± 0.1263 MMD-O 0.132±0.00450.132\, ± 0.0045 0.196±0.00030.196\, ± 0.0003 0.492±0.00210.492\, ± 0.0021 MMDAgg 0.120±0.00180.120\, ± 0.0018 0.296±0.00450.296\, ± 0.0045 0.696±0.00370.696\, ± 0.0037 Table 3: Power differences between each PBI variant (rows) and each baseline (columns) on HDGM, averaged over n∈500,1000,2000n∈\500,1000,2000\. With AmA_m the mean power of method m averaged over seeds and sample sizes, each cell reports the absolute difference 100(APBI−Ab)100(A_PBI-A_b) in percentage points, followed in parentheses by the relative difference 100(APBI−Ab)/Ab100(A_PBI-A_b)/A_b in percent. Positive values favour PBI, and the largest gain in each column is in bold. Baseline b PBI variant AutoTST C2ST-L C2ST-S MMD-D MMD-O MMDAgg PBI-Denoiser_U −1.5-1.5 (−3.6%)(-3.6\%) +0.3+0.3 (+0.7%)(+0.7\%) +12.5+12.5 (+47.5%)(+47.5\%) +14.1+14.1 (+57.0%)(+57.0\%) +11.6+11.6 (+42.4%)(+42.4\%) +1.9+1.9 (+5.0%)(+5.0\%) PBI-Denoiser_J +13.3+13.3 (+33.0%)(+33.0\%) +15.1+15.1 (+39.0%)(+39.0\%) +27.3+27.3 (+103.5%)(+103.5\%) +28.9+28.9 (+116.7%)(+116.7\%) +26.4+26.4 (+96.6%)(+96.6\%) +16.7+16.7 (+45.0%)(+45.0\%) PBI-Velocity_U +20.1+20.1 (+49.8%)(+49.8\%) +21.9+21.9 (+56.6%)(+56.6\%) +34.1+34.1 (+129.3%)(+129.3\%) +35.7+35.7 (+144.1%)(+144.1\%) +33.2+33.2 (+121.5%)(+121.5\%) +23.5+23.5 (+63.3%)(+63.3\%) PBI-Velocity_J +22.1+22.1 (+54.8%)(+54.8\%) +23.9+23.9 (+61.7%)(+61.7\%) +36.1+36.1 (+136.9%)(+136.9\%) +37.7+37.7 (+152.2%)(+152.2\%) +35.2+35.2 (+128.8%)(+128.8\%) +25.5+25.5 (+68.7%)(+68.7\%) 4.2 MNIST Data: MNIST contains greyscale images of handwritten digits [34]. We construct a controlled class-contamination alternative from a source digit and a target digit. Under ρ0 _0, all images are drawn from the source digit 33. Under ρ1 _1, a fraction 0.10.1 of the images are drawn from the target digit 88, and the remainder from the source digit. Each sample of n observations from ρ1 _1 contains exactly ⌊0.1n+12⌋ 0.1\,n+ 12 target-digit images, that is, 0.1n0.1\,n rounded half up, with the remaining images drawn independently from the source digit. We vary n∈50,100,150,180,200n∈\50,100,150,180,200\. Results: Tables 4 and 5 show that AutoTST is the strongest baseline at every sample size, and that the PBI variants become substantially more powerful as more data are available. At n=50n=50, PBI-Velocity_U ties AutoTST at power 0.1200.120, while the other variants remain slightly below it. At n=100n=100, PBI-Velocity_U reaches 0.2720.272, improving over AutoTST by 10.410.4 p (61.9%61.9\% relative), and both denoiser variants match AutoTST at 0.1680.168. From n=150n=150 onwards, all four PBI variants exceed the strongest baseline. The largest gain at n=180n=180 is obtained by PBI-Denoiser_U, which reaches 0.6040.604 and improves over AutoTST by 20.820.8 p. At n=200n=200, PBI-Velocity_U is strongest, with power 0.6680.668 and a gain of 25.225.2 p. The denoiser variants dip slightly between n=180n=180 and n=200n=200, a change within the seed-to-seed variability reported in Table 4. Overall, uniform weighting is more effective than Jeffreys weighting on this task, particularly for the velocity representation, which performs better under uniform weighting across all sample sizes. Table 4: Test power on MNIST with source digit 33 and target digit 88 at a contamination fraction of 0.10.1, at level κ=0.05κ=0.05 and with n training and n test observations per population. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds (not the standard deviation). The best result in each column is in bold, and the best baseline is underlined; tied entries share the formatting. Sample size n Method 5050 100100 150150 180180 200200 PBI-Denoiser_U 0.096±0.00330.096\, ± 0.0033 0.168±0.00650.168\, ± 0.0065 0.336±0.03390.336\, ± 0.0339 0.604±0.02270.604\, ± 0.0227 0.532±0.01810.532\, ± 0.0181 PBI-Denoiser_J 0.092±0.00090.092\, ± 0.0009 0.168±0.01290.168\, ± 0.0129 0.340±0.00920.340\, ± 0.0092 0.532±0.01570.532\, ± 0.0157 0.500±0.03060.500\, ± 0.0306 PBI-Velocity_U 0.120±0.00180.120\, ± 0.0018 0.272±0.01050.272\, ± 0.0105 0.376±0.01250.376\, ± 0.0125 0.544±0.05290.544\, ± 0.0529 0.668±0.03870.668\, ± 0.0387 PBI-Velocity_J 0.108±0.00230.108\, ± 0.0023 0.244±0.01690.244\, ± 0.0169 0.316±0.01950.316\, ± 0.0195 0.488±0.03850.488\, ± 0.0385 0.572±0.03930.572\, ± 0.0393 AutoTST 0.120¯±0.0048 0.120\, ± 0.0048 0.168¯±0.0117 0.168\, ± 0.0117 0.288¯±0.0071 0.288\, ± 0.0071 0.396¯±0.0245 0.396\, ± 0.0245 0.416¯±0.0361 0.416\, ± 0.0361 C2ST-L 0.052±0.00010.052\, ± 0.0001 0.044±0.00190.044\, ± 0.0019 0.112±0.00230.112\, ± 0.0023 0.104±0.00550.104\, ± 0.0055 0.120±0.00740.120\, ± 0.0074 C2ST-S 0.044±0.00030.044\, ± 0.0003 0.036±0.00110.036\, ± 0.0011 0.092±0.00410.092\, ± 0.0041 0.076±0.00250.076\, ± 0.0025 0.084±0.00430.084\, ± 0.0043 MMD-D 0.048±0.00090.048\, ± 0.0009 0.084±0.00210.084\, ± 0.0021 0.112±0.00350.112\, ± 0.0035 0.100±0.00020.100\, ± 0.0002 0.152±0.00170.152\, ± 0.0017 MMD-O 0.084±0.00110.084\, ± 0.0011 0.088±0.00090.088\, ± 0.0009 0.180±0.00160.180\, ± 0.0016 0.172±0.00010.172\, ± 0.0001 0.264±0.00610.264\, ± 0.0061 MMDAgg 0.100±0.00140.100\, ± 0.0014 0.136±0.00070.136\, ± 0.0007 0.240±0.01080.240\, ± 0.0108 0.248±0.00110.248\, ± 0.0011 0.308±0.00410.308\, ± 0.0041 Table 5: Power differences between each PBI variant (rows) and each baseline (columns) on MNIST, averaged over n∈50,100,150,180,200n∈\50,100,150,180,200\. Cells are computed as in Table 3. Positive values favour PBI, and the largest gain in each column is in bold. Baseline b PBI variant AutoTST C2ST-L C2ST-S MMD-D MMD-O MMDAgg PBI-Denoiser_U +7.0+7.0 (+25.1%)(+25.1\%) +26.1+26.1 (+301.9%)(+301.9\%) +28.1+28.1 (+422.9%)(+422.9\%) +24.8+24.8 (+250.0%)(+250.0\%) +19.0+19.0 (+120.3%)(+120.3\%) +14.1+14.1 (+68.2%)(+68.2\%) PBI-Denoiser_J +4.9+4.9 (+17.6%)(+17.6\%) +24.0+24.0 (+277.8%)(+277.8\%) +26.0+26.0 (+391.6%)(+391.6\%) +22.7+22.7 (+229.0%)(+229.0\%) +16.9+16.9 (+107.1%)(+107.1\%) +12.0+12.0 (+58.1%)(+58.1\%) PBI-Velocity_U +11.8+11.8 (+42.7%)(+42.7\%) +31.0+31.0 (+358.3%)(+358.3\%) +33.0+33.0 (+496.4%)(+496.4\%) +29.7+29.7 (+299.2%)(+299.2\%) +23.8+23.8 (+151.3%)(+151.3\%) +19.0+19.0 (+91.9%)(+91.9\%) PBI-Velocity_J +6.8+6.8 (+24.5%)(+24.5\%) +25.9+25.9 (+300.0%)(+300.0\%) +27.9+27.9 (+420.5%)(+420.5\%) +24.6+24.6 (+248.4%)(+248.4\%) +18.8+18.8 (+119.3%)(+119.3\%) +13.9+13.9 (+67.4%)(+67.4\%) 4.3 OrganAMNIST Data: OrganAMNIST is a MedMNIST subset consisting of grayscale axial CT slices derived from three-dimensional abdominal CT volumes and centred on annotated organ bounding boxes [61]. We use the 64×6464× 64 dimensionality and apply the class-contamination construction with Left Lung as the source class and Right Lung as the target class. Under ρ0 _0, all images are drawn from the source class, whereas each sample of n observations from ρ1 _1 contains ⌊0.1n+12⌋ 0.1\,n+ 12 target-class images, with the remainder drawn from the source class. We vary n∈100,200,350,420n∈\100,200,350,420\. Results: Tables 6 and 7 show that AutoTST is the strongest baseline at every sample size, yet every PBI variant exceeds it throughout the experiment. In the lower-sample regime, PBI-Denoiser_J is the strongest variant. At n=100n=100, it reaches 0.4880.488 power compared with 0.3120.312 for AutoTST, a gain of 17.617.6 p (56.4%56.4\% relative). At n=200n=200, the corresponding powers are 0.9680.968 and 0.6360.636, yielding a gain of 33.233.2 p (52.2%52.2\% relative). All four PBI variants attain power 1.0001.000 at n=350n=350 and n=420n=420, with zero variance across five training seeds, whereas AutoTST reaches 0.8960.896 and 0.9920.992, respectively. Averaged over all sample sizes, PBI-Denoiser_J has the largest gain against every baseline, including 15.515.5 p (21.9%21.9\% relative) over AutoTST and 25.125.1 p (40.9%40.9\% relative) over MMDAgg. Thus, OrganAMNIST favours the denoiser parameterisation in the lower-sample regime, while both parameterisations saturate once n≥350n≥ 350. Table 6: Test power on OrganAMNIST with source class Left Lung and target class Right Lung at a contamination fraction of 0.10.1, at level κ=0.05κ=0.05 and with n training and n test observations per population. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds (not the standard deviation). The best result in each column is in bold, and the best baseline is underlined. Sample size n Method 100100 200200 350350 420420 PBI-Denoiser_U 0.416±0.05430.416\, ± 0.0543 0.912±0.00030.912\, ± 0.0003 1.000±0.00001.000\, ± 0.0000 1.000±0.00001.000\, ± 0.0000 PBI-Denoiser_J 0.488±0.07370.488\, ± 0.0737 0.968±0.00130.968\, ± 0.0013 1.000±0.00001.000\, ± 0.0000 1.000±0.00001.000\, ± 0.0000 PBI-Velocity_U 0.388±0.00970.388\, ± 0.0097 0.900±0.00760.900\, ± 0.0076 1.000±0.00001.000\, ± 0.0000 1.000±0.00001.000\, ± 0.0000 PBI-Velocity_J 0.328±0.01410.328\, ± 0.0141 0.868±0.01170.868\, ± 0.0117 1.000±0.00001.000\, ± 0.0000 1.000±0.00001.000\, ± 0.0000 AutoTST 0.312¯±0.0291 0.312\, ± 0.0291 0.636¯±0.0199 0.636\, ± 0.0199 0.896¯±0.0123 0.896\, ± 0.0123 0.992¯±0.0003 0.992\, ± 0.0003 C2ST-L 0.072±0.00190.072\, ± 0.0019 0.124±0.00470.124\, ± 0.0047 0.228±0.00270.228\, ± 0.0027 0.264±0.00370.264\, ± 0.0037 C2ST-S 0.064±0.00070.064\, ± 0.0007 0.084±0.00110.084\, ± 0.0011 0.204±0.00430.204\, ± 0.0043 0.232±0.00450.232\, ± 0.0045 MMD-D 0.076±0.00110.076\, ± 0.0011 0.184±0.00830.184\, ± 0.0083 0.312±0.00690.312\, ± 0.0069 0.428±0.00530.428\, ± 0.0053 MMD-O 0.112±0.00290.112\, ± 0.0029 0.264±0.00110.264\, ± 0.0011 0.564±0.00050.564\, ± 0.0005 0.664±0.00150.664\, ± 0.0015 MMDAgg 0.192±0.00230.192\, ± 0.0023 0.440±0.01380.440\, ± 0.0138 0.860±0.00300.860\, ± 0.0030 0.960±0.00080.960\, ± 0.0008 Table 7: Power differences between each PBI variant (rows) and each baseline (columns) on OrganAMNIST, averaged over n∈100,200,350,420n∈\100,200,350,420\. Cells are computed as in Table 3. Positive values favour PBI, and the largest gain in each column is in bold. Baseline b PBI variant AutoTST C2ST-L C2ST-S MMD-D MMD-O MMDAgg PBI-Denoiser_U +12.3+12.3 (+17.3%)(+17.3\%) +66.0+66.0 (+383.7%)(+383.7\%) +68.6+68.6 (+469.9%)(+469.9\%) +58.2+58.2 (+232.8%)(+232.8\%) +43.1+43.1 (+107.5%)(+107.5\%) +21.9+21.9 (+35.7%)(+35.7\%) PBI-Denoiser_J +15.5+15.5 (+21.9%)(+21.9\%) +69.2+69.2 (+402.3%)(+402.3\%) +71.8+71.8 (+491.8%)(+491.8\%) +61.4+61.4 (+245.6%)(+245.6\%) +46.3+46.3 (+115.5%)(+115.5\%) +25.1+25.1 (+40.9%)(+40.9\%) PBI-Velocity_U +11.3+11.3 (+15.9%)(+15.9\%) +65.0+65.0 (+377.9%)(+377.9\%) +67.6+67.6 (+463.0%)(+463.0\%) +57.2+57.2 (+228.8%)(+228.8\%) +42.1+42.1 (+105.0%)(+105.0\%) +20.9+20.9 (+34.1%)(+34.1\%) PBI-Velocity_J +9.0+9.0 (+12.7%)(+12.7\%) +62.7+62.7 (+364.5%)(+364.5\%) +65.3+65.3 (+447.3%)(+447.3\%) +54.9+54.9 (+219.6%)(+219.6\%) +39.8+39.8 (+99.3%)(+99.3\%) +18.6+18.6 (+30.3%)(+30.3\%) 4.4 Fashion-MNIST Data: Fashion-MNIST contains greyscale images of apparel with the same dimensions and the same training and test splits as MNIST [60]. We apply the class-contamination construction, with Sneaker as the source class and Ankle Boot as the target class. Under ρ0 _0, all images are drawn from the source class, whereas each sample of n observations from ρ1 _1 contains ⌊0.1n+12⌋ 0.1\,n+ 12 target-class images, with the remainder drawn from the source class. We again vary n∈50,100,150,180,200n∈\50,100,150,180,200\. Results: Tables 8 and 9 reveal a different ordering of the PBI variants. MMDAgg is the strongest baseline at every sample size, and PBI-Denoiser_U is the only variant that exceeds it at every sample size. Its largest advantage occurs at n=100n=100, where it reaches 0.7320.732 against 0.5000.500 for MMDAgg, a gain of 23.223.2 p (46.4%46.4\% relative). It remains stronger at n=150n=150 and n=180n=180, with gains of 10.410.4 and 5.65.6 p respectively, and approaches ceiling power at n=200n=200 (0.9760.976 against 0.9640.964). PBI-Denoiser_J exceeds MMDAgg at n=100n=100 and is within 1.21.2 p of it at n=150n=150 and n=180n=180. Both velocity variants trail MMDAgg across all sample sizes and exhibit markedly greater variability across seeds, suggesting that velocity regression is less stable on this dataset. Thus, unlike HDGM and MNIST, Fashion-MNIST clearly favours the denoiser representation with uniform weighting. Table 8: Test power on Fashion-MNIST with source class Sneaker and target class Ankle Boot at a contamination fraction of 0.10.1, at level κ=0.05κ=0.05 and with n training and n test observations per population. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds (not the standard deviation). The best result in each column is in bold, and the best baseline is underlined. Sample size n Method 5050 100100 150150 180180 200200 PBI-Denoiser_U 0.272±0.00570.272\, ± 0.0057 0.732±0.00570.732\, ± 0.0057 0.916±0.00370.916\, ± 0.0037 0.952±0.00050.952\, ± 0.0005 0.976±0.00030.976\, ± 0.0003 PBI-Denoiser_J 0.172±0.00730.172\, ± 0.0073 0.608±0.06050.608\, ± 0.0605 0.800±0.00560.800\, ± 0.0056 0.892±0.01170.892\, ± 0.0117 0.936±0.00110.936\, ± 0.0011 PBI-Velocity_U 0.152±0.00610.152\, ± 0.0061 0.472±0.03550.472\, ± 0.0355 0.656±0.10150.656\, ± 0.1015 0.772±0.11770.772\, ± 0.1177 0.808±0.08410.808\, ± 0.0841 PBI-Velocity_J 0.132±0.00470.132\, ± 0.0047 0.360±0.07240.360\, ± 0.0724 0.516±0.13790.516\, ± 0.1379 0.700±0.12080.700\, ± 0.1208 0.712±0.13410.712\, ± 0.1341 AutoTST 0.132±0.00110.132\, ± 0.0011 0.392±0.02110.392\, ± 0.0211 0.408±0.01570.408\, ± 0.0157 0.512±0.01390.512\, ± 0.0139 0.624±0.01170.624\, ± 0.0117 C2ST-L 0.076±0.00150.076\, ± 0.0015 0.040±0.00140.040\, ± 0.0014 0.092±0.00230.092\, ± 0.0023 0.120±0.00380.120\, ± 0.0038 0.120±0.00080.120\, ± 0.0008 C2ST-S 0.036±0.00030.036\, ± 0.0003 0.020±0.00020.020\, ± 0.0002 0.080±0.00100.080\, ± 0.0010 0.104±0.00290.104\, ± 0.0029 0.092±0.00070.092\, ± 0.0007 MMD-D 0.064±0.00150.064\, ± 0.0015 0.128±0.00310.128\, ± 0.0031 0.136±0.00590.136\, ± 0.0059 0.136±0.00390.136\, ± 0.0039 0.172±0.00410.172\, ± 0.0041 MMD-O 0.100±0.00060.100\, ± 0.0006 0.208±0.00610.208\, ± 0.0061 0.380±0.00400.380\, ± 0.0040 0.384±0.00210.384\, ± 0.0021 0.512±0.00630.512\, ± 0.0063 MMDAgg 0.260¯±0.0030 0.260\, ± 0.0030 0.500¯±0.0004 0.500\, ± 0.0004 0.812¯±0.0041 0.812\, ± 0.0041 0.896¯±0.0049 0.896\, ± 0.0049 0.964¯±0.0015 0.964\, ± 0.0015 Table 9: Power differences between each PBI variant (rows) and each baseline (columns) on Fashion-MNIST, averaged over n∈50,100,150,180,200n∈\50,100,150,180,200\. Cells are computed as in Table 3. Positive values favour PBI, and the largest gain in each column is in bold. Baseline b PBI variant AutoTST C2ST-L C2ST-S MMD-D MMD-O MMDAgg PBI-Denoiser_U +35.6+35.6 (+86.1%)(+86.1\%) +68.0+68.0 (+758.9%)(+758.9\%) +70.3+70.3 (+1059.0%)(+1059.0\%) +64.2+64.2 (+505.0%)(+505.0\%) +45.3+45.3 (+142.9%)(+142.9\%) +8.3+8.3 (+12.1%)(+12.1\%) PBI-Denoiser_J +26.8+26.8 (+64.8%)(+64.8\%) +59.2+59.2 (+660.7%)(+660.7\%) +61.5+61.5 (+926.5%)(+926.5\%) +55.4+55.4 (+435.8%)(+435.8\%) +36.5+36.5 (+115.2%)(+115.2\%) −0.5-0.5 (−0.7%)(-0.7\%) PBI-Velocity_U +15.8+15.8 (+38.3%)(+38.3\%) +48.2+48.2 (+538.4%)(+538.4\%) +50.6+50.6 (+761.4%)(+761.4\%) +44.5+44.5 (+349.7%)(+349.7\%) +25.5+25.5 (+80.6%)(+80.6\%) −11.4-11.4 (−16.7%)(-16.7\%) PBI-Velocity_J +7.0+7.0 (+17.0%)(+17.0\%) +39.4+39.4 (+440.2%)(+440.2\%) +41.8+41.8 (+628.9%)(+628.9\%) +35.7+35.7 (+280.5%)(+280.5\%) +16.7+16.7 (+52.8%)(+52.8\%) −20.2-20.2 (−29.5%)(-29.5\%) 4.5 Type-I Error Figure 3 reports the empirical type-I error on all four benchmarks, with both samples drawn from ρ0 _0. The rejection rates of all four PBI variants remain close to the nominal level κ=0.05κ=0.05 at every sample size, with deviations comparable to the seed-to-seed variation shown by the shaded bands, as guaranteed by Proposition 9. Figure 3: Empirical Type-I error at level κ=0.05κ=0.05, with both samples drawn from ρ0 _0, as a function of the per-population sample size n (ntr=nte=n_tr=n_te=n). Lines show the mean rejection rate over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds for the PBI variants, and the dotted line marks κ. Baselines are styled as in Figure 2, and the HDGM axis is logarithmic. The PBI variants remain near the nominal level on all four benchmarks, as guaranteed by Proposition 9. 4.6 Samples from the Fitted Models Figures 4, 5 and 6 show samples generated by the separately trained denoiser and velocity models at the smallest sample size of each image benchmark, alongside held-out observations from both populations. The ρ0 _0 models are trained on the source class alone (digit 3, Left Lung and Sneaker), and the ρ1 _1 models on samples containing 10%10\% of the target class (digit 8, Right Lung and Ankle Boot). Although the samples are visibly noisy, the tests built from these fitted fields attain the powers reported in Sections 4.2, 4.3 and 4.4 at the same sample sizes. (a) Denoiser (b) Velocity Figure 4: Samples from the fitted MNIST models at ntr=nte=50n_tr=n_te=50 per population, for (a) the denoiser and (b) the velocity representation. Rows show, from top to bottom, held-out observations from ρ0 _0 (digit 3), samples from the ρ0 _0 model, held-out observations from ρ1 _1 (90% digit 3, 10% digit 8) and samples from the ρ1 _1 model; columns are unpaired examples, so a row need not contain a target-class image. Samples are generated with the probability flow ODE with 150150 steps. Although the samples are visibly noisy, models trained at this sample size already achieve the power in Table 4. (a) Denoiser (b) Velocity Figure 5: Samples from the fitted OrganAMNIST models at ntr=nte=100n_tr=n_te=100 per population, for (a) the denoiser and (b) the velocity representation. Rows show, from top to bottom, held-out observations from ρ0 _0 (Left Lung), samples from the ρ0 _0 model, held-out observations from ρ1 _1 (90% Left Lung, 10% Right Lung) and samples from the ρ1 _1 model; columns are unpaired examples, so a row need not contain a target-class image. Samples are generated with the probability flow ODE with 150150 steps. Although the samples are visibly noisy, models trained at this sample size already achieve the power in Table 6. (a) Denoiser (b) Velocity Figure 6: Samples from the fitted Fashion-MNIST models at ntr=nte=50n_tr=n_te=50 per population, for (a) the denoiser and (b) the velocity representation. Rows show, from top to bottom, held-out observations from ρ0 _0 (Sneaker), samples from the ρ0 _0 model, held-out observations from ρ1 _1 (90% Sneaker, 10% Ankle Boot) and samples from the ρ1 _1 model; columns are unpaired examples, so a row need not contain a target-class image. Samples are generated with the probability flow ODE with 150150 steps. Although the samples are visibly noisy, models trained at this sample size already achieve the power in Table 8. 4.7 Ablation on Model Architecture We compare four model architectures on MNIST and Fashion-MNIST: Conditioned, Separate, Common Contrast, and Restricted Common Contrast. All four architectures estimate the two fields f 0 f^\,0 and f 1 f^\,1, where f^ℓ(s,x) f^\, (s,x) estimates the denoiser ηsℓ(x)η _s(x) or the velocity usℓ(x)u _s(x), depending on the representation. They differ only in how parameters are shared between the two populations; the interpolant, the regression targets of Algorithm 2 and the test are unchanged. • Separate: two independent networks, f^ℓ(s,x)=fθℓ(s,x) f^\, (s,x)=f_ _ (s,x), each fitted on ℓtrD^tr_ alone. This is the architecture used in all other experiments. • Conditioned: a single network with a binary population input, f^ℓ(s,x)=fθ(s,x,ℓ) f^\, (s,x)=f_θ(s,x, ), which shares features between the populations while allowing the two fields to differ. • Common Contrast: a common field hθh_θ and a contrast field dϕd_φ, both of full capacity, with f 0=hθ−12dϕ f^\,0=h_θ- 12d_φ and f 1=hθ+12dϕ f^\,1=h_θ+ 12d_φ. This reparameterises rather than restricts the pair of fields, since any pair is recovered by h=(f0+f1)/2h=(f^0+f^1)/2 and d=f1−f0d=f^1-f^0; by Proposition 4, the population contrast vanishes at any single radius if and only if ℋ0H_0 holds. • Restricted Common Contrast: the same decomposition with a lower-capacity contrast network. The zero contrast of ℋ0H_0 remains representable, but under ℋ1H_1, the restricted dϕd_φ may capture only part of the population contrast, which can reduce variance at the cost of attenuating the signal. The shared architectures are fitted by minimising the sum of the two regression losses of Algorithm 2, each on its own training split. All four use the training splits only, so Proposition 9 applies unchanged: the architecture affects power but not validity. For the two common–contrast architectures, the local contrast reduces to ψs=−2dϕ(s,ωs)⋅(y−hθ(s,ωs)) _s=-2\,d_φ(s, _s)· (y-h_θ(s, _s) ), where y is the regression target, so the two-sample signal passes through dϕd_φ alone. All experiments use (ntrain=ntest=50)(n_train=n_test=50) observations per population and five training seeds (50 test sets each). We report mean empirical power and unbiased sample variance across seeds. We use ϵ=0.3ε=0.3 for the contamination, whereas the main experiment uses ϵ=0.1ε=0.1. Table 10: Architecture ablation on MNIST (source digit 3, target digit 8) at contamination fraction ϵ=0.3ε=0.3, higher than the 0.10.1 used in the main experiments, with ntr=nte=50n_tr=n_te=50 per population. Columns are the four architectures described in Section 4.7; Separate is the architecture used in all other experiments. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds; bold marks the best variant for each architecture. Method Conditioned Separate Common Contrast Restricted Common Contrast PBI-Denoiser_U 0.796±0.02890.796\, ± 0.0289 0.820±0.01140.820\, ± 0.0114 0.792±0.02970.792\, ± 0.0297 0.784±0.02470.784\, ± 0.0247 PBI-Denoiser_J 0.736±0.00750.736\, ± 0.0075 0.712±0.00210.712\, ± 0.0021 0.680±0.04220.680\, ± 0.0422 0.680±0.01500.680\, ± 0.0150 PBI-Velocity_U 0.784±0.05490.784\, ± 0.0549 0.768±0.02250.768\, ± 0.0225 0.744±0.00430.744\, ± 0.0043 0.668±0.08250.668\, ± 0.0825 PBI-Velocity_J 0.668±0.08010.668\, ± 0.0801 0.712±0.02670.712\, ± 0.0267 0.504±0.00870.504\, ± 0.0087 0.520±0.09880.520\, ± 0.0988 On MNIST, the Conditioned and Separate architectures achieve the highest mean power for all four PBI variants, as shown in Table 10. The strongest configuration is PBI-Denoiser_U with Separate models, which reaches power 0.8200.820, compared with 0.7960.796, 0.7920.792, and 0.7840.784 for Conditioned, Common Contrast, and Restricted Common Contrast, respectively. Architectural differences are more pronounced for PBI-Velocity_J: Separate models achieve 0.7120.712, exceeding Common Contrast and Restricted Common Contrast by 20.820.8 and 19.219.2 percentage points. Separate models also exhibit lower seed-to-seed variance than Conditioned models across all PBI variants. Table 11: Architecture ablation on Fashion-MNIST (source digit 3, target digit 8) at contamination fraction ϵ=0.3ε=0.3, higher than the 0.10.1 used in the main experiments, with ntr=nte=50n_tr=n_te=50 per population. Columns are the four architectures described in Section 4.7; Separate is the architecture used in all other experiments. Entries report the mean over five training seeds (50 test sets each; Appendix A) ± the unbiased sample variance across seeds; bold marks the best variant for each architecture. Method Conditioned Separate Common Contrast Restricted Common Contrast PBI-Denoiser_U 1.000±0.00001.000\, ± 0.0000 1.000±0.00001.000\, ± 0.0000 0.996±0.00010.996\, ± 0.0001 1.000±0.00001.000\, ± 0.0000 PBI-Denoiser_J 0.972±0.00110.972\, ± 0.0011 0.968±0.00370.968\, ± 0.0037 0.892±0.01070.892\, ± 0.0107 0.956±0.00310.956\, ± 0.0031 PBI-Velocity_U 0.912±0.00770.912\, ± 0.0077 0.868±0.01770.868\, ± 0.0177 0.896±0.00590.896\, ± 0.0059 0.952±0.00230.952\, ± 0.0023 PBI-Velocity_J 0.792±0.03770.792\, ± 0.0377 0.740±0.05900.740\, ± 0.0590 0.728±0.01810.728\, ± 0.0181 0.916±0.00430.916\, ± 0.0043 Table 11 shows ablation results on Fashion-MNIST. PBI-Denoiser_U achieves nearly saturated power across all architectures, ranging from 0.9960.996 to 1.0001.000, leaving little scope to distinguish their performance. Architectural choices have a larger effect on the velocity variants. Restricted Common Contrast achieves the highest mean power for both PBI-Velocity_U (0.9520.952) and PBI-Velocity_J (0.9160.916), improving over Separate models by 8.48.4 and 17.617.6 percentage points, respectively. It also yields the lowest sample variance among the four architectures for both velocity variants. Across both datasets, uniform weighting yields higher mean power than Jeffreys weighting for every matched architecture and parameterisation. PBI-Denoiser_U is the strongest variant within each architecture, whereas the preferred architecture depends on the dataset and parameterisation. These results indicate that architectural constraints can materially affect empirical power, particularly for the velocity variants. The comparisons are descriptive and do not establish statistical significance. 5 Observation-Level Asymmetry under Corruption The statistics of Section 3 are differences of group means of observation-level scores Φ^wrep(ω) ^rep_w(ω), so each held-out observation carries its own share of the evidence against ℋ0H_0. This appendix checks, on a controlled example, that these shares localise the discrepancy: when only part of ρ1 _1 departs from ρ0 _0, the observations responsible for the departure receive the largest scores. Setup: We reuse the two-dimensional construction of Figure 1. The reference distribution ρ0 _0 is an equal-weight mixture of three Gaussians with means (−1.6,−1)(-1.6,-1), (0,1.75)(0,1.75) and (1.6,−1)(1.6,-1) and covariance 0.120.1\,I_2. The second distribution is a contaminated copy, ρ1=0.9ρ0+0.1C, _1=0.9\, _0+0.1\,C, (71) where the corruption C is one of three structural perturbations of increasing difficulty: (i) a new cluster C=((0,−0.1), 0.12)C=N ((0,-0.1),\,0.1\,I_2 ) between the modes of ρ0 _0, where ρ0 _0 has almost no mass; (i) a shifted mode C=((2.2,−0.4), 0.12)C=N ((2.2,-0.4),\,0.1\,I_2 ), a copy of the right-hand mode translated by (0.6,0.6)(0.6,0.6) that partially overlaps it; (i) a deformed mode C=((0,1.75),diag(0.9, 0.02))C=N ((0,1.75),\,diag(0.9,\,0.02) ), which shares its centre with the top mode of ρ0 _0 but is stretched horizontally. Each held-out sample contains 500500 observations per group, so that about 5050 observations of 1teD_1^te are corrupted. The corruption labels are used only for evaluation and are never available during testing. Fields: We consider two sets of fields. The exact fields are the posterior-mean denoisers of ρ0 _0 and ρ1 _1, available in closed form for Gaussian mixtures. The fitted fields mimic the pipeline of Algorithm 1: for each group, we draw an independent training sample of 500500 observations, fit a Gaussian mixture by expectation–maximisation with the number of components selected by BIC over 1,…,6\1,…,6\, and use the posterior-mean denoiser of the fitted mixture. The fitted fields are therefore misspecified whenever the selected mixture does not resolve the corruption. Scores: For each held-out observation we compute the Monte Carlo score Φ^wden(ω) ^den_w(ω) of Equation 70 with K=64K=64 draws, truncation τ=0.05τ=0.05 and μ uniform on ℐτI_τ, for the uniform (wUw_U) and Jeffreys (wJdenw^den_J) weights, and report the asymmetry A(ω)=−Φ^wden(ω).A(ω)=- ^den_w(ω). (72) By Lemma 5, A is large when the channel-11 denoiser reconstructs ω much better than the channel-00 denoiser, and it has mean zero on the part of ρ1 _1 that coincides with ρ0 _0 up to reweighting. As a geometric counterpart of Figure 1, we also compute the mirror mismatch |X0−ω||X_0-ω| of an observation ω∼ρ1ω _1: we integrate the probability flow ODE from t=1−ϵt=1-ε to the bottleneck with the velocity of the ρ1 _1 half, continue to t=ϵt=ε with the velocity of the ρ0 _0 half, and record the endpoint X0X_0 (ϵ=2×10−3ε=2× 10^-3, adaptive Runge–Kutta with relative tolerance 10−610^-6). Under ℋ0H_0 the path is reflection symmetric and X0=ωX_0=ω (Proposition 4). Because it requires two ODE solves per observation, the mirror mismatch is computed on a subset of 3030 corrupted and 9090 uncorrupted observations per seed. Metrics: We rank the held-out observations from ρ1 _1 by A and report the area under the ROC curve (AUROC) for separating corrupted from uncorrupted observations, and the precision among the k highest-ranked observations (P@k), where k is the number of corrupted observations; a random ranking gives P@k=0.10k=0.10. We also report the rejection rate of the two-sided permutation test of Section 3.6 with the uniform denoiser statistic, R=999R=999 permutations and κ=0.05κ=0.05. All numbers are averaged over 1010 independent seeds. Table 12: Observation-level localisation under structural corruption, ρ1=0.9ρ0+0.1C _1=0.9\, _0+0.1\,C, with 500 held-out observations per group, averaged over 10 seeds. A¯ A is the mean asymmetry (denoiser, uniform weight) of observations from ρ0 _0 and of uncorrupted and corrupted observations from ρ1 _1. AUROC (uniform and Jeffreys weights) and P@k (uniform weight) measure how well A ranks the corrupted observations within the ρ1 _1 sample; random ranking gives AUROC 0.50.5 and P@k=0.10k=0.10. Mirror is the AUROC of the geometric mismatch |X0−ω||X_0-ω|, computed on 30 corrupted and 90 uncorrupted observations per seed. Rej. is the rejection rate of the test at κ=0.05κ=0.05 with R=999R=999 permutations. A¯ A AUROC Fields Corruption C ρ0 _0 ρ1 _1 matched ρ1 _1 corrupted U J P@k Mirror Rej. Exact new cluster −0.005-0.005 −0.005-0.005 0.1340.134 0.999 0.998 0.99 1.000 1.0 shifted mode −0.012-0.012 −0.008-0.008 0.2070.207 0.979 0.981 0.78 0.942 1.0 deformed mode −0.005-0.005 −0.006-0.006 0.1380.138 0.904 0.922 0.44 0.796 0.7 Fitted new cluster −0.006-0.006 −0.005-0.005 0.1390.139 0.982 1.000 0.89 0.990 0.9 shifted mode −0.010-0.010 −0.010-0.010 0.2040.204 0.982 0.981 0.79 0.925 0.9 deformed mode −0.006-0.006 −0.006-0.006 0.1340.134 0.828 0.867 0.38 0.750 0.4 Results: Table 12 shows that the asymmetry is concentrated on the corrupted observations. Observations from ρ0 _0 and uncorrupted observations from ρ1 _1 have a mean asymmetry close to zero, whereas corrupted observations score between 2020 and 4040 times higher, the observation-level counterpart of the broken reflection symmetry in the right panel of Figure 1. The new cluster and the shifted mode are ranked almost perfectly (AUROC between 0.980.98 and 1.001.00), and the corresponding tests reject in every or nearly every seed. The deformed mode is harder: its observations near the shared centre are individually indistinguishable from the uncorrupted top mode, so only its tails can be localised (P@k≈0.4k≈ 0.4), and the global test is also less powerful. Fitted fields reproduce these conclusions, with a loss concentrated on the deformed mode, where the selected mixture does not always resolve the corruption; this is the observation-level counterpart of the dependence of power on field accuracy in Proposition 11. The geometric mirror mismatch orders the corruptions in the same way. Still, it is less discriminative than A and requires ODE solves, thereby supporting the use of the regression-risk score for localisation. Scope: The scores rank observations but do not provide per-observation error control: the permutation p-value is a statement about the two distributions, and flagging individual observations with guarantees would require an additional multiple-testing or conformal procedure. We omit the velocity representation from Table 12 because, for fields of the affine form of Lemma 3, the velocity witness is a reweighting of the denoiser witness. Indeed, Us(ω,z)−usℓ(ωs)=as′ω+σs′z−λsηsℓ(ωs)−(σs′/σs)(asω+σsz)=λs(ω−ηsℓ(ωs))U_s(ω,z)-u _s( _s)=a _sω+σ _sz- _sη _s( _s)-(σ _s/ _s)(a_sω+ _sz)= _s (ω-η _s( _s) ), so ψsvel=λs2ψsdenψ^vel_s= _s^2\,ψ^den_s pointwise; in particular the Jeffreys-weighted velocity and denoiser scores coincide. The two representations differ in practice only because the velocity and denoiser networks are trained separately. 6 Related Work Two-sample tests with learned statistics: Kernel tests based on the maximum mean discrepancy (MMD) compare mean embeddings of the two distributions in a reproducing kernel Hilbert space [21], and their power depends strongly on the kernel and on the representation of the data. Most adaptive tests, therefore, learn part of the statistics from data, and they differ in how this learning is separated from calibration. Sample-splitting methods fit the statistic on one part of the data and evaluate it on the other, so that the learned statistic can be treated as fixed during calibration: kernels chosen to maximise a power criterion [22, 56], deep kernels [38], classifiers whose held-out accuracy or logits serve as the statistic [42, 11, 26], and real-valued witness functions whose mean difference is tested on held-out data [32], fitted with automated machine learning in AutoTST [33]. Other methods avoid the split by learning representations on auxiliary data [28], by accounting for kernel selection through selective inference [31], or by aggregating or fusing tests over a collection of kernels with jointly calibrated thresholds [52, 6]. Our tests follow the sample-splitting template, but the witness is not trained to discriminate between the samples: it is obtained from two arcs of a generative model, each fitted to a single population, and its population mean is characterised exactly in terms of the fields of a stochastic interpolant. Since the statistic is a difference of means of observation-level scores, it also localises the discrepancy in data space, which complements tests that learn interpretable test locations [14, 24]. Scores, denoisers and divergences: Score functions have long been used for model criticism, most prominently through the kernel Stein discrepancy for goodness-of-fit testing against a model with a known score [39, 13], and more recently for simple-versus-simple hypothesis testing and change-point detection when the scores of both hypotheses are known [45]. In our setting, neither distribution is known, and both fields are learned from finite samples. The link between Gaussian smoothing, score differences and the Kullback–Leibler divergence goes back to de Bruijn’s identity and its extension to relative entropy [43], and to the I-MMSE relation and its mismatched-estimation counterpart, which expresses the relative entropy between two Gaussian-smoothed distributions as an integral of excess denoising error [23, 57]. The same identities underlie likelihood bounds for weighted denoising score matching [55, 27] and information-theoretic views of diffusion models [29].Proposition 7 is the symmetrised form of the mismatched-estimation identity, written along the Gaussian-bottleneck path in terms of the two fields that the test estimates; its role here is to identify the population target of the Jeffreys-weighted statistics. Generative models as tools for inference: Beyond synthesis, diffusion and flow models have been used for likelihood evaluation through the probability flow ODE [54], for out-of-distribution detection [19, 41], and for zero-shot classification by comparing denoising errors across class-conditional models [36, 15]. The last line of work is closest in spirit to our denoising-risk contrast, which compares the held-out regression risks of two models. Still, it uses these comparisons for prediction rather than for calibrated inference. Paths of intermediate distributions have also been used to estimate density ratios, by telescoping over a sequence of bridges [49] or by integrating infinitesimal classifiers along a continuum [12]; our statistics likewise accumulate local evidence along a path, but we target a statistical test with finite-sample validity rather than a density-ratio estimate in this line of work. Finally, independently trained diffusion models and flow-based models learn remarkably similar fields [25, 62, 16, 9], which supports treating these fields as stable statistical objects. Midpoint inference from flow models: Closest to our work is a recent line of research built on the zero-flow criterion [58]: for rectified flow with independently coupled endpoints, the velocity field vanishes at t=12t= 12 if and only if the two endpoint distributions coincide. The criterion was introduced to certify conditional sufficiency in representation learning [58], has been turned into a midpoint divergence for training one-step generators, including a generalisation to symmetric stochastic interpolants [53], and, concurrent to our work, into a two-sample test that learns the midpoint velocity on a training split and evaluates its alignment with the displacements of held-out pairs, calibrated by a studentised asymptotic approximation [58]. These approaches and ours rest on the same symmetry: when the endpoints of a symmetric interpolant are exchangeable, its law is invariant under the time reflection t↦1−t 1-t. They exploit it differently. The zero-flow tests read the symmetry off a single field at a single time, which must be learned from pairs of observations drawn from both populations. For the Gaussian-bottleneck interpolant, the midpoint velocity vanishes identically under both hypotheses, and the evidence is instead spread across noise levels. Each half of the path is a Gaussian channel acting on one population, so the symmetry becomes an equality between two separately trained fields at every noise level (Proposition 4). Each field is, after rescaling its input, the denoiser of a standard diffusion model trained on a single population, which lets the test reuse diffusion training, yields observation-level scores, and admits permutation calibration with finite-sample validity. 7 Discussion and Limitations Figure 7: Observation-level localisation (ρ1=0.9ρ0+0.1C _1=0.9\, _0+0.1\,C, 500 held-out observations per group.) (a)–(c) Held-out ρ1 _1 observations coloured by the asymmetry A(ω)=−Φ^Uden(ω)A(ω)=- ^den_U(ω), computed with the exact denoisers for one seed; positive values mean the ρ1 _1 denoiser reconstructs ω better than the ρ0 _0 denoiser. Black edges mark corrupted observations, and dashed ellipses show the 2σ2σ region of C. For the deformed mode, uncorrupted points near the shared centre also score highly, so only the tails are localised. (d) AUROC for ranking the corrupted observations within the ρ1 _1 sample, averaged over 10 seeds, for A and for the geometric mirror mismatch |X0−ω||X_0-ω|, with exact (filled) and fitted (open) fields. Setup and full results: Sections 5 and 12. We showed that the generative dynamics learned by transport-based models are useful objects of statistical inference, not only tools for generation. Routing both distributions through a shared Gaussian bottleneck turns the null hypothesis into a reflection symmetry of the path, which held-out regression risks convert into observation-level scores. Permutation calibration makes the resulting tests valid in finite samples for any trained network. Across four benchmarks, they improve on the strongest baseline by up to 33 percentage points at an equal total sample budget. Sample splitting: The main cost of our approach is that the fields must be learned, which requires splitting each dataset into training and testing sets. This cost is shared by the learned baselines we compare against and by classical kernel tests when the bandwidth is selected from the data. A tempting way to avoid the split is to sample from the Gaussian midpoint and check the symmetry of the learned paths on generated data. Such a test is not calibrated: its type-I error is controlled only if the modelling error is negligible. Moreover, English et al. [17] shows that inferences drawn from a learned path agree across training and held-out data, but not with data newly sampled from the same model. Their observation that path lengths on real and generated data align when generation is accurate suggests a criterion for early stopping and model selection. Because it can be computed on the training split alone, it would leave the validity of our tests intact. Low-sample regimes and Model quality: Deep generative models might be expected to overfit when few observations are available. Recent work identifies several sources of implicit regularisation in diffusion and flow models, arising from the training dynamics, the architecture and the regression objective itself [7, 18, 5]. Consistent with this, our tests match or exceed the strongest baseline, even at the smallest sample size of every image benchmark. Overfitting cannot inflate the type-I error rate in any case, since calibration uses only held-out data; empirically, it stays close to κ (Figure 3). For modelling, we used off-the-shelf architectures and training settings without any hyperparameter search. Our goal was two-sample testing rather than generation, and we wanted to show that the tests remain powerful under modest computational budgets. Samples from our trained models (Figures 4, 5 and 6) show that a perfect generative model is not needed for a powerful test. Conversely, since power is governed by how accurately the fields are learned (Section 3.7), better-trained models, for instance, fine-tuned pretrained ones, should yield more powerful tests. Architecture also matters, mainly for the velocity variants (Section 4.7). Power, however, is limited to differences that the fitted fields can represent. If both models converge to the same misspecified approximation, the test has no power even though ℋ1H_1 holds. Choice of variant: No single combination of representation and weighting dominates across benchmarks, so the variant must be fixed in advance. Alternatively, several variants can be combined into a single valid test, for example, with a Bonferroni correction or by permuting the maximum of the standardised statistics under shared permutations. We have not yet evaluated such aggregated tests and leave them, together with the use of pretrained models, to future work. acknowledgments-disclosure-of-funding.We extend our gratitude to Prof. Masato Taki for helpful discussions. We gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 459422098, JST SPRING, Japan Grant Number JPMJSP2104, JST CREST (PMJCR2015), JSPS KAKENHI (24K02905, 25H01107), JST ERATO (JPMJER2601), the National Research Foundation, Singapore and the Ministry of Digital Development and Information under the AI Visiting Professorship Programme (award number AIVP-2024-004). Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views of the National Research Foundation, Singapore and the Ministry of Digital Development and Information. Appendix A Experimental Details and Notations Table 13 summarises the notation used throughout the paper. Table 13: Summary of notation. Symbol Meaning ρ0,ρ1 _0, _1 Unknown data densities; ℋ0:ρ0=ρ1H_0 _0= _1 ℓ=xℓ(i)i=1nℓD_ =\x_ ^(i)\_i=1^n_ Dataset drawn i.i.d. from ρℓ _ , ℓ∈0,1 ∈\0,1\ κ Significance level (often denoted α in the testing literature) xt,ρtx_t, _t Stochastic interpolant and its density at time t∈[0,1]t∈[0,1] z∼(0,d)z (0,I_d) Gaussian latent variable s=|t−12|s=|t- 12| Radial coordinate, the distance from the Gaussian bottleneck ts0,ts1t_s^0,t_s^1 Mirrored times 12−s 12-s and 12+s 12+s as,σsa_s, _s Radial signal and noise schedules ζsℓ _s Density of channel ℓ at radius s, ζsℓ=ρtsℓ _s = _t_s ; ζ0ℓ=(0,12d) _0 =N(0, 12I_d) and ζ1/2ℓ=ρℓ _1/2 = _ ηsℓ,usℓ _s ,u_s Population denoiser and outward velocity of channel ℓ η^sℓ,u^sℓ η_s , u_s Learned counterparts, fitted on ℓtrD_ ^tr λs,χs _s, _s Known positive scalars relating denoiser, velocity and score contrasts ω,ωsω, _s Held-out observation and its corrupted version at radius s ψsden,ψsvel _s^den, _s^vel Local risk contrasts at radius s DsrepD_s^rep Population signal at radius s, rep∈den,velrep∈\den,vel\ ℐτI_τ Truncated path interval [τ,12−τ][τ, 12-τ] wU,wJrepw_U,w_J^rep Uniform and Jeffreys path weights T^wrep T_w^rep Test statistic with representation reprep and weight w εℓ|k(s) _ |k(s) Error of learned field ℓ measured under channel k FC-MLP backbone: For vector-valued observations x∈ℝdx ^d, the network input is the concatenation [t,x][t,x] of the scalar time variable and the observation. The implementation supports appending a binary population indicator for a conditioned model, but the experiments reported here use separate networks for ρ0 _0 and ρ1 _1 and therefore use only [t,x][t,x]. The HDGM backbone has four hidden layers of width 128 with SiLU activations and a linear 1010 dimensional output. U-Net backbone: For image-valued observations, we use a three-level encoder–decoder U-Net. The encoder applies double-convolution blocks at channel widths c,2c,4c,8c,2c,4c,8c, with factor-two max pooling between levels. Conditional time embeddings are introduced at widths c, 2c2c, and 4c4c. The decoder uses transposed convolutions for upsampling and concatenates the corresponding encoder features, giving double-convolution inputs of widths 12c12c, 6c6c, and 3c3c. A final 1×11× 1 convolution maps the c-channel representation to the output. We set c=16c=16 for MNIST and Fashion-MNIST and c=32c=32 for OrganAMNIST. For HDGM, the fixed-budget implementation produces 1300, 1300, 2600 recorded optimisation steps for n=500n=500, 1000, and 2000, respectively. MNIST and OrganAMNIST use 1,600 at each sample size, while Fashion-MNIST uses 2,500 recorded steps. All timing records correspond to the final_fixed_budget checkpoint-selection rule. Baseline implementation details: We use the official NumPy implementation of MMDAgg [52], considering ten bandwidths and retaining its default Monte Carlo calibration parameters. AutoTST is implemented using the official codebase of Kübler et al. [33]; we use the AutoGluon tabular backend with the best_quality preset and a time budget of 150 seconds. Image observations are flattened before being passed to AutoTST. We adapt MMD-D, MMD-O, C2ST-L, and C2ST-S from the implementation of Liu et al. [38]. C2ST-S computes its statistics from the classifier’s hard class predictions, whereas C2ST-L uses its predicted class probabilities. MMD-O optimises the bandwidth of a Gaussian kernel applied directly to the observations, while MMD-D jointly learns a neural feature representation and the parameters of a deep kernel by maximising a variance-normalised empirical MMD objective. We use convolutional feature networks for image data and fully connected networks for vector-valued data. Appendix B Per-Benchmark Power and Type-I Error Curves For better readability, this appendix shows the power curves of Figure 2 and the type-I error curves of Figure 3 for each benchmark separately, with all ten methods shown individually. Figures 8, 10, 12 and 14 report power, and Figures 9, 11, 13 and 15 report the empirical type-I error at level κ=0.05κ=0.05; the numerical power values are in the tables of Section 4. Figure 8: Test power on HDGM at level κ=0.05κ=0.05 against the per-population sample size n (ntr=nte=n_tr=n_te=n), with all ten methods shown individually; a condensed version appears in Figure 2. Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; numerical values are in Table 2. Figure 9: Empirical Type-I error on HDGM at level κ=0.05κ=0.05, with both samples drawn from ρ0 _0, against the per-population sample size n (ntr=nte=n_tr=n_te=n). Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; the dotted line marks κ. The PBI variants stay close to the nominal level, as guaranteed by Proposition 9. Figure 10: Test power on MNIST at level κ=0.05κ=0.05 against the per-population sample size n (ntr=nte=n_tr=n_te=n), with all ten methods shown individually; a condensed version appears in Figure 2. Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; numerical values are in Table 4. Figure 11: Empirical Type-I error on MNIST at level κ=0.05κ=0.05, with both samples drawn from ρ0 _0, against the per-population sample size n (ntr=nte=n_tr=n_te=n). Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; the dotted line marks κ. The PBI variants stay close to the nominal level, as guaranteed by Proposition 9. Figure 12: Test power on OrganAMNIST at level κ=0.05κ=0.05 against the per-population sample size n (ntr=nte=n_tr=n_te=n), with all ten methods shown individually; a condensed version appears in Figure 2. Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; numerical values are in Table 6. Figure 13: Empirical Type-I error on OrganAMNIST at level κ=0.05κ=0.05, with both samples drawn from ρ0 _0, against the per-population sample size n (ntr=nte=n_tr=n_te=n). Lines show the mean over five training seeds (50 test sets each); shaded bands show ±1± 1 standard deviation across seeds; the dotted line marks κ. The PBI variants stay close to the nominal level, as guaranteed by Proposition 9. Figure 14: Test power on Fashion-MNIST at level κ=0.05κ=0.05 against the per-population sample size n (ntr=nte=n_tr=n_te=n), with all ten methods shown individually; a condensed version appears in Figure 2. Lines show the mean over five training seeds (50 test sets each) and shaded bands ±1± 1 standard deviation across seeds; numerical values are in Table 8. Figure 15: Empirical Type-I error on Fashion-MNIST at level κ=0.05κ=0.05, with both samples drawn from ρ0 _0, against the per-population sample size n (ntr=nte=n_tr=n_te=n). Lines show the mean over five training seeds (50 test sets each) and shaded bands ±1± 1 standard deviation across seeds; the dotted line marks κ. The PBI variants stay close to the nominal level, as guaranteed by Proposition 9. References [1] M. Albergo, N. M. Boffi, and E. Vanden-Eijnden. Stochastic Interpolants: A Unifying Framework for Flows and Diffusions. Journal of Machine Learning Research, 26(209):1–80, 2025. ISSN 1533-7928. URL http://jmlr.org/papers/v26/23-1605.html. [2] M. S. Albergo and E. Vanden-Eijnden. Building Normalizing Flows with Stochastic Interpolants. Sept. 2022. URL https://openreview.net/forum?id=li7qeBbCR1t. [3] T. Anderson. Hotelling (1931) the generalization of. Breakthroughs in Statistics, 1:45, 1991. [4] S. Barocas, M. Hardt, and A. Narayanan. Fairness and Machine Learning: Limitations and Opportunities. MIT Press, Dec. 2023. ISBN 9780262376525. Google-Books-ID: HuGwEAAAQBAJ. [5] Q. Bertrand, A. Gagneux, M. Massias, and R. Emonet. On the Closed-Form of Flow Matching: Generalization Does Not Arise from Target Stochasticity, Dec. 2025. URL http://arxiv.org/abs/2506.03719. arXiv:2506.03719 [cs]. [6] F. Biggs, A. Schrab, and A. Gretton. MMD-Fuse: Learning and Combining Kernels for Two-Sample Testing Without Data Splitting. Nov. 2023. URL https://openreview.net/forum?id=JOkgEY9os2. [7] T. Bonnaire, R. Urfin, G. Biroli, and M. Mezard. Why Diffusion Models Don’t Memorize: The Role of Implicit Dynamical Regularization in Training. Oct. 2025. URL https://openreview.net/forum?id=BSZqpqgqM0. [8] V. D. Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, Aug. 2022. ISSN 2835-8856. URL https://openreview.net/forum?id=MhK5aXo3gB. [9] R. Briq, M. Kamp, O. Fried, S. Cohen, and S. Kesselheim. The Amazing Stability of Flow Matching. Nov. 2025. URL https://openreview.net/forum?id=K6ptzIhoKy. [10] C. Bunne, Y.-P. Hsieh, M. Cuturi, and A. Krause. The Schrödinger Bridge between Gaussian Measures has a Closed Form. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, pages 5802–5833. PMLR, Apr. 2023. URL https://proceedings.mlr.press/v206/bunne23a.html. [11] X. Cheng and A. Cloninger. Classification logit two-sample testing by neural networks for differentiating near manifold densities. IEEE Transactions on Information Theory, 68(10):6631–6662, 2022. doi: 10.1109/TIT.2022.3175691. [12] K. Choi, C. Meng, Y. Song, and S. Ermon. Density Ratio Estimation via Infinitesimal Classification, Mar. 2022. URL http://arxiv.org/abs/2111.11010. arXiv:2111.11010 [cs.LG]. [13] K. Chwialkowski, H. Strathmann, and A. Gretton. A Kernel Test of Goodness of Fit. In Proceedings of The 33rd International Conference on Machine Learning, pages 2606–2615. PMLR, June 2016. URL https://proceedings.mlr.press/v48/chwialkowski16.html. [14] K. P. Chwialkowski, A. Ramdas, D. Sejdinovic, and A. Gretton. Fast Two-Sample Testing with Analytic Representations of Probability Measures. In Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015. URL https://papers.nips.c/paper_files/paper/2015/hash/b571ecea16a9824023e1af16897a582-Abstract.html. [15] K. Clark and P. Jaini. Text-to-Image Diffusion Models are Zero-Shot Classifiers, Sept. 2023. URL http://arxiv.org/abs/2303.15233. arXiv:2303.15233 [cs.CV]. [16] E. English and T. Suzuki. Path Invariance and the Robustness of Flow Matching: Beyond Architectural and Data Perturbations. Mar. 2026. URL https://openreview.net/forum?id=7XH9TUEzCI. [17] E. English, W. Huang, and T. Suzuki. Is your Flow Matching Model Really Generalising? A Path-Length Diagnostic. May 2026. URL https://openreview.net/forum?id=mG9tp027lG. [18] T. Farghly, P. Rebeschini, G. Deligiannidis, and A. Doucet. Implicit Regularisation in Diffusion Models: An Algorithm-Dependent Generalisation Analysis. Oct. 2025. URL https://openreview.net/forum?id=3lSqgESdPu. [19] M. S. Graham, W. H. L. Pinaya, P.-D. Tudosiu, P. Nachev, S. Ourselin, and M. J. Cardoso. Denoising diffusion models for out-of-distribution detection, Apr. 2023. URL http://arxiv.org/abs/2211.07740. arXiv:2211.07740 [cs.LG]. [20] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud. FFJORD: Free-Form Continuous Dynamics for Scalable Reversible Generative Models. Sept. 2018. URL https://openreview.net/forum?id=rJxgknCcK7. [21] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola. A Kernel Two-Sample Test. Journal of Machine Learning Research, 13(25):723–773, 2012a. ISSN 1533-7928. URL http://jmlr.org/papers/v13/gretton12a.html. [22] A. Gretton, D. Sejdinovic, H. Strathmann, S. Balakrishnan, M. Pontil, K. Fukumizu, and B. K. Sriperumbudur. Optimal kernel choice for large-scale two-sample tests. In Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012b. URL https://proceedings.neurips.c/paper/2012/hash/dbe272bab69f8e13f14b405e038deb64-Abstract.html. [23] D. Guo, S. Shamai, and S. Verdu. Mutual Information and Minimum Mean-square Error in Gaussian Channels, Dec. 2004. URL http://arxiv.org/abs/cs/0412108. arXiv:cs/0412108. [24] W. Jitkrittum, Z. Szabó, K. P. Chwialkowski, and A. Gretton. Interpretable Distribution Features with Maximum Testing Power. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016. URL https://papers.nips.c/paper_files/paper/2016/hash/0a09c8844ba8f0936c20bd791130d6b6-Abstract.html. [25] Z. Kadkhodaie, F. Guth, E. Simoncelli, and S. Mallat. Generalization in diffusion models arises from geometry-adaptive harmonic representations. International Conference on Learning Representations, 2024:46543–46567, May 2024. URL https://proceedings.iclr.c/paper_files/paper/2024/hash/cbaf319a4712385b5ba8a414808b5713-Abstract-Conference.html. [26] I. Kim, A. Ramdas, A. Singh, and L. Wasserman. Classification accuracy as a proxy for two sample testing, Feb. 2020. URL http://arxiv.org/abs/1602.02210. arXiv:1602.02210 [cs.LG]. [27] D. P. Kingma, T. Salimans, B. Poole, and J. Ho. Variational Diffusion Models. Nov. 2021. URL https://openreview.net/forum?id=2LdBqxc1Yv. [28] M. Kirchler, S. Khorasani, M. Kloft, and C. Lippert. Two-sample Testing Using Deep Learning. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pages 1387–1398. PMLR, June 2020. URL https://proceedings.mlr.press/v108/kirchler20a.html. [29] X. Kong, R. Brekelmans, and G. V. Steeg. Information-Theoretic Diffusion, Feb. 2023. URL http://arxiv.org/abs/2302.03792. arXiv:2302.03792 [cs.LG]. [30] Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catanzaro. DiffWave: A Versatile Diffusion Model for Audio Synthesis. Oct. 2020. URL https://openreview.net/forum?id=a-xFK8Ymz5J. [31] J. Kübler, W. Jitkrittum, B. Schölkopf, and K. Muandet. Learning Kernel Tests Without Data Splitting. In Advances in Neural Information Processing Systems, volume 33, pages 6245–6255. Curran Associates, Inc., 2020. URL https://proceedings.neurips.c/paper/2020/hash/44f683a84163b3523afe57c2e008bc8c-Abstract.html. [32] J. M. Kübler, W. Jitkrittum, B. Schölkopf, and K. Muandet. A Witness Two-Sample Test, Feb. 2022a. URL http://arxiv.org/abs/2102.05573. arXiv:2102.05573 [cs.LG]. [33] J. M. Kübler, V. Stimper, S. Buchholz, K. Muandet, and B. Schölkopf. AutoML Two-Sample Test. Oct. 2022b. URL https://openreview.net/forum?id=JLWOTZpWZzY. [34] Y. LeCun, C. Cortes, and C. J. Burges. The mnist database of handwritten digits. http://yann.lecun.com/exdb/mnist/, 1998. [35] E. L. Lehmann and J. P. Romano. Invariance. In Testing Statistical Hypotheses, pages 212–276. Springer, New York, NY, 2005. ISBN 9780387276052. doi: 10.1007/0-387-27605-X_6. URL https://doi.org/10.1007/0-387-27605-X_6. [36] A. C. Li, M. Prabhudesai, S. Duggal, E. Brown, and D. Pathak. Your Diffusion Model is Secretly a Zero-Shot Classifier, Sept. 2023. URL http://arxiv.org/abs/2303.16203. arXiv:2303.16203 [cs.LG]. [37] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le. Flow Matching for Generative Modeling. Sept. 2022. URL https://openreview.net/forum?id=PqvMRDCJT9t. [38] F. Liu, W. Xu, J. Lu, G. Zhang, A. Gretton, and D. J. Sutherland. Learning deep kernels for non-parametric two-sample tests. In H. D. I and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 6316–6326. PMLR, 13–18 Jul 2020. URL https://proceedings.mlr.press/v119/liu20m.html. [39] Q. Liu, J. Lee, and M. Jordan. A Kernelized Stein Discrepancy for Goodness-of-fit Tests. In Proceedings of The 33rd International Conference on Machine Learning, pages 276–284. PMLR, June 2016. URL https://proceedings.mlr.press/v48/liub16.html. [40] X. Liu, C. Gong, and Q. Liu. Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow. Sept. 2022. URL https://openreview.net/forum?id=XVjTT1nw5z. [41] Z. Liu, J. P. Zhou, Y. Wang, and K. Q. Weinberger. Unsupervised Out-of-Distribution Detection with Diffusion Inpainting, Aug. 2023. URL http://arxiv.org/abs/2302.10326. arXiv:2302.10326 [cs.CV]. [42] D. Lopez-Paz and M. Oquab. Revisiting Classifier Two-Sample Tests. Feb. 2017. URL https://openreview.net/forum?id=SJkXfE5x. [43] S. Lyu. Interpretation and Generalization of Score Matching, May 2012. URL http://arxiv.org/abs/1205.2629. arXiv:1205.2629 [cs.LG]. [44] B. Mlodozeniec, R. Eschenhagen, J. Bae, A. Immer, D. Krueger, and R. Turner. Influence Functions for Scalable Data Attribution in Diffusion Models, May 2025. URL http://arxiv.org/abs/2410.13850. arXiv:2410.13850 [cs.LG]. [45] S. Moushegian, T. Banerjee, and V. Tarokh. Diffusion-Based Hypothesis Testing and Change-Point Detection, June 2025. URL http://arxiv.org/abs/2506.16089. arXiv:2506.16089 [stat.ML]. [46] J. Pidstrigach. Score-Based Generative Models Detect Manifolds. Oct. 2022. URL https://openreview.net/forum?id=AiNrnIrDfD9. [47] S. Rabanser, S. Günnemann, and Z. Lipton. Failing Loudly: An Empirical Study of Methods for Detecting Dataset Shift. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL https://proceedings.neurips.c/paper_files/paper/2019/hash/846c260d715e5b854ffad5f70a516c88-Abstract.html. [48] A. Ramdas, S. J. Reddi, B. Poczos, A. Singh, and L. Wasserman. On the Decreasing Power of Kernel and Distance Based Nonparametric Hypothesis Tests in High Dimensions. Proceedings of the AAAI Conference on Artificial Intelligence, 29(1), Mar. 2015. ISSN 2374-3468. doi: 10.1609/aaai.v29i1.9692. URL https://ojs.aaai.org/index.php/AAAI/article/view/9692. [49] B. Rhodes, K. Xu, and M. Gutmann. Telescoping Density-Ratio Estimation. In Advances in Neural Information Processing Systems, volume 33, pages 4905–4916. Curran Associates, Inc., 2020. URL https://proceedings.neurips.c/paper/2020/hash/33d3b157ddc0896addfb22fa2a519097-Abstract.html. [50] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer. High-Resolution Image Synthesis With Latent Diffusion Models. pages 10684–10695, 2022. URL https://openaccess.thecvf.com/content/CVPR2022/html/Rombach_High-Resolution_Image_Synthesis_With_Latent_Diffusion_Models_CVPR_2022_paper.html. [51] T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, X. Chen, and X. Chen. Improved Techniques for Training GANs. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016. URL https://proceedings.neurips.c/paper_files/paper/2016/hash/8a3363abe792db2d8761d6403605aeb7-Abstract.html. [52] A. Schrab, I. Kim, M. Albert, B. Laurent, B. Guedj, and A. Gretton. Mmd aggregated two-sample test. Journal of Machine Learning Research, 24(194):1–81, 2023. URL http://jmlr.org/papers/v24/21-1289.html. [53] D. Shlenskii, N. Gushchin, L. Novitskiy, D. V. Dylov, and A. Korotin. Midpoint Generative Models, June 2026. URL http://arxiv.org/abs/2605.29920. arXiv:2605.29920 [cs.LG]. [54] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-Based Generative Modeling through Stochastic Differential Equations. Oct. 2020. URL https://openreview.net/forum?id=PxTIG12RRHS&utm_campaign=NLP%20News&utm_medium=email&utm_source=Revue%20newsletter. [55] Y. Song, C. Durkan, I. Murray, and S. Ermon. Maximum Likelihood Training of Score-Based Diffusion Models. In Advances in Neural Information Processing Systems, volume 34, pages 1415–1428. Curran Associates, Inc., 2021. URL https://proceedings.neurips.c/paper/2021/hash/0a9fdbb17feb6ccb7ec405cfb85222c4-Abstract.html. [56] D. J. Sutherland, H.-Y. Tung, H. Strathmann, S. De, A. Ramdas, A. Smola, and A. Gretton. Generative Models and Model Criticism via Optimized Maximum Mean Discrepancy. Mar. 2021. URL https://openreview.net/forum?id=HJWHIKqgl. [57] S. Verdú. Mismatched Estimation and Relative Entropy. IEEE Transactions on Information Theory, 56(8):3712–3720, Aug. 2010. ISSN 1557-9654. doi: 10.1109/TIT.2010.2050800. URL https://ieeexplore.ieee.org/document/5508632. [58] Y. Wang, L. Wang, S. Liu, and T. Suzuki. Zero-Flow Two-Sample Tests, July 2026. URL http://arxiv.org/abs/2607.21542. arXiv:2607.21542 [cs.LG] version: 1. [59] J. L. Watson, D. Juergens, N. R. Bennett, B. L. Trippe, J. Yim, H. E. Eisenach, W. Ahern, A. J. Borst, R. J. Ragotte, L. F. Milles, B. I. M. Wicky, N. Hanikel, S. J. Pellock, A. Courbet, W. Sheffler, J. Wang, P. Venkatesh, I. Sappington, S. V. Torres, A. Lauko, V. De Bortoli, E. Mathieu, S. Ovchinnikov, R. Barzilay, T. S. Jaakkola, F. DiMaio, M. Baek, and D. Baker. De novo design of protein structure and function with RFdiffusion. Nature, 620(7976):1089–1100, Aug. 2023. ISSN 1476-4687. doi: 10.1038/s41586-023-06415-8. [60] H. Xiao, K. Rasul, and R. Vollgraf. Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017. [61] J. Yang, R. Shi, D. Wei, Z. Liu, L. Zhao, B. Ke, H. Pfister, and B. Ni. Medmnist v2-a large-scale lightweight benchmark for 2d and 3d biomedical image classification. Scientific Data, 10(1):41, 2023. [62] H. Zhang, J. Zhou, Y. Lu, M. Guo, P. Wang, L. Shen, and Q. Qu. The Emergence of Reproducibility and Consistency in Diffusion Models. In Proceedings of the 41st International Conference on Machine Learning, pages 60558–60590. PMLR, July 2024. URL https://proceedings.mlr.press/v235/zhang24cn.html.