arXiv CSopen access

PosteriorBench: From Point Estimates to Posterior Matching in Evaluating Generative Inverse Solvers

PosteriorBench: From Point Estimates to Posterior Matching in Evaluating Generative Inverse Solvers

arXiv:2609.20794v1 [cs.LG] 17 Sep 2026

Jiachen Yao1∗ Zi-Siang Hsu2∗ Xi Deng1∗ Aditi Gupta3 Xin Ju4,5 Sally M. Benson4,5 Gege Wen6,5 Anima Anandkumar1 1 California Institute of Technology 2 National Taiwan University 3 Lawrence Berkeley National Laboratory 4 Department of Energy Science and Engineering, Stanford University 5 EarthFlow AI, Inc. 6 Department of Earth Sciences and Engineering, Imperial College London

Abstract Generative models are increasingly used to solve scientific inverse problems, but existing evaluations still focus primarily on whether a method can produce a single plausible reconstruction. This is insufficient for ill-posed problems, where multiple solutions may be consistent with the same sparse or noisy observations. In these settings, a method can achieve strong pointwise accuracy while still failing to capture the true posterior through mode collapse, overconfident uncertainty, or averaging incompatible solutions. We introduce PosteriorBench, a benchmark for evaluating the distributional accuracy of generative inverse solvers. PosteriorBench evaluates four physics-based inverse problems: Darcy flow inversion, Poisson source recovery, carbon capture and storage, and light transport material inference. For each task, we construct high-fidelity reference posteriors using computationally heavy but established procedures such as rejection sampling and Markov chain Monte Carlo, enabling direct assessment of whether solvers recover the full set of solutions rather than the single best sample. We pair these references with a five-metric posterior evaluation suite: posterior-mean error, posterior-standarddeviation error, maximum mean discrepancy, sliced Wasserstein distance, and radially averaged power-spectrum error. Together, these metrics assess pointwise accuracy, marginal uncertainty, distributional alignment, and global frequency fidelity. The benchmark spans sparse sensing, low-resolution observations, nonlinear forward models, varying noise levels, and multimodal priors, with a unified pipeline for distribution matching and uncertainty quantification. Our experiments reveal substantial distribution-matching gaps across current solvers, while showing that neural operators improve resolution robustness and that guidance weights and generation noise are key to posterior-variance calibration. The code is available at https://github.com/neuraloperator/PosteriorBench.

1

Introduction

Inverse problems are central to scientific computing: one observes indirect, partial, or noisy measurements and seeks to infer the parameter field that produced them. They arise in subsurface flow, optical imaging, fluid dynamics, and many other domains where direct measurement is expensive or impossible [1, 2]. The difficulty is not only that the forward physics may be nonlinear and expensive, but also that the inverse map is usually non-unique. Bayesian inverse problems make this ambiguity explicit by placing a prior p(x) over the unknown field x and conditioning on the observed data yobs through Bayes’ rule, p(x | yobs ) ∝ p(yobs | x)p(x). (1) ∗ Equal contribution.

Preprint.

PosteriorBench Posterior matching for generative inverse solvers Tasks Darcy flow Poisson source Carbon storage Light transport

Applications

Reference Prior ensemble Forward simulation Likelihood weighting

Binary permeability

Solvers

Evaluation

Diffusion models

Metrics Mean / std MMD / SWD Spectral error Obs. residual Runtime

Decoupled sampling Ensemble assimilation Stochastic surrogates

Smooth source fields

Sparse well observations

Ablations Noise level Guidance Resolution Time budget

Optical material properties

Figure 1: Overview of PosteriorBench. Across four scientific inverse tasks, each case pairs a fixed observation with a weighted reference posterior and a solver-generated ensemble. Evaluation compares posterior moments, distributional alignment, spatial structure, observation consistency, and computational cost. Here the likelihood p(yobs | x) is induced by the forward measurement model, while the prior p(x) encodes data distribution assumptions [3, 4]. Recent generative models, especially diffusion and flow-based models, have made this Bayesian decomposition increasingly tractable. Score-based generative models provide expressive priors for high-dimensional data [5–7], and diffusion-based sampling methods combine such priors with measurement likelihoods at inference time [8–12]. In scientific settings, this area now intersects with neural operators and physics-informed learning [13–16], leading to solvers that use joint coefficientsolution diffusion models, PDE residual guidance, or function-space formulations for physical fields [17–22]. These methods can produce ensembles, not merely point estimates, and are now being proposed as posterior samplers for PDE-constrained inverse problems. Existing scientific inverse-problem benchmarks such as InverseBench [23] evaluate plug-and-play diffusion samplers across physical inverse problems, but still center on a single solution for each case. More broadly, evaluation in this area is often anchored to single held-out ground-truth, which are insufficient for ill-posed problems. A solver can match the observations and still underestimate posterior variance or blur out incompatible modes into an unphysical mean. This gap motivates PosteriorBench, a benchmark for evaluating whether generative inverse solvers recover the posterior they claim to sample from. PosteriorBench focuses specifically on posterior matching: each benchmark case is paired with a transparent, computationally expensive reference posterior, so solvers can be evaluated as posterior samplers rather than only by their best reconstruction. Generated ensembles are compared to the reference posteriors using metrics for distributional moments, spatial structure, observation consistency, and computational cost. This framing is especially important for scientific settings, where posterior uncertainty informs decision making and downstream physical interpretation. The benchmark contains four tasks from different scientific domains: Darcy flow inversion, Poisson source recovery, carbon capture and storage (CCS), and light transport material inference (LTMI). Together they cover binary and smooth priors, sparse point and column observations, nonlinear physical maps, low-resolution measurements, and structured scientific priors. The CCS task is a high-impact subsurface monitoring setting, where permeability fields must be inferred from sparse well observations of CO2 dynamics; reducing the number of wells can substantially lower field intervention and monitoring cost. Ensemble methods such as ES-MDA remain important domain baselines for this inverse problem [24, 25]. The LTMI task adds a two-layer radiative inverse problem motivated by atmospheric imaging, optical tomography, and nondestructive material inspection, where optical measurements must be translated into plausible internal material structure. From the benchmark, we find that function-space diffusion samplers such as DDIS [22] and FunDPS [20] are strong posterior samplers across several scientific inverse tasks. We also find that posterior-mean error alone can be misleading: a solver may place the ensemble center near the reference mean while still misrepresenting posterior spread. This makes moment consistency a joint 2

requirement, where posterior mean and standard deviation errors should be interpreted together and, when possible, alongside distributional metrics such as MMD and SWD. We further compare posterior metrics against a traditional pointwise reconstruction metric. This comparison reveals that aggressively fitting a single reference field can degrade posterior structure: low pointwise error can coincide with poor posterior variance, distributional mismatch, or distorted spatial statistics. Through diagnostics based on the governing PDE constraint and the latent GRF smoothness parameter, we find that posterior metrics better reflect physically meaningful recovery than pointwise error alone. We also study how solvers respond to different observation noise levels. For a Gaussian observation model, guidance weights should, in theory, scale with the inverse observation-noise variance. Our sweeps show that learned samplers do not resolve posterior calibration by this scaling alone: stronger guidance improves observation consistency but often underestimates posterior variance, while weaker guidance preserves diversity at the cost of a biased or weakly conditioned posterior mean. This mean–variance tradeoff motivates conditioning mechanisms that can calibrate posterior mean and uncertainty jointly, rather than relying on a single weight. Section 4 details more ablation studies. Our contributions are threefold. First, we introduce a distribution-centered evaluation protocol for generative scientific inverse solvers that explicitly incorporates reference posterior construction and validation. Second, we organize four diverse inverse tasks under a common distributional evaluation suite, and we release code for reproducing these experiments and evaluating new solvers on PosteriorBench. Third, our experiments and ablations identify persistent distribution-matching gaps and insights into future method development in physics-based inverse problems.

2

Preliminaries

An inverse problem seeks an unknown physical field x ∈ X from indirect observations yobs ∈ Y produced by a forward map A : X → Y, such as a PDE solver, reservoir simulator, or radiative transport model. Because observations are often sparse or noisy and A is often many-to-one, the scientifically relevant target is usually not a single reconstruction but the posterior distribution (1). In PosteriorBench, a task or problem denotes a family of inverse problems, while a case denotes one observation-conditioned posterior recovery setting. For each case, a solver receives yobs and returns samples intended to approximate Equation (1). Diffusion models learn the prior by perturbing clean samples x0 ∼ p(x) into noisy variables xt and training a network sθ (xt , t) to approximate the score ∇xt log pt (xt ) [5–7]. For a forward noising SDE dxt = f (xt , t)dt + g(t)dwt , reverse sampling follows   dxt = f (xt , t) − g(t)2 sθ (xt , t) dt + g(t)dw̄t , t : T → 0, (2) where w̄t denotes Brownian motion in reverse time. To condition such a prior on observations, diffusion posterior sampling [8] uses the posterior-score decomposition ∇x log p(x | yobs ) = ∇x log p(x) + ∇x log p(yobs | x),

(3)

combining the learned diffusion score with an inference-time likelihood or guidance term [8–12]. For the additive Gaussian observation model η ∼ N (0, σy2 I),

(4)

1 2 ∇x ∥A(x) − yobs ∥2 . 2σy2

(5)

yobs = A(x) + η, the likelihood term is ∇x log p(yobs | x) = −

Practical diffusion posterior samplers usually apply this gradient to a denoised estimate x̂0 (xt ) during the reverse diffusion trajectory, so that each reverse step balances prior against agreement with the observed data. In scientific inverse problems, A may be differentiable, replaced by a neuraloperator surrogate or a physics residual, leading to variants such as joint coefficient-solution diffusion, function-space guidance, and decoupled prior sampling [17, 20–22]. 3

Darcy Obs (u)

Poisson Obs (φ )

CCS Obs (Dynamics)

LT Obs (u1)

LT Unknown (u2)

Posterior (a)

Posterior ( f )

Posterior (Geo-model)

Posterior (σt1)

Posterior (σt2)

0.78%

0.78%

1.00%

16×16

16×16

Figure 2: Representative cases from PosteriorBench. The top row shows the observation available to the inverse solver: sparse pressure measurements for Darcy flow and Poisson source recovery, sparse well-column measurements for CCS, and a low-resolution filter for light-transport material inference (LTMI). The bottom row shows multiple samples from the corresponding reference posterior.

3

PosteriorBench

PosteriorBench is designed to evaluate whether a generative inverse solver recovers the posterior distribution rather than only a single accurate reconstruction, as illustrated by the representative cases in Figure 2. The benchmark emphasizes settings where ambiguity is intrinsic, posterior mass is multimode, and uncertainty quantification is scientifically meaningful. At a high level, the benchmark is built around three principles: controllable inverse problems with known physics, high-fidelity reference posteriors, and evaluation metrics that compare distributions rather than point estimates. 3.1

Benchmark Task Design

PosteriorBench is organized around inverse problems in which sparse or aggregated measurements admit multiple physically plausible latent fields. We choose four tasks that vary along the main axes that affect posterior recovery: the prior over unknown fields, the observation pattern, the forward physics, and the dominant source of posterior ambiguity. Table 1 summarizes these design choices, while the following subsections describe the construction of each task. Table 1: Keyword feature comparison of the PosteriorBench tasks. Each row summarizes the prior, observation pattern, forward model, and dominant ambiguity source used in the benchmark. Task

Prior type

Observation

Forward map

Ambiguity

Darcy flow inversion Poisson source recovery CO2 capture and storage Light transp. material inference

Binary Smooth GRF Geostat. multimodal Cloud texture

Sparse sensors Sparse sensors Well columns Low-res optical

PDE solver PDE solver Reservoir sim. Radiative

Phase layout Spectral modes Plume uncertainty Material layout

3.1.1

Darcy Flow Inversion

Darcy flow is a canonical elliptic inverse problem for porous-media modeling [20]. On Ω = (0, 1)2 , we consider the steady equation −∇ · (a(x)∇u(x)) = 1,

x ∈ Ω,

u|∂Ω = 0,

(6)

where a(x) is the unknown permeability or conductivity field and u(x) is the corresponding pressure response. This setup uses constant unit forcing and homogeneous Dirichlet boundary conditions, 4

matching the standard Darcy data-generation convention. The inverse task is to recover the posterior distribution of a from sparse point observations of u. Following the standard Darcy construction used in operator-learning benchmarks [14], we sample a Gaussian random field (GRF) and threshold it into binary high- and low-conductivity phases, with representative values 12, 3. This prior creates discontinuous material interfaces and channel-like structures. Sparse observations constrain the induced flow field, but they do not uniquely determine the underlying phase layout: multiple connected high-conductivity pathways can produce similar pressure measurements at the observed locations. The Darcy task therefore evaluates whether a solver captures posterior uncertainty rather than a single plausible reconstruction. 3.1.2

Poisson Source Recovery

To assess solvers in a continuous and more spectrally variable setting, we design the Poisson source recovery task, governed by the equation −∇2 ϕ = f, (7) where f is the source term and ϕ is the potential. The goal is to recover the source term f from sparse observations of the potential ϕ. The workflow for posterior construction follows the Darcy task using sparse observation, but the prior introduces broader variation across cases. The source terms f are drawn from Gaussian random fields with varying correlation length τ and smoothness α. Sampling these hyperparameters independently for each case increases the structural diversity and spectral variability of the fields. The Poisson task therefore evaluates whether inverse solvers generalize across a broad family of prior spectra rather than overfitting to a single fixed prior. 3.1.3

Carbon Capture and Storage

Carbon capture and storage (CCS) requires reliable characterization of subsurface heterogeneity in order to forecast and monitor CO2 plume migration, pressure buildup, and storage security [26]. In PosteriorBench, the CCS task is formulated as a Bayesian inverse problem over a static permeability field. Let m denote the subsurface geomodel and let s = F (m) denote the dynamic CO2 saturation response after injection. Given sparse observations collected at monitoring wells, the goal is to recover the posterior distribution p(m | yobs ) rather than a single calibrated permeability map.

The task follows the sparse-monitoring regime used in recent function-space diffusion work for CCS [27]. Permeability realizations are generated from geostatistical priors using SGeMS-style simulation [28], and the corresponding CO2 saturation fields are produced with a high-fidelity reservoir simulator such as ECLIPSE [29]. Observations are represented as vertical strip or column patterns that mimic well measurements, so the inverse solver observes only a small fraction of the spatial domain. This setting is deliberately challenging: many geomodels can match the same well measurements, and the posterior can contain substantial spatial uncertainty away from the wells. For this task, the reference posterior is constructed by drawing a large candidate pool from the geological prior and accepting or weighting candidates according to their mismatch to the observed saturation response under the forward model or a validated neural-operator surrogate. The benchmark reports posterior matching against this reference ensemble, while also tracking observation consistency and reference-posterior validation diagnostics. Further details on the dataset, simulation, and posterior construction are provided in Section B.2. We note that ES-MDA is a widely used data-assimilation method in reservoir characterization, thus we include it in our baselines [24, 25]. 3.1.4

Light transport material inference (LTMI)

Inverse light transport through multilayer scattering materials arises in a range of applications, including atmospheric retrieval of cloud and aerosol structure, diffuse optical tomography of layered biological tissue, and nondestructive optical inspection of semitransparent materials. In these settings, the goal is to reconstruct spatially varying optical properties, such as extinction or scattering coefficients, from sparse observations. This inverse problem is challenging because the unknown field is high-dimensional, while the available measurements are often severely limited in viewpoint or acquisition time due to hardware and practical constraints. For example, in atmospheric imaging, it is often difficult to obtain simultaneous multi-view observations of the same cloud field. 5

In our dataset, each sample consists of a two-layer participating medium arranged along the viewing direction, where each layer contains a spatially varying density field with either cloud-like or cellularlike morphology. We then consider two kinds of measurement, reflectance fields and transmittance fields, each with low-resolution or pointwise observation. The goal is to infer the extinction coefficients of the participating material, which is ambiguous because different arrangements of scattering material can produce similar aggregate measurements under limited observation. 3.2

Reference Posterior Construction

A central feature of PosteriorBench is that each benchmark case is paired with a high-fidelity reference posterior. Depending on the problem structure, this reference distribution is obtained through slow but reliable procedures such as rejection sampling and Markov chain Monte Carlo. For rejectionsampling-based tasks, we draw a large candidate pool from the prior and simulate the observation process for each candidate. We specify the assumed Gaussian observation noise level σ and use it to assign soft likelihood weights   ∥A(xi ) − yobs ∥22 wi ∝ exp − . (8) 2σ 2 For efficient metric computation, we set threshold ϵ = 3σ and randomly keep 100 samples whose observation P mismatch is below the ϵ-threshold. We normalize the weights over the retained ensemble so that i wi = 1. Before using these reference posteriors as evaluation targets, we validate that they are stable and observation-consistent. The validation protocol is provided in Section C. 3.3

Evaluation Metrics

To evaluate inverse solvers as posterior samplers, we use a metric suite that compares generated ensembles with weighted reference posteriors. The reference distribution is represented as a weighted empirical distribution Pref = {(xi , wi )}M i=1 , where weights w are derived from rejection sampling or importance sampling. In contrast, the solvers under evaluation typically produce an unweighted ensemble of samples Pgen = {x′j }N j=1 . Considering the asymmetry, our metrics are as follows. 3.3.1

Marginal Moment Consistency

We first check whether the generated ensemble matches the low-order posterior statistics often used in downstream scientific analysis. The posterior mean measures the expected reconstructed field, while the marginal standard deviation measures the magnitude of pointwise uncertainty. For generated M samples Pgen = {x′j }N j=1 and a weighted reference posterior Pref = {(xi , wi )}i=1 , we define N

µ̄gen =

1 X ′ x , N j=1 j

µ̄ref =

M X

w i xi ,

(9)

i=1

and compute the relative mean error and, similarly, the std error: ∥µ̄gen − µ̄ref ∥2 ∥σ̄gen − σ̄ref ∥2 Errµ = , Errσ = . (10) ∥µ̄ref ∥2 ∥σ̄ref ∥2 We emphasize that the relative mean error differs from the mean squared error computed against a single target. These moment errors provide a straightforward check on whether a solver recovers not only the expected field, but also where the inverse problem remains uncertain. 3.3.2

Distributional and Geometric Alignment

To evaluate alignment beyond marginal moments, we use two complementary distributional measures. The first compares samples directly in the physical field space through a characteristic kernel, while the second compares one-dimensional projections after field-scale normalization. Maximum Mean Discrepancy (MMD) We use MMD to detect discrepancies in higher-order spatial statistics. Let the generated ensemble have normalized weights wj′ = 1/N and let the reference posterior have normalized weights wi . The squared MMD is X X X MMD2 = wi wi′ k(xi , xi′ ) + wj′ wj′ ′ k(x′j , x′j ′ ) − 2 wi wj′ k(xi , x′j ). (11) i,i′

j,j ′

i,j

6

The reported MMD is the square root of this quantity. The kernel is a multi-scale RBF kernel   1 X ∥x − y∥22 k(x, y) = exp − , S = {0.2, 0.5, 1.0, 2.0, 5.0}, (12) |S| ds s∈S

where d is the median squared distance between generated and reference samples. This multi-scale kernel makes the statistic sensitive to both coarse structural shifts and finer spatial discrepancies. Sliced Wasserstein Distance (SWD) To evaluate geometric proximity, we compute a sliced Wasserstein distance between generated and reference ensembles. Because the tasks have different physical units and dynamic ranges, SWD is computed after field-wise z-score normalization. We then project the normalized fields onto smooth random directions rather than i.i.d. pixel-wise Gaussian vectors to avoid local cancellation. This GRF projection distribution favors spatially coherent test functions, making SWD aligned with physical field discrepancies. For each direction vℓ , we compute the weighted one-dimensional Wasserstein distance between projected samples:   L N M X X 1X SWD = W1  wj′ δ⟨x̃′j ,vℓ ⟩ , wi δ⟨x̃i ,vℓ ⟩  . (13) L j=1 i=1 ℓ=1

We use L = 128 projections with default GRF parameters α = 2 and τ = 3. 3.3.3

Spectral Analysis

In physical inverse problems, the statistical texture and energy distribution across scales are of informative as well. We therefore compute the Radially Averaged Power Spectrum (RAPS) to evaluate spectral consistency. For each sample x ∈ RH×W , we first compute its 2D power spectrum Px (q) := |F (x)(q)|2 using the Discrete Fourier Transform (DFT). The 2D spectrum is qthen mapped to a 1D representation S(k) by averaging the power density within radial bins k = kx2 + ky2 . For the reference distribution Pref , the ensemble spectrum is defined as the weighted average: Sref (k) =

M X

wi Rxi (k),

Rxi (k) =

i=1

(14)

q∈Bk

For generated samples, we define Sgen (k) = N1 per-bin relative errors Errspec = exp

1 X Pxi (q), |Bk |

PN

1 X log |K| k∈K

j=1 Rxj (k) and report the geometric mean of



′

|Sgen (k) − Sref (k)| |Sref (k)|

! .

(15)

Here K denotes the valid frequency bins. The geometric mean is used to avoid overemphasizing high-frequency bins with low power. This metric assesses whether the solver preserves the physical energy cascade and avoids common pitfalls such as over-smoothing, which is often hard to tell from spatial-domain metrics.

4

Experiments

4.1

Methods

We evaluate eight probabilistic inverse solvers on all four PosteriorBench tasks: ECI-sampling [30], DiffusionPDE [17], FunDPS [20], Fun-DDPS [27], DDIS [22], FunDiff [21], ES-MDA [24, 25], and FNO with MC Dropout. Fun-DDPS and DDIS instantiate decoupled posterior sampling strategies that pair a learned prior with surrogate- or physics-based likelihood guidance, while ES-MDA provides a classical ensemble data-assimilation baseline and FNO with MC Dropout provides a direct neural uncertainty baseline. 7

Table 2: Main benchmark results across all PosteriorBench tasks. Lower is better for all metrics. Best and second-best results within each task and posterior-quality metric are highlighted in dark and light blue, respectively; runtime is rounded up to 0.1 min. Task

Method

Darcy

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff MC-dropout

Mean Error ↓

Std Error ↓

MMD ↓ 0.6457 0.5812 0.2040 0.2131 0.2268 0.1997 0.3229 0.6728

55.2020 48.4836 3.0778 2.7840 6.7372 2.9236 4.0833 6.0297

SWD ↓

Spectral Error ↓ 0.2007 0.3481 0.2016 0.1222 0.3252 0.1676 0.1225 0.1123

3.6 min 0.1 min 8.8 min 7.1 min 40.3 min 14.2 min 0.5 min 0.1 min

Poisson

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff MC-dropout

3.4434 1.0933 0.4720 0.1839 0.3319 0.1413 1.0529 0.3495

8.7178 5.6689 0.5862 0.4329 0.7004 0.3073 4.3679 0.8837

0.7623 0.6370 0.6692 0.3576 0.6814 0.2570 0.6363 0.7157

69.5794 41.5815 4.3900 2.4304 4.2661 1.8568 34.8920 4.4383

35.0923 7.7886 0.2303 0.0752 0.1907 0.1760 4.0411 0.2053

0.8 min 0.1 min 9.4 min 7.3 min 40.0 min 15.3 min 0.5 min 0.1 min

CCS

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff MC-dropout

0.6753 0.2457 0.1419 0.4659 0.2015 0.3277 0.1792 0.2574

1.5330 1.0903 0.2952 1.0305 0.3549 0.4332 0.5195 0.8926

0.4991 0.2941 0.2018 0.2410 0.2781 0.3007 0.3647 0.6844

10.5775 8.0556 4.2564 6.9947 5.8572 7.3551 7.1329 13.5715

88.8274 14.5449 0.6370 14.7840 0.7294 2.2043 0.6896 0.6555

0.3 min 0.1 min 6.3 min 6.2 min 34.0 min 9.2 min 1.5 min 0.1 min

LTMI

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff MC-dropout

0.3045 0.1250 0.1611 0.1699 0.1593 0.1504 0.1492 0.1554

0.3911 0.2007 0.2867 0.2374 0.2265 0.2265 0.3994 0.8542

0.4509 0.2842 0.3327 0.3326 0.3183 0.3033 0.3753 0.6707

8.8375 2.4285 2.8501 2.9389 2.9462 2.9999 3.1404 4.6755

0.1810 0.0612 0.0640 0.1828 0.1554 0.1869 0.2624 0.5060

0.1 min 0.1 min 5.4 min 4.8 min 10.8 min 4.8 min 0.4 min 0.1 min

4.2

0.5112 0.4519 0.0804 0.0767 0.1009 0.0620 0.1203 0.1990

1.4975 1.7547 0.4354 0.4178 0.7871 0.4130 0.5900 0.9509

Time ↓

Main Results

Table 2 summarizes the quantitative evaluation of posterior-generating solvers across the benchmark tasks. Section D provides full guidance-weight sweeps and case-level diagnostics, pairwise metric analyses, and standard deviations across cases. Mean-only summaries can be misleading. A central purpose of PosteriorBench is to prevent posterior evaluation from collapsing back to a single point-summary comparison. The LTMI results provide a concrete example: FNO with MC Dropout attains lower posterior-mean error than FunDPS in Table 2, yet its posterior-std error, MMD, and SWD remain high. Figure 3 shows why this is not a metric inconsistency: For a representative LTMI case, the MC Dropout sample is visibly over-smoothed relative to the reference posterior sample, while FunDPS better preserves the fine spatial texture of the material field. Thus, a low mean error can reflect agreement with a central tendency while still missing the geometry and spread of the posterior distribution. Marginal moment consistency should therefore be read jointly, together with metrics like MMD and SWD. Useful inductive bias in Function-space diffusion samplers. Across the benchmark, FunDPS, Fun-DDPS, and DDIS form a family of guided diffusion posterior samplers built on function-space score priors. At the same time, their performance is task-dependent: classical methods such as ES-MDA can be stronger in settings such as LTMI, and the best-performing function-space diffusion sampler varies across datasets. Their results suggest that function-space score priors are useful for posterior matching under sparse or low-resolution observations. The comparison with DiffusionPDE points to the value of function-space score backbones, whose spectral operator blocks are better 8

aligned with continuous physical fields than the grid-based backbone. The comparison with ECI highlights the role of soft likelihood guidance: ECI enforces observed entries directly through hard replacement, whereas guided diffusion samplers condition the sampling process on the observations while trying to stay on the prior manifold. During inference, these samplers dynamically balance learned prior against data consistency, which is central to recovering ensembles rather than only point reconstructions. The remaining variation across tasks and metrics motivates a closer comparison within this guided diffusion family. Reference

FunDPS

MC Dropout 5 4 3 2 1 0

Figure 3: LTMI case visualization for the first material field σt1 . The panels show one reference posterior sample, one FunDPS-generated sample, and one FNO with MC-Dropout-generated sample. Although MC Dropout attains lower LTMI posterior-mean error than FunDPS in Table 2, its sample is visibly over-smoothed relative to the reference structure, consistent with its high posterior-standarddeviation error, MMD, and SWD.

Darcy

Metric value

Mean error

Poisson

CCS

Std error

Light transport

MMD

SWD

Spectral error

10.0

1.5 0.4

0.6

7.5

0.4

0.4

5.0

0.2

1.0 0.2

0.5 0.2 0.25

0.50

0.75

0.25

0.50

0.75

2.5 0.25

0.50

0.75

0.0 0.25

0.50

0.75

0.25

0.50

0.75

Pointwise metric

Figure 4: Relationship between a traditional pointwise metric and the five posterior metrics. Each panel compares the pointwise metric with one posterior metric across methods and tasks. The V-shaped trends show that the distributional metrics are not monotone with respect to pointwise error. 4.3

Pointwise-Distributional Tradeoff

In addition to the five posterior metrics used in Table 2, we compare a traditional pointwise relative L2 metric against one reference sample. Figure 4 shows that the relationship is often V-shaped rather than monotone. At the low-pointwise-error side, pushing samples closer to a single reference can remove distributional structure, so distribution-sensitive errors increase even as the pointwise metric improves. At the high-error side, both pointwise and posterior errors are large, and the two families of metrics become positively associated. This pattern indicates that a pointwise metric alone can misidentify over-fitted or over-concentrated samples. The five posterior metrics are therefore useful because they expose the trade-off between single-reference reconstruction and distributional fidelity. 4.4

Qualitative Posterior Comparisons

Figure 5 compares posterior mean and standard deviation estimates on Poisson source recovery and Darcy flow inversion. Both FunDPS and DDIS recover the main posterior structures, while their error maps reveal larger discrepancies in high-gradient and high-uncertainty regions. Both models also exhibit conservative predictions. For Poisson source recovery, the mean error maps show residuals that pull extreme values toward zero, indicating that both models underestimate the 9

magnitude of the source extrema. This conservative behavior extends to uncertainty quantification: the standard-deviation error maps for both Poisson source recovery and Darcy flow inversion are predominantly negative, suggesting that the models systematically underestimate posterior variance. Beyond these shared traits, errors in Darcy flow inversion concentrate near sharp interface-like structures, highlighting a more challenging posterior landscape. Overall, DDIS produces more spatially balanced residuals and fewer large localized artifacts, though accurately capturing the full scale of posterior extremes and standard deviation remains a shared challenge. Poisson Source Recovery

Darcy Flow Inversion DDIS error

Reference

0.07

0.00

FunDPS error

DDIS error

Mean

FunDPS error

Mean

Reference

0.00

-5.3 4.0

Std. dev.

Std. dev.

-0.07 0.07

0.00

5.3

0.00

-0.07

-4.0

Figure 5: Qualitative posterior comparisons. Left: Poisson source recovery. Right: Darcy flow inversion. For each problem, the top row shows posterior mean µ and the bottom row shows posterior standard deviation σ. Within each group, we show the reference posterior statistic, FunDPS error, and DDIS error from left to right.

4.5

Out-of-Distribution Experiment

We evaluate out-of-distribution generalization on Poisson source recovery by varying the GRF parameter range used to train the Fun-DDPS prior. The prior is trained under three α regimes: single (α ≡ 2.25), narrow (α ∈ 2.25 ± 0.375), and full (α ∈ 2.25 ± 0.75). The differentiable surrogate used for likelihood guidance is trained on the full parameter range in all three settings, so the ablation isolates the effect of prior-training coverage. Table 3 reveals a counterintuitive split between pointwise and distributional behavior. The full prior gives the lowest pointwise relative L2 , but the distribution-sensitive metrics (std error, MMD, and SWD) worsen as the prior-training range expands. Thus, broader prior coverage might not automatically improve posterior matching when the prior must represent source fields from a larger, harder-to-learn range. To diagnose this behavior, Figure 6 evaluates physics and parameter consistency. (i) We compute the Poisson residual |A(ϕ̂) − fˆ| for each predicted pair (fˆ, ϕ̂). The full setting produces more high-residual outliers, indicating more frequent PDE violations. (ii) We estimate the GRF smoothness α from each generated source using a DCT-domain spectral likelihood. The estimates are consistently biased below the reference posterior; broader training ranges widen their distribution mainly toward lower values, revealing poor recovery of latent smoothness despite plausible pixel-space samples. Taken together, the PDE-residual and α-recovery diagnostics further show that lower pointwise error need not imply better posterior distribution matching. The distributional metrics align more closely with latent-parameter recovery and physical-consistency diagnostics than pointwise relative L2 alone. A second finding is that Poisson source recovery remains a demanding benchmark despite being generated from a synthetic GRF family. Variation in the latent smoothness parameter induces a sufficiently complex distribution that strong samplers do not fully recover the parameter structure. Table 3: Fun-DDPS Poisson out-of-distribution evaluation across GRF prior-training ranges. Single uses α = 2.25, Narrow uses α ∈ [1.875, 2.625], and Full uses α ∈ [1.5, 3.0]. Prior coverage Single Narrow Full

Pointwise Rel L2

Mean Rel L2

Std Rel L2

MMD

SWD

Spectral Rel L2

0.4176 0.3418 0.2956

0.1551 0.1487 0.1464

0.3235 0.3151 0.3720

0.2669 0.2879 0.3730

1.7041 1.8243 2.4415

0.0330 0.0358 0.0329

10

Single

0.0

0.1

0.2

0.3

0.4

Narrow

0.5

0.6

0.0

0.1

0.2

0.3

0.4

Full

0.5

0.6

0.0

0.1

0.2

1.4

1.6

1.8

0.3

0.4

0.5

0.6

2.4

2.6

2.8

PDE residual MAE Prediction

GT

1

case

10 20 30 40 50 1.4

1.6

1.8

2.0

2.2

2.4

2.6

2.8

3.0

1.4

1.6

1.8

2.0

2.2

2.4

2.6

2.8

3.0

2.0

2.2

3.0

alpha

Figure 6: Fun-DDPS Poisson OOD diagnostics. Top row: PDE residual diagnostic, where each point corresponds to a generated sample evaluated by using the predicted pair (fˆ, ϕ̂) and computing the residual A(ϕ̂) − fˆ under the Poisson operator. Bottom row: α-recovery diagnostic across single, narrow, and full prior-training ranges, where recovered α ranges are estimated from generated source fields using a DCT-domain GRF spectral likelihood and compared with the reference posterior ranges.

4.6

Observation Noise and Guidance Calibration

This experiment tests whether solver hyper-parameters that are often described as inverse noise scales actually calibrate to the posterior induced by different observation-noise levels. As shown in Equations (4) and (5), the Gaussian likelihood score scales with the likelihood precision σy−2 , which motivates scaling data-consistency guidance with the observation-noise level. However, the guidance coefficient used in practice may not follow this rule due to the interaction with discretization, normalization, annealing schedules, etc. For Poisson source recovery, we vary the observation-noise scale σ and report the guidance optimum σ,λ σ λ⋆m (σ) = arg minλ m(Pgen , Pref ). Figure 7 shows that the optimal guidance weight generally increases with the inverse observation-noise variance σ −2 for almost all metrics. Spectral error exhibits larger fluctuations, but still follows the same broad positive relation. This behavior is consistent with the Gaussian likelihood assumption. The full guidance-weight sweep and case-level spatial diagnostics are reported in Section D.1. The same sweep also shows that the optimum is not metric-invariant. At a fixed noise scale σ, the values of λ⋆m (σ) are different across metrics, especially between posterior-mean and posterior-std errors. Because the mean and standard deviation are complementary summaries of the same posterior distribution, this separation indicates that the guidance strength that best fits posterior center need not be variance-calibrated. Thus, a single scalar guidance weight faces trade off between mean accuracy and uncertainty calibration. 4.7

Resolution Ablation

The resolution ablation study in Table 4 shows that FunDPS consistently outperforms DiffusionPDE across all training settings when evaluated at 128 × 128. FunDPS achieves lower mean and standarddeviation errors, as well as improved MMD, SWD, and spectral metrics, with the best distributional alignment under multi-resolution training. In contrast, DiffusionPDE shows limited gains with increased resolution. These results indicate that function-space training improves robustness to train-test resolution changes. 11

Optimal guidance

400 Mean rel. L2 Std rel. L2 MMD SWD Spectral rel. L2

300 200

2 ref.

100 0

0.2

0.4

0.6

Inverse noise variance

0.8

2 (×108 )

1.0

Figure 7: Relationship between inverse observation-noise variance σ −2 and the optimal guidance weight for FunDPS on Poisson source recovery, with posterior threshold set to ϵ = 3σ. Each curve reports the guidance value that minimizes one evaluation metric at each inverse noise variance. The optimal guidance weights generally increase with σ −2 , indicating a positive association between likelihood precision and preferred guidance strength. The metric-specific optima also differ at the same σ −2 , notably between posterior-mean and posterior-standard-deviation errors. Table 4: Resolution generalization study on Darcy flow inversion. The table compares DiffusionPDE and FunDPS at each training resolution setting. The inference is conducted at 128 × 128 resolution. Resolution

Method

64 × 64

DiffusionPDE FunDPS

128 × 128

DiffusionPDE FunDPS

Multi-res

DiffusionPDE FunDPS

5

Mean error ↓

Std error ↓

MMD ↓

SWD ↓

Spectral error ↓

0.1838 ± 0.0433 0.1070 ± 0.0191

0.9324 ± 0.1902 0.4540 ± 0.0935

0.3361 ± 0.0506 0.2570 ± 0.0474

1.6860 ± 0.4621 0.8440 ± 0.2320

0.0690 ± 0.0478 0.0250 ± 0.0179

0.1833 ± 0.0423 0.0886 ± 0.0237

0.9324 ± 0.1860 0.4010 ± 0.0821

0.3368 ± 0.0494 0.2210 ± 0.0451

1.5465 ± 0.5450 0.7400 ± 0.1980

0.0661 ± 0.0485 0.0154 ± 0.0110

0.1745 ± 0.0319 0.0926 ± 0.0182

0.6708 ± 0.1016 0.3570 ± 0.0396

0.3306 ± 0.0425 0.2030 ± 0.0292

1.2667 ± 0.4895 0.8090 ± 0.2820

0.0506 ± 0.0348 0.0158 ± 0.0129

Discussion

PosteriorBench evaluates generative scientific inverse solvers as posterior samplers rather than as single-reconstruction methods. Across Darcy flow inversion, Poisson source recovery, carbon capture and storage, and light-transport material inference, the benchmark makes the residual ambiguity under partial observations explicit by comparing generated ensembles against reference posterior distributions. The experiments highlight two main findings. First, PosteriorBench exposes concrete regimes in which pointwise and distributional evaluation disagree. Second, function-space diffusion samplers provide strong posterior-matching performance across the benchmark, yet still leave gaps in jointly matching posterior means and standard deviations, recovering latent parameters, and satisfying PDE constraints. These findings highlight that a solver can match observations or point estimates while still misrepresenting posterior uncertainty, spatial structure, or distinct meaningful modes. The present benchmark is necessarily limited by the cost of constructing high-fidelity reference posteriors and the number of solvers evaluated so far. Our broader vision is for PosteriorBench to serve as an evolving evaluation platform expanding through community effort, while ultimately shifting the field away from single point reconstructions and toward calibrated, reproducible posterior recovery for scientific decision making under uncertainty.

Acknowledgments and Disclosure of Funding Anima Anandkumar is supported in part by Bren endowed chair, ONR (MURI grant N00014-231-2654), and the AI2050 senior fellow program at Schmidt Sciences. Jiachen Yao is supported in part by the Naren and Vinita Gupta Fellowship. The authors thank Modal for providing part of the compute credits. 12

References [1] Albert Tarantola. Inverse problem theory and methods for model parameter estimation. SIAM, 2005. [2] Jennifer L Mueller and Samuli Siltanen. Linear and nonlinear inverse problems with practical applications. SIAM, 2012. [3] Simon L Cotter, Massoumeh Dashti, James Cooper Robinson, and Andrew M Stuart. Bayesian inverse problems for functions and applications to fluid mechanics. Inverse problems, 25(11):115008, 2009. [4] Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin. Bayesian data analysis. Chapman and Hall/CRC, 1995. [5] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International conference on machine learning, pages 2256–2265. PMLR, 2015. [6] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840–6851, 2020. [7] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456, 2020. [8] Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. arXiv preprint arXiv:2209.14687, 2022. [9] Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz. Pseudoinverse-guided diffusion models for inverse problems. In International Conference on Learning Representations, 2023. [10] Zihui Wu, Yu Sun, Yifan Chen, Bingliang Zhang, Yisong Yue, and Katherine Bouman. Principled probabilistic imaging using diffusion models as plug-and-play priors. Advances in Neural Information Processing Systems, 37:118389–118427, 2024. [11] Bingliang Zhang, Wenda Chu, Julius Berner, Chenlin Meng, Anima Anandkumar, and Yang Song. Improving diffusion inverse problem solving with decoupled noise annealing. In Proceedings of the Computer Vision and Pattern Recognition Conference, pages 20895–20905, 2025. [12] Giannis Daras, Hyungjin Chung, Chieh-Hsin Lai, Yuki Mitsufuji, Jong Chul Ye, Peyman Milanfar, Alexandros G Dimakis, and Mauricio Delbracio. A survey on diffusion models for inverse problems. arXiv preprint arXiv:2410.00083, 2024. [13] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. [14] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020. [15] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1–97, 2023. [16] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1–27, 2024. [17] Jiahe Huang, Guandao Yang, Zichen Wang, and Jeong Joon Park. Diffusionpde: Generative pde-solving under partial observation. arXiv preprint arXiv:2406.17763, 2024. 13

[18] Dule Shu, Zijie Li, and Amir Barati Farimani. A physics-informed diffusion model for highfidelity flow field reconstruction. Journal of Computational Physics, 478:111972, 2023. [19] Christian Jacobsen, Yilin Zhuang, and Karthik Duraisamy. Cocogen: Physically consistent and conditioned score-based generative models for forward and inverse problems. SIAM Journal on Scientific Computing, 47(2):C399–C425, 2025. [20] Jiachen Yao, Abbas Mammadov, Julius Berner, Gavin Kerrigan, Jong Chul Ye, Kamyar Azizzadenesheli, and Anima Anandkumar. Guided diffusion sampling on function spaces with applications to pdes. In Advances in Neural Information Processing Systems, 2025. [21] Sifan Wang, Zehao Dou, Tong-Rui Liu, and Lu Lu. Fundiff: Diffusion models over function spaces for physics-informed generative modeling, 2025. [22] Thomas YL Lin, Jiachen Yao, Lufang Chiang, Julius Berner, and Anima Anandkumar. Decoupled diffusion sampling for inverse problems on function spaces. arXiv preprint arXiv:2601.23280, 2026. [23] Hongkai Zheng, Wenda Chu, Bingliang Zhang, Zihui Wu, Austin Wang, Berthy Feng, Caifeng Zou, Yu Sun, Nikola Borislavov Kovachki, Zachary E Ross, Katherine Bouman, and Yisong Yue. Inversebench: Benchmarking plug-and-play diffusion models for scientific inverse problems. In The Thirteenth International Conference on Learning Representations, 2025. [24] Alexandre A Emerick and Albert C Reynolds. Ensemble smoother with multiple data assimilation. Computers & Geosciences, 55:3–15, 2013. [25] Seungpil Jung, Kyungbook Lee, Changhyup Park, and Jonggeun Choe. Ensemble-based data assimilation in reservoir characterization: A review. Energies, 11(2):445, 2018. [26] Stephen Pacala and Robert Socolow. Stabilization wedges: solving the climate problem for the next 50 years with current technologies. Science, 305(5686):968–972, 2004. [27] Xin Ju, Jiachen Yao, Anima Anandkumar, Sally M Benson, and Gege Wen. Function-space decoupled diffusion for forward and inverse modeling in carbon capture and storage. In AI&PDE: ICLR 2026 Workshop on AI and Partial Differential Equations, 2026. [28] Nicolas Remy, Alexandre Boucher, and Jianbing Wu. Applied geostatistics with SGeMS: A user’s guide. Cambridge University Press, 2009. [29] Schlumberger. Eclipse reservoir simulation software: Reference manual, 2009. [30] Chaoran Cheng, Boran Han, Danielle C. Maddix, Abdul Fatir Ansari, Andrew Stuart, Michael W. Mahoney, and Bernie Wang. Gradient-free generation for hard-constrained systems. In The Thirteenth International Conference on Learning Representations, 2025. [31] Charles W Groetsch and CW Groetsch. Inverse problems in the mathematical sciences, volume 52. Springer, 1993. [32] James V Beck, Ben Blackwell, and Charles R St Clair. Inverse heat conduction: Ill-posed problems. James Beck, 1985. [33] Marco A Iglesias, Kody JH Law, and Andrew M Stuart. Ensemble kalman methods for inverse problems. Inverse Problems, 29(4):045001, 2013. [34] Gabriel Cardoso, Yazid Janati El Idrissi, Sylvain Le Corff, and Eric Moulines. Monte carlo guided diffusion for bayesian linear inverse problems. arXiv preprint arXiv:2308.07983, 2023. [35] Zehao Dou and Yang Song. Diffusion posterior sampling for linear inverse problem solving: A filtering perspective. In The Twelfth International Conference on Learning Representations, 2024. [36] Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193, 2019. 14

[37] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755–6766, 2020. [38] Qianying Cao, Somdatta Goswami, and George Em Karniadakis. Laplace neural operator for solving differential equations. Nature Machine Intelligence, 6(6):631–640, 2024. [39] Zijie Li, Kazem Meidani, and Amir Barati Farimani. Transformer for partial differential equations’ operator learning. arXiv preprint arXiv:2205.13671, 2022. [40] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388):1–26, 2023. [41] Roberto Molinaro, Yunan Yang, Björn Engquist, and Siddhartha Mishra. Neural inverse operators for solving pde inverse problems. arXiv preprint arXiv:2301.11167, 2023. [42] David J. C. MacKay. A practical bayesian framework for backpropagation networks. Neural Computation, 4(3):448–472, 1992. [43] Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In Proceedings of the 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059. PMLR, 2016. [44] Liu Yang, Xuhui Meng, and George Em Karniadakis. B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data. Journal of Computational Physics, 425:109913, 2021. [45] Guang Lin, Christian Moya, and Zecheng Zhang. B-DeepONet: An enhanced bayesian DeepONet for solving noisy parametric PDEs using accelerated replica exchange SGLD. Journal of Computational Physics, 473:111713, 2023. [46] Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet, and Mohammad Norouzi. Palette: Image-to-image diffusion models. In ACM SIGGRAPH 2022 conference proceedings, pages 1–10, 2022. [47] Yusuke Tashiro, Jiaming Song, Yang Song, and Stefano Ermon. Csdi: Conditional score-based diffusion models for probabilistic time series imputation. Advances in neural information processing systems, 34:24804–24816, 2021. [48] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. Advances in neural information processing systems, 35:23593–23606, 2022. [49] Yinhuai Wang, Jiwen Yu, and Jian Zhang. Zero-shot image restoration using denoising diffusion null-space model. arXiv preprint arXiv:2212.00490, 2022. [50] Jiaming Song, Qinsheng Zhang, Hongxu Yin, Morteza Mardani, Ming-Yu Liu, Jan Kautz, Yongxin Chen, and Arash Vahdat. Loss-guided diffusion models for plug-and-play controllable generation. In International Conference on Machine Learning, pages 32483–32498. PMLR, 2023. [51] Morteza Mardani, Jiaming Song, Jan Kautz, and Arash Vahdat. A variational perspective on solving inverse problems with diffusion models. arXiv preprint arXiv:2305.04391, 2023. [52] Jan-Hendrik Bastek, WaiChing Sun, and Dennis M. Kochmann. Physics-informed diffusion models, 2025. [53] Yi Zhang and Difan Zou. Physics-informed distillation of diffusion models for pde-constrained generation, 2025. [54] Zeyu Li, Hongkun Dou, Shen Fang, Wang Han, Yue Deng, and Lijun Yang. Physics-aligned field reconstruction with diffusion bridge. In The Thirteenth International Conference on Learning Representations, 2025. 15

[55] Peiyan Hu, Rui Wang, Xiang Zheng, Tao Zhang, Haodong Feng, Ruiqi Feng, Long Wei, Yue Wang, Zhi-Ming Ma, and Tailin Wu. Wavelet diffusion neural operator, 2025. [56] Marc Amorós-Trepat, Luis Medrano-Navarro, Qiang Liu, Luca Guastoni, and Nils Thuerey. Guiding diffusion models to reconstruct flow fields from sparse data. Physics of Fluids, 38(1), January 2026. [57] Aliaksandra Shysheya, Cristiana Diaconu, Federico Bergamin, Paris Perdikaris, José Miguel Hernández-Lobato, Richard Turner, and Emile Mathieu. On conditional diffusion models for pde simulations. Advances in Neural Information Processing Systems, 37:23246–23300, 2024. [58] Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking simulation-based inference. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 343–351. PMLR, 2021. [59] Antoine Wehenkel, Juan L. Gamella, Ozan Sener, Jens Behrmann, Guillermo Sapiro, JoernHenrik Jacobsen, and Marco Cuturi. Addressing misspecification in simulation-based inference through data-driven calibration. In Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research, pages 65949– 65980. PMLR, 2025. [60] Martin Zach, Youssef Haouchat, and Michael Unser. A statistical benchmark for diffusion posterior sampling algorithms. arXiv preprint arXiv:2509.12821, 2025. [61] Yinhao Zhu and Nicholas Zabaras. Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics, 366:415–447, 2018. [62] Guido Di Federico and Louis J Durlofsky. Latent diffusion models for parameterization of faciesbased geomodels and their use in data assimilation. Computers & Geosciences, 194:105755, 2025. [63] Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. U-fno—an enhanced fourier neural operator-based deep-learning model for multiphase flow. Advances in Water Resources, 163:104180, 2022. [64] Gege Wen, Zongyi Li, Qirui Long, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. Real-time high-resolution co2 geological storage prediction using nested fourier neural operators. Energy & Environmental Science, 16(4):1732–1741, 2023. [65] Gabriel S. Seabra, Nikolaj T. Mücke, Vinicius L. S. Silva, Denis Voskov, and Femke C. Vossepoel. Ai enhanced data assimilation and uncertainty quantification applied to geological carbon storage. International Journal of Greenhouse Gas Control, 136:104190, 2024. [66] Su Jiang and Louis J Durlofsky. History matching for geological carbon storage using data-space inversion with spatio-temporal data parameterization. International Journal of Greenhouse Gas Control, 134:104124, 2024. [67] Yifu Han, François P Hamon, Su Jiang, and Louis J Durlofsky. Surrogate model for geological co2 storage and its use in hierarchical mcmc history matching. Advances in Water Resources, 187:104678, 2024. [68] Wenchao Teng and Louis J Durlofsky. Likelihood-free inference and hierarchical data assimilation for geological carbon storage. Advances in Water Resources, 201:104961, 2025. [69] Zhongzheng Wang, Yuntian Chen, Wenhao Fu, Mengge Du, Guodong Chen, Xiaopeng Ma, and Dongxiao Zhang. Generative inverse modeling for improved geological co2 storage prediction via conditional diffusion models. Applied Energy, 395:126071, 2025. [70] Zhao Feng, Xin-Yang Liu, Meet Hemant Parikh, Junyi Guo, Pan Du, Bicheng Yan, and JianXun Wang. Generative latent diffusion model for inverse modeling and uncertainty analysis in geological carbon sequestration. arXiv preprint arXiv:2508.16640, 2025. 16

[71] Miguel Liu-Schiaffini, Julius Berner, Boris Bonev, Thorsten Kurth, Kamyar Azizzadenesheli, and Anima Anandkumar. Neural operators with localized integral and differential kernels. arXiv preprint arXiv:2402.16845, 2024.

17

A

Related Work

Bayesian scientific inverse problems. Inverse problems are classically formulated as the recovery of unknown parameters or fields from indirect observations through a forward physical model [1, 2, 31, 32]. The Bayesian formulation treats the unknown as a random function and combines a prior with the observation likelihood to obtain a posterior distribution [3, 4]. This view is essential when observations are sparse, noisy, or non-identifying, since posterior uncertainty can remain large even when the observations are matched. Classical data-assimilation and ensemble methods, including ensemble Kalman approaches, provide scalable approximations for some high-dimensional systems but can be limited by Gaussian or linearized update assumptions [33]. In high-dimensional scientific settings, asymptotically reliable samplers such as MCMC, sequential Monte Carlo, or rejection sampling are often too expensive to use as practical solvers, which motivates amortized or generative approximations [34, 35]. PosteriorBench uses these slow methods not as deployment algorithms, but as reference procedures for evaluating whether faster learned samplers recover the intended posterior. Neural operators and PDE inverse solvers. Physics-informed neural networks and neural operators have become standard tools for learning PDE solution maps and solving PDE-constrained inverse problems [13–15]. Operator-learning architectures such as FNO and DeepONet learn mappings between function spaces and can generalize across discretizations more naturally than fixed-grid networks [14, 15, 36]. Subsequent architectures extend this idea with multipole, Laplace, transformer, and geometry-aware operator designs [37–40]. Physics-informed neural operators further add PDE residuals or weak supervision to improve generalization when paired data are limited [16]. For inverse problems, deterministic neural solvers can be accurate when the target is a point estimate, but they do not by themselves represent the full posterior over plausible fields. For example, neural inverse operators directly learn maps from observations to unknown coefficients [41]. PosteriorBench treats these models both as components of generative solvers and as important baselines or surrogates. Uncertainty-aware neural samplers. Uncertainty-aware neural predictors provide another route to posterior sampling without training a full generative inverse model. Bayesian neural networks place distributions over network weights and infer a weight posterior, so predictive variation reflects model uncertainty [42]. MC dropout offers a scalable approximation by keeping dropout active at inference time and using repeated stochastic forward passes as approximate Bayesian predictions [43]. In scientific machine learning, Bayesian PDE solvers such as B-PINNs [44] and B-DeepONet [45] combine Bayesian neural networks or Bayesian operator learning with PDE constraints to quantify uncertainty in forward and inverse settings. These approaches are complementary to PosteriorBench, but their inverse examples typically recover single- or low-dimensional quantities, whereas PosteriorBench asks solvers to sample posteriors over 64×64 or 128×128 fields. This high dimensionality is one reason generative inverse samplers are attractive, and hence why the benchmark focuses on them. B-PINNs also require fitting a separate model for each observation case, which does not scale to the amortized multi-case evaluation used here. Diffusion posterior sampling. Diffusion and score-based models were first developed as powerful unconditional generators [5–7], then adapted to inverse problems by conditioning a learned prior on measurements. Conditional models learn task-specific conditional distributions [46, 47], whereas plug-and-play posterior samplers reuse an unconditional prior with an observation model at test time. For linear or image-domain inverse problems, methods such as DDRM, DDNM, pseudoinverse, loss guidance, and RED-diff provide different mechanisms for combining denoising priors with data consistency [9, 10, 48–51]. Sequential Monte Carlo and filtering perspectives seek stronger posterior correctness guarantees for some classes of inverse problems [34, 35]. DAPS reduces approximation error by decoupling denoising and likelihood updates [11]. These advances motivate distributional evaluation, but many reported results still emphasize reconstruction quality, perceptual quality, or observation consistency rather than calibrated posterior matching. Diffusion models for scientific applications. Scientific inverse problems add challenges that are muted in natural-image restoration: the forward map may be a PDE solver, observations may live in a different physical field than the unknown, and the state is more naturally a function than a fixedresolution image. DiffusionPDE models paired physical fields under partial observation [17]; physicsinformed diffusion methods add residual or constraint guidance [18, 52, 53]; and CoCoGen constructs conditioned score-based models for forward and inverse physical problems [19]. Function-space 18

approaches such as FunDPS and FunDiff address discretization dependence by defining generative modeling over functions rather than only arrays [20, 21]. DDIS and related decoupled designs separate prior learning from the physics-induced likelihood using a neural operator surrogate [22]. Other recent work explores diffusion bridges, wavelet diffusion operators, sparse flow-field reconstruction, and conditional PDE simulation [54–57]. PosteriorBench provides a common setting for comparing these design choices as posterior samplers. Benchmarks for scientific inverse problems. Scientific inverse-problem benchmarks increasingly extend beyond natural-image restoration. InverseBench tests plug-and-play diffusion priors with physical forward models, but is limited to single-reference reconstruction [23]. Simulation-based inference benchmarks more directly assess posterior estimation, but often consider lower-dimensional parameters and tractable likelihoods [58]; learned posterior targets can also make scores depend on model architecture, training coverage, and calibration [59]. Zach et al. enable precise distributional tests for one-dimensional Bayesian linear inverse problems with 64 discretization points and Lévyprocess priors, whose posteriors admit efficient Gibbs sampling [60]. PosteriorBench instead targets scientific function-space posteriors on 64×64 or 128×128 fields under sparse, low-resolution, or column observations generated by expensive physics. It constructs reference distributions through prior sampling, physical simulation, observation matching, and importance weighting, enabling posterior-level evaluation across multiple physics domains, including CCS and LTMI. Carbon capture and storage. Carbon capture and storage is a major climate-mitigation technology, but safe deployment requires uncertainty-aware monitoring of subsurface CO2 migration and pressure buildup [26]. Reservoir characterization and data assimilation have long relied on ensemble Kalman methods and ensemble smoothers, including ES-MDA, because they can update high-dimensional geomodel ensembles from sparse monitoring data [24, 25]. These methods are practical and domainrelevant but can struggle with strongly non-Gaussian geological priors and channelized or facies-like structures. Deep generative parameterizations and learned geomodel priors have been explored as a way to represent non-Gaussian reservoir structure within data-assimilation workflows [61, 62]. Neural-operator surrogates have accelerated geological CO2 storage simulation [63, 64], and recent work combines generative models with data assimilation, history matching, or inverse modeling for geological carbon storage [65–70]. PosteriorBench includes CCS to test posterior recovery in a realistic sparse-well setting, with ESMDA retained as a CCS-specific baseline.

B

Data Generation Details

B.1

Rejection sampling

To evaluate posterior-generating inverse solvers, we establish high-quality reference posterior samples using an accelerated rejection sampling scheme. Since querying the numerical forward solver sequentially during inference is computationally prohibitive, particularly for complex physical systems, we adopt an offline-to-online sampling strategy that leverages spatial invariances. Offline Prior Pool Generation. We first construct a massive offline dataset, denoted as the prior pool Dpool = {(x(i) , y(i) )}N i=1 , where N denotes a sufficiently large, task-specific pool size. Here, x(i) ∼ p(x) represents a generalized input physical parameter field (e.g., permeability, forcing terms, scattering coefficients, or geological models) sampled from the prior distribution, and y(i) = F(x(i) ) is the corresponding physical system state obtained via the forward numerical solver F.

The configuration of the prior, the numerical solver F, and the computational backend depend heavily on the specific partial differential equation (PDE) task: • Darcy Flow and Poisson Equation (JAX-accelerated): The input priors are constructed based on Gaussian Random Fields (GRFs). For Darcy flow, a binary field is generated by thresholding a GRF; for the Poisson equation, continuous GRFs with varying length scales and smoothness parameters (τ, α) are used to ensure a diverse distribution of spatial frequencies. Because these setups allow for highly parallelizable synthetic generation, both solvers are explicitly optimized using JAX, enabling massive batched execution on GPUs. • Light Transport and Carbon Capture and Storage (CCS): Unlike strictly synthetic GRF setups, these tasks involve highly complex physical structures (such as varying scattering 19

media for light transport and realistic geomodels for CCS). The respective prior pools are generated using their domain-specific, high-fidelity physical simulators, which are necessary to accurately capture the complex, non-linear forward dynamics governing light scattering or multiphase fluid flow. Observation Operators (H). During the online sampling phase, we define a target pair (xgt , ygt ) and obtain simulated observations ogt = H(ygt ). We denote i and j as the discrete spatial indices across the computational grid. To test the benchmark solvers across measurement regimes, we consider three distinct observation operators H tailored to different physical scenarios: 1. Sparse Random Observation (Darcy, Poisson): Sensors are randomly scattered across the spatial domain. The operator is defined as Hsparse (y) = {y(im , jm )}M m=1 , where M is the total number of sensors and (im , jm ) denotes the discrete spatial coordinates of the m-th. 2. Low-Resolution Observation (Darcy, Poisson, Light Transport): The system provides a coarse-grained view of the full solution, modeled via an average pooling operation: Hlow-res (y) = AvgPool(y, s), which downsamples the original high-resolution field to a coarser grid of scale s × s (e.g., 16 × 16). 3. Column Observation (CCS): Consistent with well-log data in geophysics, observations are only available along specific vertical columns. The operator extracts data strictly along these vertical indices: Hcol (y) = {y(ic , j) | ∀j, c ∈ C}, where C is the set of observable column indices ic . Rejection Sampling and Importance Weighting. To efficiently obtain posterior samples from Dpool given ogt , we employ a rejection sampling scheme based on the observation discrepancy. For tasks defined on a Cartesian grid with isotropic physical properties, specifically Darcy Flow and the Poisson Equation, we further exploit their discrete rotational equivariance. For these tasks, we apply discrete spatial rotations Rk ∈ {0◦ , 90◦ , 180◦ , 270◦ } to each candidate pair (x(i) , y(i) ), effectively quadrupling the pool size without additional solver calls. In contrast, for Light Transport and CCS, samples are utilized in their original orientation to preserve task-specific physical constraints. A candidate (x(i,k) , y(i,k) ) is accepted if the Root Mean Square Error (RMSE) between its observation and the target observation falls below a predefined threshold ϵ:   2 1/2 1 RMSE(i,k) = H(y(i,k) ) − ogt ≤ϵ (16) |Ωobs | 2 where |Ωobs | denotes the total number of observation points. We continue the search until K valid posterior samples (e.g., K = 100) are collected for the given target.

Finally, to account for the continuous nature of the posterior probability, we assign an importance weight w(i,k) to each accepted sample based on a Gaussian likelihood formulation. The unnormalized weights are computed as:  2  1  (i,k) (i,k) w̃ = exp − 2 RMSE (17) 2σ To ensure the likelihood properly reflects the acceptance criterion, we set the standard deviation to σ = 13 ϵ, ensuring that samples near the acceptance boundary ϵ are appropriately down-weighted. P The final weights are normalized such that w̃(i,k) = 1, yielding a weighted empirical posterior distribution that rigorously approximates the true Bayesian posterior. B.2

Carbon Capture and Storage

The CCS benchmark uses a synthetic dataset of supercritical CO2 injection into a radially symmetric deep saline aquifer. Each realization pairs a heterogeneous permeability field m ∈ R64×200 with the corresponding CO2 saturation field s = F (m) ∈ R64×200 after 30 years of continuous injection. Governing equations. The forward model solves the conservation of mass for a two-phase (CO2 – water) system in porous media. For phase α ∈ {w, g} (water and gas), the mass balance reads ∂ (ϕ ρα Sα ) + ∇ · (ρα uα ) = qα , (18) ∂t 20

Table 5: Rejection-sampling hyperparameters for Darcy flow and Poisson source recovery. Parameter

Darcy flow

Shared settings Resolution Threshold (η) Noise (σ) Random sensors Prior specification Smoothness (α) Correlation length (τ )

Poisson

128 4 × 10−4 1.33 × 10−4 128 2.0 3.0

1.5 ∼ 3.0 2.0 ∼ 4.0

where ϕ is porosity, ρα is density, Sα is saturation, and uα is the Darcy velocity krα K uα = − (∇Pα − ρα g). (19) µα Here K is the absolute permeability tensor (the unknown field), krα is relative permeability, µα is viscosity, and Pα is phase pressure. The system is closed by the saturation constraint Sw + Sg = 1 and the capillary pressure relation Pc = Pg − Pw . Simulation domain. The domain represents an infinite-acting aquifer with a radius of 100 km and a thickness of 135 m, discretized on a 2D radial grid (64 depth × 200 radial cells). No-flow conditions are enforced at the top and bottom caprock boundaries. CO2 is injected at a constant rate of 0.36 Mt/year for 30 years through a single vertical well. Geomodel prior. Permeability realizations are generated from geostatistical priors using sequential Gaussian simulation (SGeMS) [28]. The geostatistical hyperparameters (mean permeability µkr ∼ U[10, 500] mD, standard deviation σkr ∼ U[1, 500] mD, radial correlation length ∼ U[359, 35,900] m, and vertical correlation length ∼ U[14, 56] m) are sampled independently for each realization, creating a diverse prior that spans a wide range of geological scenarios. The forward simulations are performed with the industry-standard reservoir simulator ECLIPSE (e300) [29]. The dataset comprises 12,000 training pairs and 1,390 test pairs. Further details are provided by [27]. Observation Pattern Observations mimic realistic well-monitoring data: vertical column measurements are collected at two locations corresponding to an injection well (x = 0) and a monitoring well (x = 50, approximately 491 m away). Each column provides 64 saturation measurements along the depth axis, yielding 128 observed values out of 12,800 total grid cells (1% spatial coverage). Gaussian observation noise with σobs = 0.04 is added to the ground-truth saturation values at the observed locations. Reference Posterior Construction The reference posterior for each test case is constructed via rejection sampling from a pool of 2 million prior geomodel samples. Each candidate is evaluated through a pre-trained Local Neural Operator (LNO) surrogate [71] that maps permeability to saturation, and the observation mismatch is computed at the monitored well locations. Candidates are accepted with probability proportional to the Gaussian likelihood:   ∥Mobs ⊙ (Lϕ (m) − yobs )∥2 L(m) = exp − , (20) 2 · |M 2σobs obs |

where Mobs is the binary observation mask and |Mobs | is the number of observed points. This procedure yields approximately 26,000 accepted samples per case (acceptance rate ∼1.3%), providing a dense empirical approximation to the true posterior. The use of the neural-operator surrogate rather than the full simulator makes this large-scale rejection sampling computationally feasible.

C

Reference Posterior Validation

Each metric is an estimator evaluated against a finite reference ensemble and therefore carries a variance component attributable to the reference set alone. We isolate this component by constructing 21

two independent reference ensembles, GT1 and GT2 , at equal sampling budget and scoring the same solver outputs against each ensemble. Under reference convergence, GT1 and GT2 are exchangeable and should agree up to Monte Carlo error. For each metric m, we compute   mGT1 mGT2 rm = max , . (21) mGT2 mGT1 We report the worst-case ratio over all evaluated solvers in Table 6. A value close to 1 indicates that the metric is insensitive to the particular finite reference ensemble used. Table 6: Reference-ensemble consistency check. Each entry reports the worst-case ratio rm over evaluated solvers when the same solver outputs are scored against two independent reference ensembles. Task Poisson Darcy Flow Light Transport Carbon Capture and Storage

Mean Rel. L2

Std Rel. L2

MMD

SWD

Spectral Rel. L2

1.0295 1.0272 1.0132 1.0767

1.0198 1.0269 1.0141 1.0505

1.0112 1.0196 1.0090 1.0366

1.0431 1.0490 1.0302 1.0738

1.0674 1.1146 1.0589 1.1266

All ratios are close to 1, and all non-spectral metrics remain below 1.08. These values are worst-case ratios rather than average-case ratios, making the check conservative. The results suggest that the reference ensembles are sufficiently converged for the reported solver comparisons and make the empirical evaluation noise floor transparent.

D

Additional Results

D.1

Guidance Weight Calibration

Figure 8 reports the full guidance-weight sweep for FunDPS on Poisson source recovery across posterior thresholds ϵ ∈ {3, 4, 5, 6, 7, 8} × 10−4 and guidance weights λ ∈ [0, 1000]. Each panel visualizes one posterior metric as a function of λ under each threshold. Across metrics, very small guidance weights under-condition on the observation, whereas overly large weights increasingly distort posterior calibration. The location of the optimum generally shifts toward smaller λ as the threshold increases, consistent with the interpretation that noisier observations should exert weaker likelihood guidance. However, the best guidance strength is not metric-invariant: posterior-standarddeviation error is often minimized at substantially smaller λ than posterior-mean error, while MMD, SWD, and spectral error select intermediate regimes. This separation indicates that a single scalar guidance coefficient cannot simultaneously optimize both mean tendency and uncertainty calibration. Figure 9 shows case-level error fields for a representative Poisson source-recovery case at three guidance weights, λ ∈ {25, 250, 1000}. These fields provide a spatial diagnostic of the two main failure modes observed in the sweep. With weak guidance, the posterior mean retains coherent residual structure because the sampler remains insufficiently conditioned on the observation. With overly strong guidance, the correction becomes spatially uneven and can introduce localized overcorrection, indicating that stronger observation consistency alone does not guarantee a calibrated posterior field. The standard-deviation panels show the corresponding effect on uncertainty: both under-guidance and over-guidance can distort the spatial distribution of posterior spread.

22

GT threshold threshold 3e-4 threshold 4e-4 threshold 5e-4

GT threshold

threshold 6e-4 threshold 7e-4 threshold 8e-4

Std relative L2

Mean relative L2

1.0 0.8 0.6 0.4

threshold 3e-4 threshold 4e-4 threshold 5e-4

1.5

1.0

0.5

0.2 0

100 200 300 400 500 600 700 800 900 1000

0

100 200 300 400 500 600 700 800 900 1000

Guidance

Guidance

(a) Posterior mean error 0.8

(b) Posterior standard-deviation error GT threshold

GT threshold threshold 3e-4 threshold 4e-4 threshold 5e-4

0.7

threshold 6e-4 threshold 7e-4 threshold 8e-4

threshold 3e-4 threshold 4e-4 threshold 5e-4

40

0.6

SWD

MMD

threshold 6e-4 threshold 7e-4 threshold 8e-4

0.5

threshold 6e-4 threshold 7e-4 threshold 8e-4

30 20

0.4 10

0.3 0

0

100 200 300 400 500 600 700 800 900 1000

0

100 200 300 400 500 600 700 800 900 1000

Guidance

Guidance

(c) MMD

(d) SWD

Spectral relative L2

GT threshold threshold 3e-4 threshold 4e-4 threshold 5e-4

0.6

threshold 6e-4 threshold 7e-4 threshold 8e-4

0.4

0.2

0

100 200 300 400 500 600 700 800 900 1000

Guidance (e) Spectral error

Figure 8: Metric-specific guidance-weight sweeps for FunDPS on Poisson source recovery. Each panel plots one evaluation metric as a function of the guidance weight λ, with separate curves for posterior thresholds ϵ ∈ {3, 4, 5, 6, 7, 8} × 10−4 . The curves expose both the broad decrease in preferred guidance strength as observation noise increases and the disagreement between metricspecific optima at the same threshold.

23

GT threshold 3e-4

4e-4

5e-4

6e-4

7e-4

8e-4

0.6

25

Guidance

0.4

0.2

250

0.0

?0.2

?0.4

1000

?0.6

GT threshold 3e-4

4e-4

5e-4

6e-4

7e-4

8e-4

0.15

25

Guidance

0.10

0.05

250

0.00

?0.05

?0.10

1000 ?0.15

Figure 9: Case-level delta fields for the Poisson guidance sweep. The top panel visualizes posteriormean error fields and the bottom panel visualizes posterior-standard-deviation error fields for the same representative case at representative weak, intermediate, and strong guidance weights. These spatial diagnostics illustrate the localized errors induced by under-guidance and over-guidance.

24

D.2

Metric-Pair Diagnostics

We visualize pairwise relationships among the five main posterior metrics. Each panel reports raw metric values on log-log axes; consistent trends indicate agreement between two metrics, while scattered panels highlight metric-specific failure modes. Mean vs Std

Mean vs MMD

100

10 1

10 1

100

Mean vs SWD

100

Mean vs Spectral

101

100 101

10 1

10 2 10 1

10 1

100

Std vs MMD

100

Std vs SWD

101

100

100

25

101

Std vs Spectral

MMD vs SWD

100

101

10 1

10 2 100

101

MMD vs Spectral

101

100

100

10 1

10 1

10 2

10 2

SWD vs Spectral

101

Figure 10: Raw-value metric-pair diagnostics for the five posterior metrics.

26

D.3

Standard Deviations for Main Result Table Table 7: Standard deviations across cases for the main benchmark metrics in Table 2.

Task

Method

Mean error

Std error

MMD

SWD

Spectral error

Darcy flow inversion

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff

0.1395 0.0339 0.0156 0.0096 0.0172 0.0191 0.0395

0.3115 0.3445 0.0605 0.0660 0.1412 0.0584 0.0457

0.1027 0.0352 0.0241 0.0330 0.0242 0.0283 0.0474

14.2912 3.6436 0.8464 0.5116 1.3175 1.5804 1.6540

0.1008 0.1899 0.0820 0.0653 0.0822 0.0446 0.0379

Poisson source recovery

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff

1.4779 0.1754 0.1693 0.0931 0.1825 0.0546 0.1372

0.6601 0.3814 0.0640 0.0469 0.1459 0.0311 0.3464

0.0118 0.0247 0.0477 0.0472 0.0401 0.0265 0.0319

7.2965 4.4578 1.3222 0.2433 0.4421 0.3378 4.9226

21.5409 5.1581 0.1573 0.0567 0.1710 0.0946 3.1165

Carbon capture and storage

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff

0.7666 0.2960 0.1093 0.7437 0.1608 0.3706 0.1205

2.0413 1.8626 0.2839 1.7852 0.2593 0.5733 0.1492

0.0637 0.1371 0.0795 0.1121 0.0995 0.1083 0.0980

8.1664 4.8350 2.0215 4.0355 3.0374 4.5436 4.1860

257.6435 48.8525 0.9662 56.4968 1.3179 7.7058 0.2191

Light transport material inference

ECI ES-MDA FunDPS Fun-DDPS DiffusionPDE DDIS FunDiff

0.0278 0.0199 0.0158 0.0170 0.0166 0.0193 0.0219

0.0691 0.0235 0.0238 0.0201 0.0201 0.0327 0.0318

0.0362 0.0366 0.0235 0.0234 0.0231 0.0334 0.0291

1.4149 0.4211 0.6200 0.5861 0.5411 0.6891 0.7305

0.0905 0.0302 0.0216 0.0891 0.0873 0.0895 0.0527

E

Experiment Details

This appendix records the training and evaluation protocol for the baseline inverse solvers used in PosteriorBench. All method–dataset pairs are run through the unified repository on NVIDIA B200 GPUs, and measured under the same evaluation pipeline. E.1

Baseline Training Protocol

For learned baselines, the training data for each task consists of paired prior samples and forward observations generated by the task-specific simulator or surrogate described in Section B. Each method is trained on the official training split for that task and is evaluated only on held-out benchmark cases. When a method requires a learned prior, score model, flow model, neural operator, or differentiable surrogate, that component is fit without access to the reference posterior samples used for evaluation. Baseline training and inference configurations are summarized in Section E.3. E.2

Baseline Evaluation Protocol

At test time, each baseline receives the same observation for a given benchmark case and returns an ensemble of posterior samples. The reported metrics compare this ensemble with the corresponding reference posterior for that case, using posterior mean error, posterior standard deviation error, maximum mean discrepancy, sliced Wasserstein distance, and radially averaged power-spectrum 27

error. For fair distributional comparison, methods with larger generated ensembles are subsampled to the common evaluation size used by the task before computing metrics. All subsampling and metric computation are performed in the physical parameter space of the inverse problem. E.3

Baseline Solver Configurations

The remaining subsections summarize the eight posterior solvers used in the benchmark. For traceability, learned methods are reported with training and, where applicable, inference tables, while ES-MDA is reported by its inference-time assimilation parameters. The tables focus on solver and architecture parameters; field choices, grid resolutions, training budgets, observation operators, and normalization constants are not treated as method hyperparameters. E.3.1

FunDPS

FunDPS [20] trains a joint diffusion prior over the unknown field together with the corresponding observable field. During posterior sampling, DPS guidance is applied to the observable channel while the joint state is sampled. Table 8: FunDPS training configuration by benchmark task. Category

Parameter

Prior

Architecture Learning rate LR ramp-up EMA half-life Dropout Channel base Channel multipliers Attention resolutions UNO blocks Operator rank

Darcy

[8]

Poisson

CCS

Light transport

ddpmpp-uno 1.0×10−4 5M images 0.5M images 0.13 64 [1, 2, 4, 4] [16] [16] 4 0.1

[16]

Table 9: FunDPS inference configuration by benchmark task. Category

Parameter

Sampling

Initial latent family Reverse steps σmin σmax ρ

Guidance

E.3.2

Loss Field weight

Darcy

0.002 80

4000

Poisson

CCS

Light transport

RBF random field 500 0.01 0.002 0.002 10 80 80 7 1600

MSE 300

3000

Fun-DDPS

Fun-DDPS [27] decouples the target-field diffusion prior from the forward map used for likelihood guidance. The observation map is represented by a task-specific FNO surrogate. E.3.3

DDIS

DDIS [22] also separates the target prior from the differentiable forward map, but uses a DAPS-style decoupled sampling update with annealing, diffusion, and Langevin correction phases. The forward map is represented by a padded FNO surrogate in the canonical configurations. E.3.4

DiffusionPDE

DiffusionPDE [17] trains a joint grid-based score model over the unknown field and the associated physical field. In our unified runs, observation guidance is applied to the observable channel during 28

Table 10: Fun-DDPS training configuration by benchmark task. Category

Parameter

Prior

Architecture Batch size Learning rate Attention resolutions Channel base Channel multipliers UNO blocks

Surrogate

Architecture Fourier modes Hidden channels Layers Epochs Batch size Learning rate

Darcy

[8]

[16, 16]

Poisson

CCS

Light transport

ddpmpp-uno 32 1.0×10−4 [16] [16] 64 [1, 2, 4, 4] 4 [16, 16]

[16]

FNO [16, 32] 48 4 100 32 1.0×10−3

[16, 16]

Table 11: Fun-DDPS inference configuration by benchmark task. Category

Parameter

Darcy

Sampling

Initial latent family Reverse steps σmin σmax ρ

RBF random field 500 0.002 80 7

Guidance

Loss Surrogate model Field weight

MSE FNO 10.0

1000

Poisson

CCS

100

Light transport

100.0

Table 12: DDIS training configuration by benchmark task. Category

Parameter

Prior

Architecture Batch size Learning rate Attention resolutions Channel base Channel multipliers UNO blocks

Surrogate

Darcy

[8]

Architecture Fourier modes Hidden channels Layers Epochs Batch size Learning rate

Poisson

CCS

ddpmpp-uno 32 1.0×10−4 [16] [16] 64 [1, 2, 4, 4] 4 FNO-pad [64, 64] 64 4 500 40 1.0×10−4

29

Light transport

[16]

Table 13: DDIS inference configuration by benchmark task. Category

Parameter

Sampling

Initial latent family

Annealing

Steps σmax σmin ρ

Diffusion

Correction steps Correction σmin Correction ρ

Langevin

Steps Learning rate LR minimum ratio LR ρ η τ

Guidance

Darcy

Poisson

CCS

Light transport

RBF random field 200

100

10 0.01 7

100

100

5 0.001 7 10 1.0×10−4

35 7.0×10−5

0.1

0.25

Loss Surrogate model Field weight

10

2.0

20 1.0×10−4 0.01 1.0 0.1 0.001

20 1.0×10−4 0.1

MSE FNO-pad 0.25

1.0

reverse diffusion. We modified the original implementation to fix the normalization issue and support batched inference, which actually improves performance over the stock implementation. Table 14: DiffusionPDE training configuration by benchmark task. Category

Parameter

Prior

Architecture Batch size Learning rate LR ramp-up EMA half-life Dropout

Darcy

Poisson

CCS

Light transport

ddpmpp 32 1.0×10−3 10M images 0.05M images 0.13

Table 15: DiffusionPDE inference configuration by benchmark task.

E.3.5

Category

Parameter

Sampling

Initial latent family Reverse steps σmin σmax ρ

Guidance

Loss Observation fraction Late-observation multiplier Field weight

Darcy

Poisson

CCS

Light transport

Gaussian 2000 0.002 80 7

40000

20000

L2 0.8 0.1 1000

400

FunDiff

FunDiff [21] represents target functions with function autoencoders and learns a conditional latent rectified flow with a DiT backbone. Unlike guidance-based diffusion samplers, FunDiff conditions on the observation representation directly and does not use a scalar likelihood-guidance weight at inference. 30

Table 16: FunDiff training configuration by benchmark task. Category

Parameter

Darcy

Poisson

FAE

Encoder patch size Embedding dimension Latent tokens Depth Attention heads Max steps Batch size Peak learning rate Target query count

16 × 16

16 × 16

DiT

Embedding dimension Depth Attention heads Max steps Batch size Peak learning rate

CCS

Light transport

8 × 25 256 64 8 8 100k 16 1.0×10−3 4096

8×8

384 16 8 117188 128 1.0×10−3

Table 17: FunDiff inference configuration by benchmark task.

E.3.6

Category

Parameter

Darcy

Poisson

CCS

Light transport

Sampling

Flow steps Decode chunk

20 1024

20 1024

100 2048

20 1024

ECI-sampling

ECI-sampling [30] trains an FNO-based flow-matching model for the joint task state. At inference, observation consistency is imposed through operator-specific conditioning, including hard replacement for sparse or low-resolution observations where supported by the method profile. Table 18: ECI-sampling training configuration by benchmark task.

E.3.7

Category

Parameter

Flow model

Backbone Fourier modes Hidden channels Layers Embedding channels Base noise

Optimization

Max steps Batch size Learning rate Weight decay Adam (β1 , β2 )

Darcy

Poisson

CCS

Light transport

FNO [32, 32] 128 6 32 Matern 58594 256 3.0×10−4 0 (0.9, 0.999)

ES-MDA

ES-MDA [24, 25] is a non-neural ensemble smoother baseline. It does not train a generative model; instead, it updates a task-specific ensemble drawn from the prior using repeated Kalman-style assimilation steps. 31

Table 19: ECI-sampling inference configuration by benchmark task. Category

Parameter

Sampling

Steps Mixture count nmix Resampling interval

Darcy

Poisson

CCS

Light transport

5 1

1 1

200 1 5

1 5

Table 20: ES-MDA inference configuration by benchmark task. Category

Parameter

Darcy

Poisson

CCS

Light transport

Ensemble

Ensemble size

128

512

256

1024

Assimilation

Number of updates Inflation factors αk

2 [2, 2]

1 [1]

4 [4, 4, 4, 4]

1 [1]

Observation error

Std. model Std. multiplier Minimum std Jitter

prior-predictive standard deviation multiplier 10 100 2 1 1.0×10−6 1.0×10−8

E.3.8

FNO with MC Dropout

The MC-dropout baseline trains a direct inverse FNO that maps masked observations to the target field. At inference, dropout layers remain active and posterior samples are obtained from repeated stochastic forward passes of the same trained model. Table 21: FNO with MC Dropout training configuration by benchmark task. Category

Parameter

Model

Backbone Fourier modes Hidden channels Layers Channel-MLP dropout Domain padding

Optimization

Darcy

Poisson

[32, 32]

[32, 32]

Max images Batch size Learning rate Weight decay Scheduler

CCS

FNO [16, 64] 128 6 0.2 0.1, one-sided 15M 64 5.0×10−4 2.0×10−2 One-cycle

32

Light transport [32, 32]