Diffusion posterior sampling for simulation-based inference in tall data settings
Authors: Julia Linhart
Gabriel Cardoso
Alexandre Gramfort
Sylvain Le Corff
Pedro L. C. Rodrigues
Venue: Accepted by TMLR ()

Abstract: Identifying the parameters of a non-linear model that best explain observed data is a core task across scientific fields. When such models rely on complex simulators, evaluating the likelihood is typically intractable, making traditional inference methods such as MCMC inapplicable. Simulation-based inference (SBI) addresses this by training deep generative models to approximate the posterior distribution over parameters using simulated data. In this work, we consider the tall data setting, where multiple independent observations provide additional information, allowing sharper posteriors and improved parameter identifiability.
Building on the flourishing score-based diffusion literature, F-NPSE (Geffner et al., 2023) estimates the tall data posterior by composing individual scores from a neural network trained only for a single context observation. This enables more flexible and simulation-efficient inference than alternative approaches for tall datasets in SBI. 
However, it relies on costly Langevin dynamics during sampling. We propose a new algorithm that eliminates the need for Langevin steps by explicitly approximating the diffusion process of the tall data posterior. Our method retains the advantages of compositional score-based inference while being significantly faster and more stable than F-NPSE. We demonstrate its improved performance on toy problems and standard SBI benchmarks, and showcase its scalability by applying it to a complex real-world model from computational neuroscience.

Published in Transactions on Machine Learning Research (02/2026)
Diffusion posterior sampling for simulation-based inference in tall data settings
Julia Linhart
julialäibart@iinrio.f
Université Paris-Saclay
Inria, CEA
Gabriel Victorino Cardoso
CMAP, Éoole Polglechnique
Institut Polytechnique de Poris
Alexandre Gramfort
Université Paris-Saclay Inria, CEA
Sylvain Le Corff
sylvain. le corfflsorbome--uniierrssseeee
LPSM, Sorbonne Université UMR CNRS 8001
Pedro L. C. Rodrigues
pedro.rodrrigesshhnrii..
Université Grenoble Alpes
Inria, CNRS Grenoble INP, LJK
Reviewed on OpenReview: https: // openreview. net/ forum? id= cobfoS60yo
Abstract
Identifying the parameters of non-linear model that best explain observed data is task across scientific fields. When such models rely complex simulators, evaluating the likelihood is typically intractable, making traditional inference methods such as MCMC inapplicable. Sinulation-based inference (SBI) addresses this by training deep gvnerative models to approximate the posterior distribution over parameters using simulated data. In this work, we consider the tall data setting, where multiple independent observations provide additional information, allowing sharper posteriors and improved parameter identifiability. Building the flourishing score-based diffusion literature, F-NPSE. (Geffner al., 2023) estimates the tall data posterior by composing individual scores from neural network trained only for a single contert observation. This enables more flexible and simulation- efficient inference than alternative approaches for tall datasets in SBI. However, it relies on costly Langevin dynamics during sampling. We propose a new algorithm that climinates the need for Langevin steps by explicitly approximating the diffusion process of the tall data posterior. Our method retains the advantages of compositional score-based inference while being significantly faster and more stable than F-NPSE. We demonstrate its improved performance toy problems and standard SBI benchmarks, and showcase its scalability by applying it to a complex real-world model from computational neuroscience.

---

Published in Transactions on Machine Learning Research (02/2026)
=
n=
64
0
analytic
estimated
0
0
6
Figure 1: The posterior distribution of a model with a Gaussian simulator and Gaussian prior concentrates around the true parameter l as the number of observations z p(x: (*) increases. The analytic posterior is compared to the posterior estimated with score-based proposal (Algorithm GAUSS).
1 Introduction
Inverting non-lincar models that describe natural phenomena is fundamental problem in many scientific domains (Gonçalves et al.) 2020 Dax et al., 2023). We adopt Bayesian perspective, where the goal is to infer the posterior distribution p(0 x) that relates input parameters E IRm to output observations E IRd. Sampling from this posterior becomes particularly challenging when the model's outputs are generated via complex simulations (e.g. based on non-linear stochastic differential equations). In such settings, the likeli- hood p(x () is often intractable, making classical Bayesian approaches such as Markov Chain Monte Carlo (MCMC) (Hastings 1970) either inapplicable or computationally prohibitive. Simulation-basee (SBI), also known as likelihood-free inference, circunvents this limitation by relying on model simulations instead of explicit likelihood evaluations. Here, novel deep learning techniques can be used to accurately approximate arbitrarily complex posterior distributions (Cranmer et al. 2020 The standard procedure as- X(0) encoding scientific knowledge about plausible parameter values and uses the simulator to draw samples from the joint distribution:
inference
sumes a prior distribution
θ X(0). X; p(r O.), (i,, p(0.x), i= 1...,,,,
with N, the simulation budget. From this simulated dataset, a neural network is trained to estimate statistical quantities of interest defined for every (0,z). Once trained, it can be evaluated in any new observation x= approximate the posterior distribution p(0 x=)), a property known amortized inference, Existing SBI methods differ in the target of the neural estimator: Neural Posterior Estimation (NPE) learns the posterior directly (Greenberg
a.. 2019 Neural Likelihood Estimation (NLE) approximates the likelihood (Papamakarios et al. 2019), and Neural Ratio Estimation (NRF) estimates the likelihood- to-evidence ratio (Hermans et al.] 2020). While NPE provides posterior samples directly via conditional
normalizing flows Papamakarios et al. 2021 NLE and NRE require additional MCMC sampling.
In this work, we consider an extension of the above Bayesian inference framework the tall data setting (Bar- denet et al. 2015), where multiple i.i.d. observations z (#i,..., .., =4) available The corresponding tall data posterior is expected to provide more precise information about how to invert the simulator-mode with increasingly sharper posterior densities the number of observations grows (see Fig
p(0 ximm) X(0)
Il j=1
p(0 xt))
While this framowork crucial for practical applications, the extension of existing SBI algorithms to tall data settings remains chall challenging and no satisfactory solution has yet been proposed. In Rodrigues el al. (202]
I.i.d.: independent and identically distributed. example, simulations (#)) ~ p(a 8*) that generated with the set of simulator parameters i.d.d,, conditionally 8*.
2

---

Published in Transactions on Machine Learning Research (02/2026)
the authors merge fixed number of extra observations via deepset (Zaheer et al. 2017 and fall back to NPE
trained on an augmented dalasel. DD, [(O,, X1,:=)), with simulations X (X...... Xin)) ~ p(1 Ox) per prior sample O; ~ X(0). While NPE leverages expressive normalizing flows to model the posterior directly, it imposes architectural constraints and requires fixed input dimensions on the neural network, limiting efficient and flexible conditional data modeling. This leads to two major drawbacks: (1) the potentially heavy simulation cost, as the simulation budget scales with the number of observations (Naae N, x n) and (2) the lack of flexibility at inference time when more observations become available, new model tailored to the new context size must be retrained from scratch. In Hermans et al. (2020), the authors show how the amortization of NRE allows to handle the tall data setting without reguiring new simulations or retraëning. By factorizing over the set of multiple observations, they reformulate the "tall data ratio" as the composition of individual ratios from a network trained only for a single context observation. An equivalent extension can be done for NLE (Geffner et al. 2023). However, these approaches still require MCMC sampling to obtain posterior samples, which remains a computational bottleneck and requires carcful hyperparameter tuning for convergence guarantees to hold.
Recent developments in score-based generative modeling (SBGM) (Ho et al. 2020 Song et al. 2021b) offer promising alternative. These methods target the gradient of the log probability- known as the score function used to reverse a diffusion process that transforms random noise into structured samples from the target distribution. A more detailed overview is provided in Section 2.1 SBGM rivals state-of-the-art generative modeling approaches such as generative adversarial networks (GANs) (Goodfeloow et al. 2014) and normalizing powssaaaaaaaaaiios et al., 2021 on challenging tigh-dimensional datasets, without requiring adversarial training special network architectures. This has led to the adoption and increasing popularity of scoro-based methods in SBI, introduced as Neural Postcrior Score Estimation (NPSE) by SSharrock et al. (2022). NPSE is now among the most flexible SBI algorithms and achieves state-of-the-arl performance, especially when paired with transformer architectures (Gloeckler et al. 2024h.
Crucially, NPSE dircctly targets the postcrior (unlike NLE/NRE), but also supports compositional inference (unlike NPE). Specifically, from the factorization in Equation ( the "tall posterior score" can be constructed as a sum of individual posterior scores. F-NPSE, recently proposed by Geffner et al. (2023), leverages this idea by composing scores from network trained on single-observation contexts. Like Hermans et al. (2020)
for NRE, it exploits the amortization of NPSE to approximate the tall data posterior without requiring
an augmented dataset, and naturally adapts to variable context sizes. However, a major drawback of this approach is that the diffusion process from the composed score is unknown, reintroducing the need of MCMC to sample from the posterior via an annealed Langevin procedure. We discuss this further in Scction 2.3
In this work, we propose a new sampling algorithm that relains the amortized and compositional benefits of F-NPSE, while eliminating the need for Longevin dynamics. By explicitly approximating the diffusion process of the tall-data. posterior, our method enables faster and more stable inference. We demonstrate its superiority over F-NPSE- in both, accuracy and numerical two Gaussian toy models for which all quantities of interest are known analytically, and multiplc examples from the SBI benchmark. Finally, we apply our method to invert a complex model from computational neuroscience, demonstrating its scalability to challenging real-world problems.
thhliityythruuhh several numerical experiments:
Section 2 reviews basee-based generative models, their application to SBI, and the F-NPSE method. Section 3 introduces the mathematical loundations of our approach and presents the proposed algorithms. Section reports experimental results, and Section 5 concludes with discussion of findings and future perspectives.
2 Background
2.1
Score based generative models (SBGM)
The goal of SBGM is to estimate an unknown target data distribution Pdata from i.i.d. samples with the help of a forward diffusion process that adds noise to the training data. The main idca is to approximately sample the target distribution by solving the associated backward diffusion process. This procedure requires estimating the score functions of the diffused data. distributions for different levels of added noise.
3

---

Published in Transactions on Machine Learning Research (02/2026)
Forward diffusion process (data noise). Formally, the idea is to construct sequence of distribu- tions (Et)r that defines increasingly noisy versions the larget distribution Pdata, by convolving the data distribution with known forward kernel 9/jo
/ (ejo(0r(00))aaan (66) 1 E [1.
Po = Pdatas (c(),)
(2)
In this work, we focus on the variance preserving (VP) framework (Ho et al. 2020 Yang et al. 2023 in which
the forward kernel defined by (1))(0 N (0,0 var0o, vmmm with (at E [0, 1|² decreasing
sequence of time-dependent scale factors and 1/2 1 a determines the amount of added noise 2
Score estimation. Following [Vincent 2011
learn the scores Ve, (..): of cach noisy distri- bution pt via neural network s5, by minimizing the denoising score-matching (DSM) loss
DOSM (Ø) = =- r 0r [Isa(Or,)) t) Ve, log ge (0, (0.)
(3)
where o? is a weighting function. This loss is completely tractable, as the forward kernel can be easily sampled from and it's score function is available in closed form: given a sample from the training set Oo ~ Pdata and C ~ NOO Im), we have O₁ = Va + Vuree ~ 9:j0 and Ve, log (:o(O,
(θ, Va)
Backward diffusion process (noise - data). Once the score estimator (b,,,, t) = Ve, Po is trained, the goal is to draw backwards starting from the noise distribution (0 Im) approximating pr to obtain samples approximately distributed according to Pdata at time == O:
Or ~ N10, Im) (0(.T.T)
θr-1 (0 1,7-1)
(0 On.
(4)
This represents the "gencrative" part in SBGM: the score directly drives the reverse dynamics, turning noisy inputs into progressively cleaner samples. Different ways to perform this backward sampling include the use
of annealed Langevin dynamics (Song & Ermon 2019), stochastic differential equations (Song et al., 2021b).
or ordinary differential equations (Karras al. 2022
In this work, we follow the approach proposed in Song al. 2021a). which yields the denoising diffusion implicit models (DDIM) sampler. DDIM provides a deterministic update rule to directly map Pt. to pt-1 using the learned score sy from 3 Formally, this consists in Gaussian transition kernels of the form
g (o.--1((0~1100)) N (0..1: Im))
(5)
where (or E (0, 1/2 - and the mean depends explicitly on the score network:
1 (, + v,sp((r.tt)
(0.) = Vaa
(6)
Finally, composing these kernels yields the DDIM backward Markov chain:
T
Po.o.((0..1 = (O II
yo. (0-111)00-11))
Note that this only approximates the true reverse diffusion dynamics: the exact backward kernel involves intractable expectations over Oo, which DDIM replaces with acore-based estimates (as shown in Appendix C Empirically, DDIM has been shown to produce high-quality samples with fewer steps than alternative stochas- tie samplers, making it particularly attractive settings where efficiency and stability are crucial.
Under weak conditions (Puaza has finite second moment), one can show that (T.(r]|VV Im)) as grows.
4

---

Published in Transactions on Machine Learning Research (02/2026)
2.2 Neural Posterior Score Estimation (NPSE)
Neural Posterior Score Estimation (NPSE) Sharrock et al. 2022) adapls score-based generative modeling to simulation-based inference. The goal is now to generate new samples from the posterior Pdata := p(0[x). This achieved by extending the score network So(0,, z. t) to take I as an additional input, trained to approximate the score Ve, log (o) Training follows the same denoising score matching principle used in unconditional settings, but now minimizes the loss expectation over both and æ:
T (N) = Σ hii.eee...rrppee....e.e.eee (150 (O,, - Ve, log (0, x) t=1
(8)
The key insight is that this only requires to draw training pairs from the joint distribution p(0, x) p(0)p(z]0), which is readily available in SBI via simulation. Convergence of this objective to the true conditional score function has been shown by Batzolis et. al. (2021). Once trained, NPSE enables amor- tized inference by generating samples from the posterior p(0] x for any new observation xt by plugging the learned score &6 into any score-basod sampler (c.g. DDIM).
2.3 Factorized Neural Posterior Score Estimation (F-NPSE)
The F-NPSE method proposed in Geffner et al. (2023) defines a sequence of distributions
Vo, log De (0, xi:.) (1 n1 V), log X(O.) + S (0,, 1),
(9)
1 Note, however, that this sequence does not correspond to the true diffusion process, as defined in (2. i.e. er æōnn) # Pr(0 =)..), for >> 0. Its appeal lies in the compositional structure of the score: it enables inference from single amortized score network (here 8&), trained on individual observations, eliminating the need for augmented datasets or retraining. F-NPSE has shown to outperform competing SBI methods for tall data settings (based on NPE, NLE, or NRE), achieving the best trade-off between sample efficiency and error accumulation as n grows.
which, at / O, coincides with the score of the tall data posterior from Equation
The main drawback is sampling: since the composed score does not correspond to any known forward diffusion, deterministic samplers such DDIM camnot be applied. F-NPSE therefore relies annealed Langevin dynamic47] which is sensitive to hyperparameters and can require many iterations for convergence. This makes Langevin the bottleneck of F-NPSE, and motivates the need for alternatives that enable faster and more stable sampling.
3 Diffusion Posterior Sampling for tall data
In this section, we propose a new algorithm that approximately samples the tall data posterior, using only individual scores obtained for each observation from a previously trained amortized NPSE. The novelty of our method lies in the tractable computation of the scores associated with the diffussonn process of the factorized tall posterior, thereby eliminating the need of costly and unstable Langevin steps used in F-NPSE.
(Section 3.1 1
First, we derive a formula for the score of the diffused factorized tall data posterior.
We then show how compute this score, by introducing new approximation of the considered backward diffusion process.
(Section 3.2
(Section 3.3
We propose two different algorithms to efficiently implement this approximation.
The resulting approximate score can then directly be plugged into deterministic score- score-based samplers, such as DDIM, to infer the tall data posterior.
Langevin dynamics provide MCMC procedure approximate the backward diffusion process and sequentially gencrate samples each 2t-1 from ør (Song & Ermon 2019). It only requires access the score function, closed form expressions of the backward kernels deterministic samplers like DDIM.
5

---

Published in Transactions on Machine Learning Research (02/2020)
3.1 Exact computation of the tall data posterior score
Let zim (ai 24,) be i.i.d. observations. Our goal is to sample from the tall data posterior p(0 zim) via the DDIM backward Markov chain defined in Equation (34). while only relying on score estimate s0 t) Vo, (., x), i.c. the diffused posterior for a single contert observation x. To do so, we need a closed-form expression of the diffused tall data posterior score Ve, log (r iim) that writes as a function of the individual posterior scores Vo, 1o xf), for every j E [1,
Using the factorized expression from
1 in the definition of the diffuused tall posterior from 2 we can write
Pt (0, winn) = / p(0o xi
(10)
TI / (((0)) p(0 xz)) ( 2: (0:000(0))d0..
x
Given the diffused prior (0.) ((O)) (0.)0) )000 and the diffused individual posterions pc(O, xz), we now introduce the following backward transition kernels obtained via Bayes' rule:
d (0 0,) ((00))e and Toje (0) p(0o x))00(0.100) (0,)
(11)
(, x)
Rearranging terms in Equation 10 and using 11 yields
Pr wim) / (((000))0(011000)
p(o (0))00((01000)00
j=1
La x )p? (0.)
II P. zj),
with Lx(0,, vi...) = - ITy= (00100 )d06 The corresponding score writes
Vo, logp,(0) = -1 (
+ Σ Ve, log logp,(0(l2)) + Ve, log (0,,
(12)
The first two terms in the above equation are tractable: the prior score can be computed analytically in most case and the single posterior scores are approximated by evaluating the learned score model so x, t) at every zj. This leaves us with the last term: the score of (x ætm),, which involves the backward kernels of the prior and each individual posterior. It intractable and needs to be approximated. Note that this term is missing in the score formula from F-NPSE and is the reason why Langevin corrector steps required. More intuition on the influence of this correction term is given in the following section and Appendix B
3.2 Second order approximation of the backward diffusion process
In the previous section, we derived closed-form expression of the tall data posterior score. We now show how to compute it efficiently. The difficulty lies in the last term of 12), the score of Lx(0, xi.)) It involves integrating over all backward kernels, and cannot be computed in closed form. To handle this, we approximate these kernels using Gaussian distributions, following the Tweedie framework (Boys et al. 2023):
t (fo 0t) (fo: (1,x (o.) and (o xz)) = NNo:: pr, (0.)).
(13)
where the means and covariance matrices the ones of the backward processes respectively associated with the difflused prior pd and each diffused individual posterior Pr(O, xz)) Specifically, Boys et al. (2023) show
'See Appendix for the Gaussian and Uniform cases. If not analytically computable, the prior score can be learned the classifier-free guidance approach, at the same time the posterior score, shown in Appendix I..I

---

Published in Transactions on Machine Learning Research (02/2026)
that (0:) = E [Bu 0e] and ((()) Cou(Oo (b) are functions of the score of pt and its derivatives, which are computationally tractable. Hence, this choice has two advantages: (i) it preserves the local mean/ mean/variance variance structure of the true backward diffusion process, and (ii) it yields a tractable formula for the correction term.
In what follows, Lemma 3.1 shows that, under this approximation, the correction term Lx simplifies to a product of Gaussian distributions. After taking logarithms, this product becomes a tractable sum that integrates seamlessly into the score expression of Equation 122.. Further simplifications then yield our final expression of the approximate tall posterior score formula, stated in Lemma 3.2 which forms the basis of an efficient and stable implementation. Full proofs are postponed to Appendix A
Lemma 3.1. Let A(0) دز= (7) + (1 ET) and assume it is positive definite. The approzimate log-corvection term is defined using
13 and can be wrillen as a linear combination of Gaussian log-factors
log LA xi.n) = log / (0000)) 1 II go (ool0,, x3)
(14)
(00)) + (1 - n))x(()) - (On)
(15)
where Sk (Hu, kx,k) - mlog 2x log | + u (ek for to E Kj= I...., and Call (1,-1,,, A ) with 17 = Σÿ= En Has + + (1- n ML.
Lemma 3.2. Under the same assumptions as in Lemma 3.1 and using the approzimate log-correction term
from [8] the tall posterior score in 12) can be approzimated as
Vo, 10g (in)) A(,))-
Σ 2 (o.) lo x2)) (1 10, log +F, (16)
where F= 2((,,21..1 = if for all 15J n, V (e) = and V Vo,Ex, = 0.
Formula 16 defines our approximation of the tall posterior score and depends solely on the individual prior and posterior scores (which are known analytically or estimated via NPSE), and F (explicited in the proof in Appendix A.. By enforcing constant covariance matrices Σ,,x and Σr,j, the residual F-term vanishes. This gives us a tractable and practical formula, which reduces to a precision-weighter version of the score formula from F-NPSE, enabling the use of deterministic samplers such as DDIM.
Remark. An alternative deterministic sampler is proposed in Geffner et al. (0223 (Appendix D), which also avoids Langevin dynamics by composing Gaussian reverse transitions 4t-1/e- However, imposes a shared isotropic variance, whereas our approximation preserves the covariance structure of every single posterior and the prior. In particular, our approach exact when individual posteriors and the prior are Gaussian, while theirs is not unless all covariances isotropic. A detailed thcoretical comparison is provided in Appendix M.1 highlighting the advantage of our proposal in faithfully approximating the true reverse dynamics.
3.3 Algorithms
We now turn to the practical side: computing the covariance matrices (0,) in a way that ensures they remain constant. We present two different strategies, which lead to two algorithms implementing our approx- imate tall posterior score from (16): JAC (Algorithm [] and GAUSS (Algorithm 2). For some prior choices (e.g. Gaussian) Σ,, is by construction constant and can be computed analytically (see Appendix D). otherwise, we use the same strategy for Σnj, represented by the prior fn function in both algorithms.
TTe Gaussian log- factors (denated by the seta functions Gj. Cx Cat) represent the (approximate) Gaussian contributions of the backward diffusion each single posterior and the prior, and their global contribution.
7

---

Published in Transactions on Machine Learning Research (02/2026)
Jacobian approximation (JAC, Algorithm 1 Following Boys et al. (2023 it can be shown that (o)( = -ar Va, (o,,4, with (+,( defined via the posterior score as in (6 Et, can therefore be approximaled by taking the Jacobian of the learned score 5,(0,, 23:1) () (in orange). As in (Boys el al. 2023). we do not propagate gradients through Ee,, rendering F(0c, drawbacks. First, need to calculate take the derivative w.r.t. the inputs of the score neural network, which is known to be unstable.
O. But this approach has two main x matrix which is prohibitive for large m. Second, we have to
Gaussian approximation (CAUSS, Algorithm 2 As an alternative, we consider a constant Gaussian approximation the covariance matrix, which automatically gives F(O,, ( = O. Indeed, in the case where p(0,|z)- = N (un, Eo); we have that Σ, (2i²² + Ste Im ) as derived Appendix D Note that choosing =0= Im results in the approximation proposed by SSng et al. (2023). The idea behind GAUSS is to use this formula approximation of the real covariance Σt.. To do so, first estimate Eo,5 for cach zj. by running DDIM with small number of iteralions (>v 100) for each and computing the empirical covariance matrix of the resulting samples (represented in blue).
Algorithm JAC
Algorithm 2 GAUSS
Input: 9, Fisns prior fn
Input: 0,, #im 1, Eim,, prior fn
Output: $1:n
Output: 8im
& 82 & -prior ff,,
o2 82 + -prior fn(Or,()
for j - 1 to n do
for j +
do
Sj f sp (b,, 2))
8j f- s (000,, Σ + E + 2
$-1 + 2i (Im + Ut Ve
[.j
end for
end for
A
(1 - n) Eド= $-1
+ (1 11) n + Eド= D-1
k
$1m
(1 n) ni + 1== 2
1in + (1 - n) 82 22- Σド= Dij
81m
LinSolve A, 81.m)
81:m LinSolve (A, 81:n)
The approximate score obtained via JAC or GAUSS can now directly be plugged into deterministic score- based samplers, such as DDIM, to infer the tall posterior. This should enable faster and more stable inference compared to Langevin sampling. We verify this statement with experimental results in the following section.
4 Experiments
We investigate the performance of JAC and GAUSS compared to F-NPSE, referred to LANGEVIN, on different tasks with increasing difficulty:
(Section 4.1
two Gaussian toy models for which all quantities interest are known analytically
(Section 4.2
2
various SBI benchmark examples, where the score has to be learned via NPSE
(Section
4.3
the Jansen and Rit Neural Mass Model, challenging real-world example
For JAC and GAUSS, infer the tall posterior by plugging the scores obtained via Algorithms andd into the DDIM sampler defined in Section 2.1 For F-NPSE, we use the unadjusted Langevin algorithm (ULA) from (Geffner et al. 2023) with I= 5 Langevin steps per time sime and 7 = 0.5. In all cases, we use a uniform time schedule (4: = T- We evaluate the quality of the inferred tall posteriors using distance- based metrics, comparing them, when possible, to samples from the true tall posterior. Each subsection details task-specific evaluation setups. Further implementation details can be found in Appendix reproducing all experiments is available at http:////gituuuu.ooo///aaaaaae http:///wooo.ooooooooaoo//oooonn wttpt//ffffffosoofooto
D and the code
Remark. We additionally compare with the deterministic sampler from Geffner et al. 2023 mentioned
in the Remark from Section 3.2 Full empirical results are provided in Appendix M.2 and complement the
thcoretical comparison from Appendix M.1
8

---

Published in Transactions on Machine Learning Research (02/2026)
4.1 Gaussian toy models
We consider two toy examples for which the analylic form of the posterior and corresponding score functions are known: a multivariate Gaussian p(z 0) [N (1 - + plm) with correlation factor 0.8, and a Gaussian Mixture Model (GMM) p(z () = 0.5 N (z; 0 2.255) + 0.5 N (x; 0. Σ/9), where E a diagonal matrix with values increasing linearly from 0.6 to 1.4, as in Geffner et al. (2023). Both examples are carried out with Gaussian prior (()) = (0; O, Im) with known score. See Appendix for all analytical formulas.
The following experiments evaluate the speed and robustness to noise of each tall posterior sampling algo- rithm, for increasing number of observations We quantify accuracy using the sliced Wasserstein (sW) distance between estimated and true tall posterior samples, over 5 random seeds. Each seed corresponds to a different parameter 0* ~ X(0), used to simulate the conditioning observations Fi. for the tall posterior p(0 zi.n). The true tall posterior is known analytically for Gaussian and sampled via Metropolis Adjusted Langevin (MALA) (Roberts & Tweedie,
[1, 100] and parameter dimensions m E {2, 4, 8, 10, 16, 32).
1996) for GMM.
We model the noise by considering the score estimator 20(00, Ve, logl((0()z) + Evei ro(O,, x an) where Vo, log pa (0.(zz) is the known posterior score and re randomly initialized neural net with outputs in range [-1, 1]. This construction leads to a controlled error of E for the noise predictor (0,, x,tt/ ver, which is what we actually optimize when training a score model (sce Appendix Eb. The input 6g denotes the forward-diifused parameter, 02 = v 0 + vI - E with 0 ~ x(0) and E ~ (O,
Table [ displays the total running time and sW for each algorithm on the Gaussian example. For the same number of time steps T, our algorithm yields smaller sW than the Langevin sampler while accounting for L = 5 times less neural work evaluations (one per Langevin steps). Table 4 in Appendix G shows similar results for GMM. The average speed-up of our algorithms over all considered number of time steps T' and different choices of e and n are shown in Table 2 It shows approximately the same values for m = 2, 8, 10, 16, 32, meaning that the speed-up is not impacted by the data dimension. GAUSS is the most efficient algorithm, being 2.5 times faster than the Langevin sampler. JAC is approximately 1.8 times faster.
Runtimes.
Based on Table D we consider an equivalent time setting with 400 and 1000 steps for JAC/ LANGEVIN and GAUSS respectively. This setting will be used in all subsequent experiments.
Algorithm
steps
At (s)
m
Speed up GAUSS
Speed up JAC
GAUSS
50
0.45 +/- 0.00
0.17 +/- 0.08
2
0.39 + 0.01
0.56 + 0.00
JAC
50
0.41 +/- +.00.00
3.14 +/-4.12 4.12
4
0.39 = 0.01
0.55 = 0.00
LANGEVIN
50
0.83 +/-0.00 0.00
nan +/- nan
8
0.39 = 0.01
0.56 = 0.00
GAUSS
150
0.90 +/- 0.00
0.17 0.06
10
0.38 1 0.01
0.52 1 0.00
JAC
150
1.22 +/- +/-0.00
1.57 +/3
16
0.38 1 0.01
0.54 1 0.00
LANGEVIN
150
2.50 0.01
0.65 0.42
32
0.37 1 0.01
0.52 1 0.00
GAUSS
400
2.04 +/-0.00 0.00
0.20 +/1 0.11
JAC
400
3.26 0.00
0.85 4/2 1.20
Table 2: Ratio between the runtime for GAUSS and JAC w.r.t. LANGEVIN for the Gaussian example in different dimensions m. Averaged over the number of time steps T E {50, 150, 400, 1000), different noise levels E {0, 3 10 1.1013 and the number of observations E (1, 100). Mean and std over 5 different seeds.
LANGEVIN
400
6.65 +/- +/00.02
0.65 +/40.33 0.43
GAUSS
1000
4.77 0.01
0.22 0.10
JAC
1000
8.18 +/- 0.03
0.00.0
0.25 +/-
LANGEVIN
1000
16.65 +/- 0.03
0.65 +/- 0.42
Table 1: Sliced Wasserstein (sW) and runtime At for T time steps for the Gaussian example with m = 10, n = 32 and = 10-2 Mean and sid over 5 different. seeds.
Robustness to noise. Figure 2 portrays the effect of the perturbation in the posterior approximation for cach algorithm, across n E [1, 100 observations. In the Gaussian case, GAUSS performs best in all settings as expected, since the second-order approxinations in our method are exact in this case. For GMM however, the approximation from Section 3.2 is not exact, ie. the backward kernel from [I]] is not Gaussian. We use this to analyze the effect of our approximation while controlling the score estimation error. GAUSS remains competitive with LANGEVIN, outperforming it for small n and matching it for n > 90. In both models, JAC is extremely accurate in the noise-free setting (e 0), but quickly becomes unstable noise increases.
9

---

Published in Transactions on Machine Learning Research (02/2020)
e= 0.0
< 0.001
c= 0.01
c 0.1
1.00
0.75
Gaussian
SW
0.50
0.25
GAUSS (1889)
JAC (489)
LANGEVIN (499)
0.00
1.00
0.75
GMM sW
0.50
0.25
0.00
10
10a
101
102
101
n
n
n
Figure 2: Sliced Wasserstein (sW) distance as function of n and for increasing noise levels Results are shown for both Gaussian loy examples with m 10. Mean and std over different seeds.
We include Figures and in Appendix complete analysis for the Gaussian example, across all dimensions m. They confirm the results of Figure 2 highlighting the precision of JAC in the non-perturbed and its instability heeeeeiee. They also show the superior robustness of GAUSS in high dimensions.
Overall, the above experiments show that our proposal outperforms the Langevin sampler from F-NPSE in terms of speed and robustness to noise. Furthermore, they suggest that GAUSS offers good trade-off between precision and robustness, hence a better algorithm choice than JAC in SBI settings, where the posterior score is unknown and has to learned. We investigate the validity of this statement the following section.
4.2 Benchmark SBI examples
We consider three examples from the popular SBI benchmark (Lueckmann el. al. 2021 with different. and z-space dimensions (resp. m and d), and different posterior structures:
SLCP (m 5, = 8): Uniform prior and Gaussian simulator, whose mean and covariance non-linear functions of the inpul parameters 0. Multi-modal posterior.
SIR (m = 2, d 10): Log-Normal prior and simulator based on a set of differential equations that outputs samples from a Binomial distribution. Uni-modal posterior.
Lotka-Volterra (m d 20): Log-Normal prior and simulator based on set of differential equations that outputs samples from a Log-Normal distribution. Uni-modal posterior.
The score corresponding to each prior is analytically computable, as done in Appendix D However, contrarily to the toy models from Section 4.1 the analytical posterior score is not available and has to be learned via score-matching. The goal is to evaluate the robustness of our sampling algorithms in this learning setting, We train simple MLP on Nerain samples over 5 5000 epochs using the Adam optimizer (sce Appendix E According to the equivalend time selting from Section 4.1 we use T 1000 steps for GAUSS and T 400 steps for JAC and LANGEVIN. We also consider clipped versions, ensuring that the samples every step stay within the high probability region of standard Gaussian (truncated to [-3, 3)). The goal is to stabilize JAC and LANGEVIN, that were found less robust than GAUSS. This may however slow sampling and introduce bias.
To assess the performance each sampling algorithm in this learning setting, we evaluate their accuracy for increasing n E [1, 8, 14, 22, 30) under varying Nurain {103, 3. 3.103. 10 30113.11") Our empirical evaluation sam- ples 25 ground-truth parameters 8* X(0) from the prior, used simulate the conditionning observations 2 p(æ 0*) for 1,..., for the tall posterior p(0 t. For all three tasks, reference samples from
10

---

Published in Transactions on Machine Learning Research (02/2026)
1
8
14
11 22
11 30
2.0
GAUp
1.5
JAC-clip LANGEVIN/-CLiD
SLCP
sW
1.0
0.5
0.0
Lotka-Von Volterra
0.04
sW
0.03
0.02
0.01
0.008
sW SIR
0.006
0.004
0.002
0.000
103
104
10a
10
10a
10a
10a
10
10a
10a
Nizaio
Nirain
Nirain
Nusin
Nizain
Figure Sliced Wasscrstein (sW) distance as a function of Nirain and for increasing n, between samples obtained by GAUSS, JAC and LANGEVIN, and the true tall posterior p(0 wi.)). Mean and std over 25 seeds.
the true tall posterior are obtained via MCMC using numpyro (Phan et al. 2019). and used to compute bissence-based evaluation metries. Outliers with values above the 99th percentile (or NaN) are excluded.
Figure 3 portrays the sliced Wasserstein (sW) distance for each task function of the size Nirain of the training set for the score model. Larger values of Nerain are expected to yield better score estimates and thus more accurate posterior approximations. This is represented by decreasing tendency of the sW. Overall, we observe that GAUSS outperforms all other algorithms. It yields consistently lower distance values and and
scales to high n values, compared to LANGEVIN/-clip that diverges for n 2 14 for the SIR examples. Note that JAC diverges as soon as n. >> 1, which why we didn't include it the plots. On the other hand, JAC-clip is more less equivalent to GAUSS, with slightly better results for SLCP.
We complement these results with two additional metrics: Maximum Mean Discrepancy (MMD) and Clas- sifier Two-Sample Test. (C2ST), both from the sbibm package. Their relevance is motivated by the differ- ent posterior structures of the benchmark examples (sce Figure 10 in Appendix SLCP is multimodal, Lotka-Volterra is high-dimensional but unimodal, and SIR is low-dimensional and smooth. Each metric captures different aspects the posterior approximation, making them complementary tools for assessing statistical validity. This summarized Table B full numerical results are in Appendix H They highlight that, across all metrics and all tasks, our methods significantly outperforms the Langevin sampler.
Metric
What it measures Geometric alignment of distribu- tions based on sliced ID projections
Strengths and limitations
Best method for n > GAUSS (all tasks), JAC-clip (SLCP)
sW
Sensitive to multi-modality (SLCP); less et- fective in high-dimensions (LV) and in the absence of geometric structure (SIR)
MMD
comment-based global differences using kernel methods
Effective in detecting smooth deviations and higher order correlations; known to fail multi-modal settings (SLCP)
GAUSS (SLCP, LV), JAC-clip (SLCP, SIR)
Classification accuracy between C2ST
Fast and general validation tool, capturing any inconsistency; can be "too discrimina- tive" (LV, SLCP) and is harder to interpret
GAUSS (LV),
true and approximate samples
JAC-clip (SLCP, SIR)
Table 3: Metric description and results summary. Each metric provides complementary information on the statistical validity of the approximate posterior and more less appropriate for the considered benchmark task. The last column displays which method performs best in the tall data setting (n 1).
11

---

Published in Transactions on Machine Learning Research (02/2026)
The results of this section show that our methods outperform the LANGEVIN baseline, particularly in scaling to large observation contexts, even in a setting where the posterior score is learned from data. While GAUSS is overall the most accurate, the strong performance of JAC-clip suggests it is a viable alternative when stability can be ensured. However, a notable trend across tasks is the performance degradation as n increases- as highlighted by in figures of Appendix H. Indeed, the compositional approach introduces a key limitation: approximation errors accumulate as we sum over n cvaluations of the learned score model to obtain the tall posterior score. We explore a potential solution to this issue in Appendix L.2 via partially factorized methods, offering trade-off between simulation cost and number of network evaluations.
Finally, note that we chose to limit our comparison to the F-NPSE bascline, as Geffner et al. 2023) already demonstrated its strong performance relative to other SBI methods. This supports our motivation to improve upon F-NPSE directly, and our results show that the proposed approach is indeed promising. We include results for additional benchmark tasks in Appendix [] confirming the trends observed this section.
4.3 Inverting a non-linear model from computational neuroscience
In this Section, we illustrate the benefit of the tall data setting for challenging real-world applications and consider the Bayesian inversion of the Jansen and Rit Neural Mass Model (JRNMM) (Jansen & Rit 1995). The output x(t) of this simulator is a time series obtained by taking as input a set of four parameters 0 = (C, μu, o,g). Appendix J.1 gives a full description of the JRNMM. Importantly, there exists a coupling- effect of parameters g and (H, on the amplitude of the output signal æ(t), meaning that the same observed signal could be generated for different pairs of and (H, This indeterminacy makes the inference problem ill-posed (Rodrigues et al. 2021]. We now demonstrate how a tall data posterior can give more precise information on how to invert the simulator, considering two cases:
3D JRNMM: a simplified setting in which we fix g 0 to lift the indeterminacy on (W, The Iall posterior is expected to concentrate around the true parameter values.
4D JRNMM:
the full JRNMM, with coupled parameters g and (H,o). The tall posterior should reveal the resulting indeterminacy on (J,
In each case, we estimate the tall data posterior of the JRNMM using GAUSS, JAC and LANGEVIN, according to the equivalent time setting from Section 4.1 The posterior score was estimated via NPSE using a MLP trained on Nurain 50 000 samples from the joint distribution, with a uniform prior placed over the range of physiologically meaningful values of the simulator parameters. See Appendices
E and J.2 for further details.
We first evaluate the accuracy of each tall posterior sampling algorithm. This time, samples from the true posterior are available and previously considered metrics cannot be computed. We therefore use the local Classifier- Thassifeer-Tw-SSampll Test. (4-C2ST) {Linhart et al. 2023) and the default implementation from the sbi Python package (Alvaro Tejero- Cantero et al. 2020). C-C2ST compares conditional distributions by training a classifier on scparate calibration set from the joint distribution. The validation results are given and explained in detail in Appendix J.3] They show that GAUSS is the only method that passes the test in both, the 3D and 4D cases, confirming the reliability of GAUSS in this challenging setting and motivating its use in the remainder of this section.
We now investigate the ability of our tall posterior to concentrate around the true parameters 6* used to simulate the observations zi.n as n increases. This is quantified with the MMD between the postcrior marginals and the Dirac distribution dg- at the true parameters 0*, computed on 10 samples obtained with GAUSS. Figure 4 shows the results for the simplified 3D JRNMM case. We observe as expected a progressive concentration around 0*, with sharper posterior marginals and decreasing MMD. Figure the results for the full 4D JRNMM. Here, the tall postcrior takes 'sharpened banana' shape: we observe a progressive convergence of the inferred posterior mean lowards 0* (black dots and dashed lines), sharpening around (C, 9) and a sustained dispersion along the dimensions of (H, which explains the non-decreasing MMD, This gives evidence about the indeterminacy caused by the coupling between g and (n, a), which is consistent with the results from [Rodrigues et al. 2021
12

---

Published in Transactions on Machine Learning Research (02/2026)
MMD Dirac
60
GAUp
==
и= 14
40
n=
==
1
8
14 22
30
130
140
200
300
400
500
1000
2000
3000
C
H
6
Figure 4: Inference on the 3D JRNMM (fixed g = 0) with GAUSS. (Left): MMD between the marginals of the approximate posterior and the Dirac of the true parameters 0" (black dashed lines). (Right)): Histograms of the 1D marginals of the inferred posterior for 30 single observations (n 1) and sets 21., of increasing size.
n1
n=
n= 14
n=
400
30
200
MMD to Dirac
GAUSS/-clip
400
3500
200
1000
0
1
14
22
30
20
71
50
-20
100
250
200
400
1000
3500
-20
20
#
6
g
Figure 5: Inference on the 4D JRNMM with GAUSS. (Left): MMD between the marginals of the approximate posterior and the Dirac of the true parameters 0* (black dots and dashod lines). (Right): Histograms of the ID and 2D marginals of the inferred posterior for observation sets rin of increasing size.
5 Conclusion
We presented new sampling approach for tall posterior inference in iiatulaaioo-bbsee settings, eliminating the need for Langevin dynamics as used in F-NPSE (Geffner et al. 2023), while retaining its amortized and compositional benefits. By explicitly approximating the diffusion process of the factorized tall posterior, our method enables faster and more stable inference through established score-based samplers such as DDIM.
Across toy models and benchmark tasks, our method consistently outperforms the Langevin bascline in terms of speed, accuracy, and scalability with increasing observation set sizes. In particular, GAUSS emerges as a robust and consistent choice, while JAC achieves highly accurate results when score estimates reliable.
13

---

Published in Transactions on Machine Learning Research (02/2026)
These findings validate our central motivation: improving compositional inference by modeling the diffusion dynamics of the tall posterior directly. We further demonstrated that GAUSS scales effectively to real-world problems, as illustrated by our neuroscience application. This example underscores a key benefit of the tall data setting: aggregating observations can uncover parameter dependencics that are hidden at the single- observation level. This insight aligns with recent hierarchical extensions of amortized inference (e.g. HNPE), suggesting that structured modeling across observations can improve parameter identifiability.
Looking ahead, several research directions emerge. JAC has strong potential, but requires better stabilization. Clipping samples improves robustness at the cost of slower, biased inference. Recent work by Gloeckler et al. (2025) suggests promising alternatives. Another interesting direction is the refinement of our socond-order approximation. For example, one could try to combine our fixed-covariance strategy with data-driven covari- ance estimation methods like those from (Rissanen et al. 2025). to better approximate local curvature with- out sacrificing stability. Finally, while partially factorized methods help mitigate error accumulation -he main limitation of compositional inference methods-furtiee complementary strategics are needed make our method more scalable, beyond hundreds, thousands of context observations. Arruda al. (2025) suggest one way to do that.
Together, our results demonstrate the power of score-based methods for amortized and compositional infer- ence. Beyond tall posteriors, we believe this framework opens new possibilities for modular, hierarchical, and other structured inference tasks across scientific domains.
14

---

Published in Transactions on Machine Learning Research (02/2026)
References
Markus Ableidinger, Evelyn Buckwar, and Harald Hinterleitner. A stochastic version of the Jansen and Rit neural mass model: Analysis and numeries. The Journal of Mathematical Newroscience, 7(1), August 2017. doi:
Jonas Arruda, Vikas Pandey, Catherine Sherry, Margarida Barroso, Xavier Intes, Jan Hasenauer, and Ste- fan T. Radev. Compositional amortized inference for large-scale hierarchical bayesian models, 2025. URL
attps://arxiss http://rwwvv.orraa////aa//// wtt.oor//////////...o 14429
Rémi Bardenet, Arnaud Doucet, and Chris Holmes. On Markov chain Monte Carlo methods for tall data. Journal of Machine Learning Research, 18:1-43, 131143,44 5 2015. ISSN 15337928.
Georgios Batzolis, Jan Stanezuk, Marola-Bibiane Schönlieb, and Christian Etmann. Conditional image generation with score-based diffusion models, 2021.
Rajendra Bhatia. Positive definite matrices. Princcton Scries in Applied Mathematics. Princcton University Press, Princeton, NJ, December 2006.
Benjamin Boys, Mark Girolami, Jakiw Pidstrigach, Sebastian Reich, Alan Mosca, and O Deniz Akyildiz, Tweedie moment projected diffusions for inverse problems. arXiv preprint orXi::2310.0677 2023.
Nicolas Brosse, Alain Durmus, Eric Moulines, and Sotirios Sabanis. The tamed unadjusted langevin algo- rithm. 2017.
Evelyn Buckwar, Massimiliano Tamborrino, and Irene Tubikanec. Spectral density-based and measure- preserving ABC for partially observed diffusion processes. an illustration on hamiltonian SDEs. Statistics and Computing, 30(3):627 648, November 2019. doi: 10000//11222201--999999
Kylc Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences (PNAS), 117:30055 30062, 2020. ISSN 0027-8424. doi: 10.1073/ pnas. 1912789117.
Maximilian Dax, Stephen R. Green, Jonathan Gair, Michael Pürrer, Jonas Wildberger, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf. Neural importance sampling for rapid and reliable gravitational-waay wave inference. Phys. Rec. Lett., 130:171403, Apr 2023. doi: 10.1103/ 10.1103/PhysRer// 1hysRerLeett.13. 171403.
Tomas Geffner, George Papamakarios, and Andriy Mnih. Compositional Score Modeling for Simulation- Based Inference. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th International Conference Machine Learn- ing, volume 202 of Proceedings of Machine Learning Research, pp. 11098-11116 PMI.R, 23-29 Jul 2023.
Manuel Gloeckler, Michael Deistler, Christian Weilbach, Frank Wood, and Jakob H Macke. All-in-one simulation-basie inference. 2024.
Manuel Gloeckler, Shoji Toyota, Kenji Fukumizu, and Jakob H. Macke. Compositional simulation-based inference for time series, 2025. URL https:://axiv.oog/aabbssse wwtpw/////.aa......//./// 02728
Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad, Tim P Vogels, David S Greenberg, and Jakob H Macke. Training deep neural density estimators to identify mechanistic models of neural dynamics. eLife, 9ce56261, sep 2020. ISSN 2050-084X. doi: 10.7554 / LeLiee5666..
Ian Goodfellow, Jean Pouget-Abade Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshma Bongio. Generative Adversarial Nets. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger (eds.), Advonces in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
15

---

Published in Transactions on Machine Learning Research (02/2026)
David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. In Kamalika Chaudhuri and Ruslan Salakhut (eds.), Proceedinggs of the 36th International Conference on Machine Learning, volume 97, pp. 2404 2414. PMLR, 09 15 Jun 2019.
W. K. Hastings, Monte carlo sampling methods using markov chains and their applications. Biometrika, 57: 97-109, 1970. ISSN 0006-3444. doi: 10.1093/ /10./III/E///......
Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-free MCMC with amortized approximate ratio estimators. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 4239 4248. PMLR, 13-18 Jul 2020. doi: 00448550/arxxxx1103...000
Jonathan Ho and Tim Salimans. Classifier-free diffusion guidance. In NewrIPS 2021 Workshop on Deep Generalive Models and Downstream Applications, 2021. URL https://openreviior httpp//eerevvi...eeeefforrrrrooo net/forummidd
qW8AKxfYbI
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 6840-688 6801, 2020.
Ben H. Jansen and Vincent G. Rit. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics 1995 73:4, 73:357-366, 1995. ISSN 1432-0770. doi: 10.1007//BF00199477.
Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 26565-26577, 2022.
Julia Linhart, Alexandre Gramfort, and Pedro L. C. Rodrigues. L-c2st: Local diagnostics for posterior approximations in simulation-based inference. 2023. URL https://arxiv. hwtp:////g/ab/aaba/eeeee 03580
Jan Mathhss Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking simulation- based inference. In Arindam Banerjee and Kenji Fukumizu (eds.), Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pp. 343-351. PMLR, 13-15 Apr 2021.
George Papamakarios, David Sterratt, and Iain Murray. Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In Kamalika Chaudhuri and Masashi Sugiyama (eds.), Proceedings of the Twenty-Sean Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pp. 837-848. PMLR, 16-18 Apr 2019.
George Papamakarios, Eric Nalisnick, Danilo Jimencz Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22:1-64 2021. ISSN 15337928, doi: 10.48550/amda arxiv. 191...9991..0....
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjic Bai, and Soumith Chintala. Pytorch: An imperative style, high-perforrance deep learning library. CoRR, abs/ Ab1/1192.011703, 2019. URL http://aaaii...rr//a//9991...1. 01703
Du Phan, Neeraj Pradhan, and Martin Jankowiak. Composable effects for flexible and accelerated proba- bilistic programming in NumPyro. arXix:1912. 11554, 2019.
Python Software Fundation. Python Language Reference, version 3.6, 2017.
Severi Rissanen, Markus Heinonen, and Arno Solin. Free hunch: Denoiser covariance estimation for diffusion models without extra costs, 2025. URT. hhttps://arx:: htt///aawww..ora////////0///
11149
16

---

Published in Transactions on Machine Learning Research (02/2026)
Gareth O Roberts and Richard L Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pp- 341-363 1996.
Pedro L. C. Rodrigues and Alexandre Gramfort. Learning summary features of time series for likelihood
free inference, 2020. URL
Pedro Luiz Coelho Rodrigues, Thomas Moreau, Gilles Louppe, and Alexandre Gramfort. HNPF: Leveraging Global Parameters for Neural Posterior Estimation. In Neural Information Processing Systems, 2021.
Louis Sharrock, Jack Simons, Song Liu, and Mark Beaumont. Sequential neural score estimation: Likelihood- free inference with conditional score based diffusion models. 2022.
Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics, 2015. URL https://arxiv..rg/ htt///bbs/1/0/ 03585
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In International Conference on Leurning Representations, 2021a.
Jiaming Song, Arash Vahdat, Morteza. Mardani, and Jan Kautz. Pseudoinverse-guided diffusion models for inverse problems.. In International Conference on Learning Representations, 2023.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Н. Wallach, H. Larochelle, A. Beygelzimer, F. d' B'Ac E. Fox, and R. Garnett (eds.), Advances in Newrul Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Yang Soug, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based gencrative modeling through stochastic differential equations. In International Conference Learning Representations, 2021b.
Pascal Vincent. A connection between score matching and denoising auloencoders. Neural compulation, 23 (7):1661 1674, 2011.
Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: comprehensive survey of methods and applications. ACM Computing Surveys, 56(4):1 39, 2023.
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexan- der J Smola. Deep sets. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Álvaro Tejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueekmann, Conor Durkan, Pedro J Gonçalves, David S Greenberg, Jakob H Macke, Computational Neurocngineering, and T U Munich. Sbi loolkit for simulation-basee inference. Journal of Open Source Software, 5:2505, 7 2020. ISSN 2475-9066. doi: 10.21105/joss.02551 02505.
17

---

Published in Transactions on Machine Learning Research (02/2026)
A Proofs
A.1 Proof of Lemma 3.1
Lemma 3.1
and assume it is positive definite. The approzimate log-correction term defined using 19 and can written as a linear combination of Gaussian log- tog-factors:
Let A() - (
0 (1-
log Ix xi.)) = log (010,) II
a7)00
(17)
Σ (y(0,) + +(1 11- n)) - San(
(18)
where Sk = ((1,, Tx,x) = - mm 2x - log 1r + N1 (i for to E {j 1,..., n; )} and Sall C(A ((A-In,A~)) with 17- - 2 Muj + (1 n ,,,
Proof. Let (UR)ISKSK (Rmnk and (Ex) be covariance matrices in R=xm, Denote by pk the Gaussian pdf with mean μe and covariance matrix Σg. Note that
Ph 00 + exp (a + (Ex ) TT)
1 00 Ex 00
2
where
&o = - (m log 2x log 2 + ur 2ixe pmx
Therefore,
K II (00)
Σ 54
Eall
"ull + 17 -TT00 0u
1
A0o
exp
exp
2
k=1
K
Σ
Čo -all
.o(0;; 11,7,A
= exp
with ij = Σ j pk, Ā Σx=1 Dx and
Eall - (mlog 2 - log A + n "/>A'1n)/2
We can apply this
result to Equation
[7] with
(0))
Σ (0)
j=1
A(0) = Σ
E + (1- +11m))=71(0)
2=1
since by assumption A(0) is definite positive. This provides the following reformulation of Equation 17
log L1 (f,, = log f exp Σ ((0)) (1-m)(((0)) - (l,( N(66: ((0e:A)
Finally, because the sum in the above formula docs not depend on the integrator variable 9o; we can take it out of the integral and separate the expression of L,x(0, ai.)) into two terms:
log (A(00, = log NV(06: A (0))0)0)000)) + Σ (0.) + + (1- - (arr(ft) (19)
The first term is zero the integral of the Gaussian p.d.f. over the whole space is equal to 1. This leaves with the second term and concludes the proof.
18

---

Published in Transactions on Machine Learning Research (02/2026)
Positive Definiteness of A()) The approximation from Lemma 3.1 is only valid if the matrix A(0) is symmetric positive definite (SPD), i.e. is valid covariance matrix. In what follows, we investigate what conditions on the prior distribution will ensure such property. Using the partial ordering defined by the convex cone of SPD matrices BBhtia 2006), it follows that:
A(0) 0
21 + +1- (1-))
(20)
==1
R
1-1(0) t.
(0) (n - ()=1)=
(21)
==1
1 (n 1 2 (0)
n Σ
Eax(8)> >
(22)
Note that n increases, the right-hand side of the inequality in Equation 22 converges the harmonic mean of the (0), which helps building intuition for the correct choice of (0), the covariance of the backward diffusion kernel associated to the prior distribution X(0). For instance, a sufficient choice for A(0) to be SPD, would be to have a (0) whose associated ellipsoid covers the ellipsoids generated by the covariance matrices (O), associated to the posteriors p(0 (z)) for every single observation x considered in the tall data inference task. Intuitively, this corresponds to choosing a prior that is broader than the posterior distribution, which is normally the case Bayesian approaches.
A.2 Proof of Lemma 3.2
Lemma 3.2 Under the same assumptions as in Lemma 3.1 and using the approzimate log-correction term from 18 the tall posterior score in 12) can be approzimated as
Va, (0: xi.n) ~ A(0.)-1
(1-n)2, t. (0c) Ve, log (z)) + (1 - n) Vo. logp)(()) + F.
where F((0e, xim) = is Vo, (n,( = for all 1 isisw n and Vo, Ex.(O)))
Proof. Remember that the full tall posterior score writes
Ve, (0 ) = (1 - Vo. (0)) + Σ Ve. log (0, 23) + Vo, log (x(0,, xin))
(23)
1=1
We can replace the log-correction term log LA by its estimator log Lx from 17 in Lemma 3.1 compute that term explicitly. Taking the gradient directly gives us:
Now let's
Ve, log (10, = Σ (oL)((1) (1- (o. Ca.
(24)
where we recall that --)=)-=-((( 2x log 1ー + 1 /2.
6Te ellipsoid EA associated with SPD matrix A defined as EA (x xT x 1).
19

---

Published in Transactions on Machine Learning Research (02/2026)
First, compute the terms (0). The chain rule gives us
(o) Vos(uu, (0),
= Vus(mu)((), (0)(),,(00)_T](0)) T (0) + 1=-((4,(10),,,,))000 (0),
where
VE'C(H,-1) s -1111
Therefore, we obtain
1 (E) (p.,(0)=,))(() Vol,,((0) (0). Volz(0) = (1, (1))),)))0 Veper + 2
Note that for (0.(6))) = N( va,00. mm 2023 which leads to
we have that (o.Hi.(01) (var/ (0, [ooss et al.
Vau (0.))
Me, ( + 1 (0t) - Me, (e, (o)) (0:) Vo, Dī
Vi
(25)
1 -- 1 Va, 10g I mz)) + (Er)((1) Mc. (4)) Vo, (.))
Formula 25 also applies to (2.0) (with pi, Me,,, Eea)- Replacing these expressions in 24 yields
Va log (0r, (i.,) -(1 (1- - n) Vo, logr (0.)
Va (0 | xz)) vz 0. Vo, Con(
n
[Et, pt.j (0)((0)))e_0000))) Vo, [ (0.)
j=1
n
(E))(01) M1, (0)) Vo. (0,,
2
Replacing the above formula in 23 the first two terms compensate each other, which gives the following approximation for the score:
Vo, logp, (0: wi.n)
- -u 02
(O.)
n
+
(Er)(Ot) - (M.,(0))) Vo.
(26)
2
j=1
1
(Er,)(0)) (.x(0:))xxx00²))) Vo, (0+
2
We compute VoCau((). By noting that Call before and obtain
can follow the same steps
(oCan( =
+ 2 II A(0))~1(000)(0) -V0))
(27)
1 Im A(0)='n(0)m(0)) A(0) - 'VA((0).
= - ( VA(0)-
((0) +
Varraane preserving (VP) framework introduced Section
20

---

Published in Transactions on Machine Learning Research (02/2026)
Note now that
Von(0) Σ
+ 10t -n) (Zixt()Vo)):1(0) + m(a(0) VoL
at Im
M1 (10 (oI))(().
Uz
j=1
and that
a vai
Σ
2 (0) v (111)) - (1==)),;(() (0 ve ve (0))
==1
(28)
= A(000 +
( 2z) + (1 n) log
We can use the above formulas to replace (en(0))(() in the term A(0) of Equation (27 The first term of this new expression - 0, which (setting = 0)) compensates with the first term in the score formula from 26 and finally leaves us with the wanted score approximation:
Vo, log xin)) A(6.)-1
E (o.) (0, a0)) + Vo, logp)(()) F(0r,
where
F(0, x (0) (V(()) -n()) 1
2 [Im A(0) (()))= A(0) -VA()))))())
+ (E) pr (0) (0) (
+
- (1)) (,5,))))))
F(O, æi.n) contains all the terms that depend on the gradients the covariance matrices Ee, (0) and Σx,
and is 0 if they are considered constant (soe definition (()) Lemma 3.1
21

---

Published in Transactions on Machine Learning Research (02/2026)
B Influence of the correction term
In this Appendix we provide some intuition on the influence of the correction term Vo, log L Lx(0,, Finm) in the backward diffusion process. Following the definition of the backward kernels in equation (11
For t 1: the forward kernel and the diffused data distribution approach the noise distribution: pe(O, 2) : 0, and (eo((e 0o) = N/ 0, Im). The backward kernel is thus equivalent to the target data distribution:
Pojo(8o x2, p(00 x) x)) (0. 0o)
1-11 p(Oo x))
(29)
(0 x))
Therefore the backward kernels vary very little with 8, and because they define log (0,,, (see eq- (17)), its gradient is close to zero. In other words, Vo, log (x, rim) has no significant impact at the beginning of the backward diffusion (aka. sampling or generative process).
For + 0: the denominator is the diffused distribution that gets close to the target data distribution Pe(6: x² 3o p(0b z). Therefore the backward kernel is approximately
(o Oz, a) p(0) x (2) (0: (o)
uoj((0, 0o) : 60, (0.)
(30)
Pr(O, | x))
Here, the dependence on convergence a Dirac function. This means that the gradient of log Lx(0,, will increase during the sampling process and finally explode when t approaches O. The correction term therefore plays an important role as we approach the target tall data posterior distribution, at the end of the sampling process.
22

---

Published in Transactions on Machine Learning Research (02/2026)
C Denoising Diffusion Implicit Models (DDIM)
In this section we derive the DDIM backward Markov chain mentioned in Section 2
DDIM rntroduces a set of inference distributions defined for t E [1: TT as
11111,,00 (0: 02+1) = (0.-1; Oz) of
(31)
-e( Note that f 2021a la Lemma 1, Appendix B). is chosen that (0 (06) 2: (0: 00) = N0ci va,lo. vmlm) (Song et al.
with pe (00; 0z)
10 (vz - ( /oret))/ /v/ //2² and for E (0, 1/
1/2
This property allows us to write (0--() = [I of 11-111,000-111 0c. 80 2oje (0 8, pp 0 d00d0,
Even though Equation 31
suggests a way of passing from Pt to Pt-1, it involves an intractable kernel J ( lr. 2ojo (o dd0g. DDIM approximates the marginal distribution with
(0--)) -f f 61-111.000-1104. (e
(32)
where (e()))= Ego [0)], and verifies (0,) u/ Be m score from Equation 3 we can now define the following approximation of μr (0,):
V log 9njo (0, Using the learned
(0,)
var ( + (rs,
(33)
We finally obtain the following backward Markov chain for DDIM:
Po.o.((Uu:)) = PT((O)) II ga. (1((0/-10)))
(34)
where (0r) pr(Or), (0,--11(01-110))
L (0 (o,,(()) and
40.011 (8)18))
(0)) m)), with ag > 0 free parameter.
23

---

Published in Transactions on Machine Learning Research (02/2026)
D Analytical formulas for score and related quantities
D.1 Gaussian case
The considered Bayesian Inference task is to estimate the mean E I of Gaussian simulator model p(r 0) (æ; 0, )), given a Gaussian prior A(0) = (0; μ^, Ex). For a single observation x", the true posterior also a Gaussian obtained using Bayes formula, the product of two Gaussian distributions:
p(0) x²) = ((;; Mpost (x*); Expstt
(35)
with poost (post(x*)) Epost 1--1x* +Ex MA) and Epost (E ---1+ ) can be written as
Note that in this case, the full posterior
p(0 xiin) = N0: Mpost (z(in); Epost,
(36)
with Epsst,n = (nE-1 + + Dー and Moost (ii.n)) = Epost,n 2j=1 Σ- + +E px).
Assume that (0 (0) = NN(0; va,to, 2yIm)); see Section 2 Using standard results (e.g. Equation 2.115 in Bishop, 2006), we can derive the analytic formula of the diffused prior
(0) / ((00)))0000 0o)d0o
(37)
(0,5 yauua, oz Ex + vm
and of the diffused posterior
f p(00 x² x=)000((0: 0o)d0 (0a; Vostppott (x*), a poott + veTm).
(0; |x*)=
(38)
The corresponding Fisher scores
Va, logr; (0.) (ar22s + vzlm) -0: v(rp)))
(39)
Vo, (0, x²) = (pa22ipont + vilmm) (0: Vau Post (x"))
(40)
Replacing the score from 40 in the formulas from (Boys et al. the mean and covariance matrix of the backward kernel (oje(00 0,,x*) x
2023), we get the following expressions for
x(), = 2t 1 2s Epost + vrIm
)-1)
= E,, (constant)
x+) = Vat
1 u1 Epest + +etmm)
Oe + ve aet + v((m) Mpost (x*)
Vot E,Ot + (prEp + vetm Mpoot (x*).
=
The same goes for the prior backward diffusion kernel 16 0):
Vt
-1 )
(o,)
1 - v2
2x + vrlm
= Σx,
(constant)
Or
1
-1
-1)
(0.)
1
(a + vm
0r + Us (atEs + vtIm) 1x
Vau
Vai
En +
(a+D>> + vzImm
p1.
21
24

---

Published in Transactions on Machine Learning Research (02/2026)
D.2 Mixture of Gaussians case
In the case of the Mixture of Gaussians, the prior is À = (0, Im), the simulator is given by p(x10) [x 0, 21) + [2 0, 1/922) where Σ1 = 2.252, De = 1/9 and Σ is diagonal matrix with values increasing lincarly between 0.6 and 1.4. In this case, the posterior density writes
p(0 = wi E1.,) w u2, E2,))
(41)
where for = 1,2, Mi, (3, + Im) n4 = Dr x, a.(x) = W(a; Im) and wi = wi/ Σ3 Wej. Therefore, the difffused marginals are
= (0 M1 + (1 ac)[m)
(42)
(W) / 12, + (1 at mm):
from which the score is
Vo, (x.ll( = (a) + (1 - 04) (n) - -a 1/2 m1)
(43)
(2( + (1 as) Im) 1(: - aj u22
1/2
D.3 Score of the diffused Log-Normal prior
The Log-Normal distribution easily be transformed into Gaussian distribution:
LogNormal(m, s) log O N(m,s) s).
We can therefore directly apply the Gaussian case from the previous section.
D.4 Score of the diffused Uniform prior
Consider the case where the prior a Uniform distribution A(0) Ld(0;, a, b), with (a, b) E Rm R", the lower and upper bounds respectively. Assume that (0, | 0o) No(;;; vanlo, up veIm); see Section with Ut 1 as.. We get the following analytic formula for the diffused prior:
/ (0(00) 0o)d0o
(0.)
(44)
t (be - as) / (a1,b))x-- xxx h...
1
N (Oes varlo, vmIm)
(45)
m
I
( (Vanbs $ (Vara;00ss,
(46)
Vau I (b: - ax)
where b2, the ith coordinate of 6g and variance a2. The score of the above quantity is then simply obtained by computing the score of each one- dimensional element the above product. For i m], the ith coordinate of Ve, log lo (0.) can be written
is the c.d.f. of univariate Gaussian with mean and
as Vo., log.ooffooo f(0.,)) = Vo,, ((0,.) with
1(0,,))
[(0r,)) = (D(aa Oes vi) (Varas 02,; u))
(47)
Vo., f(,,)) =
(W 01,,, vz)) - N (Vara,; & v)))
(48)
Varer
25

---

Published in Transactions on Machine Learning Research (02/2026)
E Experimental Setup
Code and compute resources. All experiments are implemented with Python (Python Software Funda- tion 2017) combined with PyTorch (Paszke et al. 2019). The code to reproduce all numerical experiments This includes specific information about the code environment and installation requirements, mentioned in the mdd file. In this url, we also provide pregenerated data and precomputed results that are necessary to quickly generate all figures of the paper without extensive computation time.
is provided in the following repository: https://gittin hwt.///aiiaaiiimaaraaaaffssiiooa
Methods. We choose F-NPSE GGefnne et al.) 2023) as a baseline for comparisons, which uses unadjusted Langevin dynamics with I. = 5 Langevin steps and a step size of 6, = 7(1 - o/)/ with r 0.5, to approx- imately sample from the tall postcrior. We do not compare to other SBI methods basod on NPE/NLE/ /RRE, as the score-based framework has shown to demonstrate clear superior performance (Sharrock et al. 2022: Geffner et al. 2023 Gloeckler et al.] 2024). Of course, this choice reduces the range of comparisons of our experimental section, but allows us to focus solely on the best alternative method from current literature.
For GAUSS and JAC, we infer the (all posterior by plugging the approximate tall posterior scores obtained via Algorithms [ and into the DDIM sampler defined in Section 2.1 Specifically, we sample from the backward Markov chain from Equation (34 with of 7 (1 a a111)/e(11 - az/a,-1) where rj = 0.2, 0.5, 0.8, for a number of steps T 50, 150, 400, 1000 respectively. We use a uniform scheduling {t: i/T)
The score network. Except for the two toy models considered in Section 4.1 the posterior scores have to be learned via NPSE. The same score model architecture and training hyper parameters are used for all tasks. Our implementation of the score model (0, x. t) is an MLP with layer normalization and 3 hidden layers of 256 hidden features. It takes as input the variables (0,x,1) I) and outputs a vector of the same size as 0 E RT. Before being passed to the score network, a positional embedding is computed for t E [0, [0.1] as per: (cos(tix), sin(ti))) isies, where we chose = 3 for the number of frequencies. Optionally, 0 and æ can each be passed to an embedding network. Note that for the benchmark experiment in Section 4.2 embedding networks were used, as we chose to use the same score model architecture across all tasks, each with different data dimensions. For the more complex neuroscience application in Section 4.3 0 and x were each embedded with a MLP with 3 hidden layers of 64 hidden features and an output dimension of 32.
For more stability, following Song et al. (2021a) and Ho et al. (2020): we actually learn the noise distribution -V s (0, x, f), with v, = 1 a the variance of the VP forward kernel, defined in Section 2.1 Given training dataset with N Nuran samples (Q,,i, Xs) p(0, x), the empirical loss function for the noise predictor' NPSE writes:
22,) Ea, which is related to the score function via (0, f)
1 N
CNPSE (d) = Σ (O.,, X2, 2,113, with O₁,, = Var, 0o,t + Ut, Z,,
(49)
where Z; ~ (O. and t; (0, 1). Training is done using the Adam optimizer over 5 5.00 epochs and early stopping on 20% of the training data. Note that early stopping requires to be done with care, since the loss function is stochastic. For each example and training set Nurin we trained several models for different learning rates and batch sizes, and selected the one that corresponds to the best trade-off between smallest validation loss and latest stopping epoch, as detailed in Appendix FF The training dataset is always normalized to zero mean and standard variance (for variables 0 and z). The same transformation has to be applied to the observations ri., and the prior score function Ve log (0) during sampling. After sampling, the inverse transformation is applied to the obtained samples of the approximate tall posterior distribution.
26

---

Published in Transactions on Machine Learning Research (02/2026)
Validation. To evaluate the accuracy of the sampling algorithms, compute the distance between 10* samples of the estimated and true posteriors considering three different metries: the sliced Wasserstein (sW) distance from the POT python package with 102 projections, the Maximum Mean Discrepancy (MMD) and the Classifor--To--ammpp Test (C2ST), both from the sbibm package, with default parameters.
The true posterior samples are obtained via MCMC using numpyro with jax, except for the Gaussian Mixture Toy Model from Section 4.1 for which we consider the Metropolis-Adiie Adjusted Langevin Algorithm (MALA);
see (Roberts & Tweedie 1996..
For the JRNMM in Section 4.3 no samples from the true posterior are available and the above metrics cannot be computed. We therefore use the local and the implementation provided in the sbi toolbox, with default parameters.
Test (E-C2ST) (Linhart et al. 2023
27

---

Published in Transactions on Machine Learning Research (02/2026)
F Loss functions, training strategy and model selection
In this section we explain our model selection strategy, namely the choices of learning rate and batch size used for the training of the score models in the experiments from Sections 4.2 and 4.3 As detailed in Appendix F. the model architecture (i.e. number of MT.P layers, hidden features, etc.) is fixed and training is done using the Adam optimizer over 5 000 epochs and early stopping using 20% of the training data to compute the validation loss. Note that carly stopping requires to be done with care, since the loss function is stochastic. For each example, we therefore trained several models for different learning rates and batch sizes, and selected the one that corresponds to the best trade-off between smallest validation loss and latest stopping epoch.
Figure 6 displays the train and validation losses for all SBI benchmarks examples and training set sizes Nurain For small Nerain, over fitting bohavior can be detected, which is particularly visible for the SLCP example. In these cases, a higher batch size of 256 is chosen, since it yields more stable loss functions, with delayed over fitting and thus also delayed early stopping. The learning rate is then chosen correspond the smallest validation loss when the gap is clear (e.g. for SLCP). In cases best validation loss in each sctting is similar and the learning rate then be chosen to correspond to the latest early stopping epoch (to mitigate varialions due to the stochasticity of the loss function and ensure longer training).
obvious overfitting is detected, the
The same strategy was then used to select the best model for the experiments on the JRNMM, for which the loss functions are displayed Figure 7 Here, the score models were trained on a large training set of size Nurain = 50 000. We therefore directly chose larger batch size of 256 for more stability, leading to equal early stopping times in cach setting. The chosen models then correspond to the smallest validation loss, here obtained for learning rate of 1
28

---

Published in Transactions on Machine Learning Research (02/2026)
Nunin 1000
Menin 3000
Nuai 10000
Nurai 30000
SLCP
0.3
0.2
0.1
0.0
Lotka-Volle Aakka-Vottera
0.4
0.3
0.2
0.0
0.4
train val (lr=0.0001)
SIR
0.3
train val (Ir=0.001)
0.2
0.1
2000
4000
0
2000
4000
0
2000
4000
2000
4000
epochs
epochs
epochs
epochs
(a) Batch size 64.
Miin
1000
Irann 3000
Nirin 10000
Niria 30000
0.5
0.4
SLCP
0.3
0.2
0.1
SIR Lotka-Voltonn Volterra
0.3
0.2
0.1
0.0
0.4
train
val (lr=0,0001)
train
val
(Ir= 0.001)
0.0
2000
4000
2000
4000
0
2000
4000
2000
4000
epochs
epochs
epochs
epochs
(b) Batch size 256.
Figure 6: Train and validation losses for the SBI benchmark examples obtained for score networks trained over 5 5000 epochs with the Adam optimzer and learning rates le-3 (orange) and Ir 1e-4 (blue). The two Figures respectively corresponds to batch size of 64 and 256. The validation loss was computed held-out validation set of 20% of the training dataset of size Nerain E [103,,.10,, 3.103. Dashed lines indicate the epoch at early stopping time. A higher batch size of 256 is chosen, since it yields more stable loss functions, with delayed over fitting and thus also delayed early stopping (particularly visible for the SLCP example). The learning rate then chosen to correspond to the smallest validation loss or, in cases were no obvious overfitting is detected, the latest early stopping epoch (to mitigate small variations due the stochasticity of the loss function and ensure longer training).
29

---

Published in Transactions on Machine Learning Research (02/2026)
3D
4D
0.5
train val (Ir=0.0001)
0.4
0.3
train val (Ir=0.001)
0.2
0.1
0.0
0
2000
4000 0
2000
4000
epochs
epochs
Figure 7: Loss functions for the Jansen and Rit Neural Mass Model obtained for score networks trained over 5 000 epochs with the Adam oplimzer, batch size of 256 and learning rates Ir = le (orange) and Ir = Le (blue). The validation loss was computed on a held-out validation set of 20% of the training dataset size Ntrin = 50 000. Dashed lines indicate the epoch at carly stopping time. As carly stopping in both settings occurs at the same time, the chosen model corresponds to the smallest validation loss, here obtained for learning rate of le
30

---

Published in Transactions on Machine Learning Research (02/2026)
G
Additional results for the toy models
Table 4 displays the results obtained for the GMM example, which are comparable to the oncs from Table 1 used to define the equivalent time setting in Section
Algorithm
N steps
Á (s)
sW
GAUSS
50
0.99 +/- +/- 0.00
0.30 +/- 0.05
JAC
50
0.89
0.00
82.73 +/-67.41
Langevin
50
1.77 4/-
0.21 4/-
GAUSS
150
1.83 0.00
0.31 4/1+//-00 0.04
JAC
150
2.66 +/-
4/- 1.07
Langevin
150
5.28 0/0
+/-- 0.02
GAUSS
400
3.91 +/-0.01
0.31 4/- 0.04 0.04
JAC
7.13 +/-0.02
3/-1 1.04
Langevin
400
14.06 0.0
0.43 4/- 0.02
GAUSS
1000
8.85 0.03
0.31 +/-0.33 0.03
JAC
1000
17.69 +/-
7.33 +/- 5.41
Langevin
1000
35.08
+/- 0.15
0.47 +/- 0.02
Table 4: Sliced Wasserstein (sW) distance and total sampling time At for the GMM toy problem with m = 10, n 32 and 10-2 and number T of sampling steps. Mean and std over 5 different seeds.
g and results obtained for the Gaussian example across all considered dimensions m [2, 4, 8, 10, 16, 32), number observations n, and perturbations e. Figure confirms the results of Figure 2 and additionally shows that JAC is less robust higher dimensions than GAUSS: for m 4, it yields sW values outside the plot boundaries (and LANGEVIN already for m > 2). Figure shows the same results but compares the results obtained for the different noise levels c, for each algorithm separately. For one, it highlights the precision of JAC in the non-perturbed case (e 0) and its instability otherwise. On the other hand, it emphasizes the superior robustness of GAUSS to noise, especially for medium noise level (e = 0.01) and in high dimensions.
To complete the robustness analysis, we include in Figures
31

---

Published in Transactions on Machine Learning Research (02/2026)
e 0.0
<= 0.001
c 0.01
0=0.1
1.00
0.75
m=
sW
0.50
0.25
0.00
1.00
0.75
sW
0.50
0.25
0.00
1.00
0.75
sW
0.50
0.25
1.00
0.75
10
sW
0.50
0.25
0.00
1.00
0.75
M 16
SW
0.50
0.25
0.00
1.00
0.75
m= 32
sW
0.50
GAUSS (1000)
0.25
JAC (400)
LANGEVIN (400)
0.00
10
102
10
102
10
n
n
N
Figure 8: Sliced Wasserstein (sW) distance function of the number of observations for each algorithm
(GAUSS, JAC and LANgEVIN) and with different levels of Results are shown for the Gaussian example in several dimensions E [2, 4, 8, 10, 16, 32]. Mean and std over 5 different seeds.
32

---

Published in Transactions on Machine Learning Research (02/2026)
1.00
2
m
#0
16
m1 =
0.75
GAUSS
SW
0.50
0.25
0.00
1.00
M M
M
JAC
SW
0.50
0.25
0.00
LANGEVIN
0.75
SW
0.50
000
0.001
0.25
0.00
'00
10a 10a 10a
10a
10a
10a
10a
10
10a
Figure 9: Sliced Wasserstein (sW) distance as a function of the number of observations n for different levels of e, for each algorithm (GAUSS, JAC and LANGEVIN). Results shown for the Gaussian example in several dimensions E 2,4, 8, 10, 16, 32). Mean and std over 5 different seeds.
33

---

Published in Transactions on Machine Learning Research (02/2026)
H Additional results for SBI benchmarks
We here include more dctailed results for the three SBI benchmark tasks from Section 4.2
As a complement
to the sliced Wassersetin (&W), we report results obtained for the Maximum Mean Discrepancy (MMD) and
the Classfien-TTo--Samppll Test (C2ST) metrics from the sbibm package, respectively in Appendix H.2 and
H.3
The results for sW are also included in Appendix H.1
Tables and r summarize the results obtained for all metrics, for the three benchmark tasks respectively, In bold are displayed the smallest values for the considered metric and number of observations n. This highlights that the approximate posterior samples obtained with our methods (GAUSS and JAC) are closest to the true posterior samples, in particular for n > 1. In red we mark the highest value for each metric across all methods and number of i..d. observations n. This highlights that LANGEVIN performs significantly worse, especially as the number of observations n increases (n 30), and this for all considcred metrics.
The results in this section highlight an interesting point: increasing distance values (i.e. degradation of performance) for higher n, which is particularly well shown for the MMD and C2ST metric in Figures [4| and 16 respectively. This can be explained by the accumulation of approximation errors, as sum over n evaluations of the score model to obtain the tall posterior score for observations. We refer the reader to Appendix
L.2 which investigates possible solution to this issue: partially factorized score-bases posterior sampling algorithms, such as PF-NPSE proposed by (Geffner et al. 2023).
The lables also point out the differences between the metrics, each providing complementary information about the statistical validity of the approximate posterior. For instance, sW known capture well the geometric features of the underlying distribution but struggles in higher dimensions. MMD allows to analyze correlations in higher dimensions, but is known to fail for example in modiil settings. C2ST is a measure that captures all kinds of inconsistencies, but this means that it can quickly yield very high values and can be "too discriminative". It also means that C2ST is less informative if the goal is to analyze specific statistical features. The three benchmark examples each have different posterior structures, as shown in Figure 10 which enables to illustrate these differences in a more detailed analysis of the results from Tables 5 and
SLCP Table 5 This is a hard example for which the postcrior has 4 modes a 5-dimensional paramcter space. The general difficulty of this task is shown with C2ST, with high scores across all methods and mumber of observations n (C2ST= 0.7 for 1 for all methods and increases rapidly for > 1). Here MMD is uninformative (low scores everywhero), while sW effectively coptures the approzimation errors due to the multi-modal structure of the posterior.
Lotka- Volterra (LV)- Table G The posterior has a single mode and lives in a 4-dimensional parameter space. The difficulty of this task is to be accurate in every dimension and to concentrate around the true parameters as increases. Again, C2ST is high for all methods and all n. MMD completes the analysis by effectively capturing the accumulation of errors as n increases. This time sW is uninformative (small scores everywhere) as there is no specific geometric structure like multi-modality in the SLCP example.
The posterior has only one mode and lives in a 2-dimensional parameter space. This task is simpler as shown by the C2ST scores that start low ( 0.6 for n = 1). Like in IV, MMD reflects the accumulation ervor as n increases and s W is uninformative given the lack of difficult geometry in the posterior distribuution
SIR Table 7
34

---

Published in Transactions on Machine Learning Research (02/2026)
SW
MMD
C2ST
14
30
14
30
1
14
30
GAUSS
0.31
0.37
0.40
0.02
0.02
0.02
0.70
0.92
0.95
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.10
0.11
0.13
0.02
0.02
0.01
0.08
0.06
0.05
JAC-clip
0.31
0.47
0.51
0.02
0.02
0.02
0.70
0.89
0.94
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.09
0.08
0.10
0.02
0.01
0.01
0.07
0.07
0.06
LANGEVIN-clip
0.31
0.56
0.81
0.02
0.05
0.10
0.70
0.92
0.97
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.11
0.17
0.26
0.02
0.05
0.07
0.08
0.06
0.03
Table 5: Results for SLCP for Nuain 3.10
Mean and std over 25 parameters 0 X(0).
s sW
MMD
C2ST
71
14
30
1
14
30
1
14
30
GAUSS
0.005
0.004
0.005
0.05
0.20
0.31
0.77
0.93
0.96
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.004
0.005
0.005
0.07
0.23
0.21
0.10
0.04
0.03
JAC-clip
0.005
0.01
0.007
0.05
0.26
0.40
0.76
0.94
0.97
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.004
0.007
0.005
0.07
0.47
0.38
0.10
0.07
0.04
LANGEVIN
0.006
0.03
0.38
0.15
0.61
0.63
0.85
0.99
1.00
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.005
0.02
0.18
0.10
0.23
0.12
0.07
0.02
0.00
Table 6: Results for LV for Nurain 3.10 Mean and std over 25 parameters l ~~ X(0).
S W
MMD
C2ST
71
14
30
1
14
30
1
14
30
GAUSS
0.001
0.002
0.002
0.007
0.23
0.38
0.57
0.78
0.86
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.001
0.001
0.001
0.007
0.16
0.25
0.03
0.07
0.07
JAC-clip
0.001
0.002
0.002
0.006
0.08
0.21
0.58
0.75
0.84
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.001
0.001
0.001
0.005
0.09
0.09
0.03
0.07
0.06
LANGEVIN
0.001
0.01
0.07
0.018
0.56
0.73
0.63
0.93
1.00
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
+/-
0.001
0.006
0.02
0.018
0.25
0.03
0.04
0.12
0.00
Table 7: Results for SIR for Ntrii 3.10
Mean and std over 25 parameters 0* ~~ A(0).
35

---

Published in Transactions on Machine Learning Research (02/2026)
SLCP
8
8
V
8
V
8
g
WV
W
AA
A
AA
LV
6
SIR
TRUE
GAUSS
Figure 10: Pairplots for SLCP, Lotka Volterra (LV) and SIR. 1D and 2D marginals for the true (green, left) vs. approximate tall posterior samples obtained with CAUSS (blue, right) for a score model trained on Ntrain = 3 samples and for n = 1, 14. 30 (dark to light colors) i.i.d. observations simulated for a given parameter (black dots and dashed lines).
36

---

Published in Transactions on Machine Learning Research (02/2026)
H.1
Sliced Wasserstein distance
1
8
14
22
= 30
2.0
GAUSS/-clip
cAC-c1ip
SLCP
sW
1.0
LANGEVIN/-CITP
0.5
0.0
Lotka-Vol Volterre
0.04
SW
0.03
0.02
0.01
0.00
0.008
SIR
sW
0.006
0.004
0.002
0.000
10
10
10a
10
10a
10
10a
10
10
10a
Nunin
Numi
Nuran
Nirain
Nirain
Figure 11: Sliced Wasserstein (sW) as function of Nurai E [ 3.10 10 3.10 between the samples obtained by each algorithm and the true tall posterior distribution p(0 ai,n) (for E [1,8, 30). Mean and std over 25 different parameters 0* ~ A(0).
Nirain 1000
Nurain 3000
Nurain 10000
Nirain 30000
2.0
GAUSS/~clip
1.5
JAC-c11p
SLCP sW
LANGEVIN/-clip
1.0
0.0
Lotka-Voltoen Volte-Vol sW
0.04
0.03
0.02
0.01
0.00
0.008
SIR
sW
0.006
0.004
0.002
0.000
14
30
22
301
14 71
22
30
8
30
7
Figure 12: Sliced Wasserstein (sW)+ function of [1,8, 14, 22, 30] between the samples obtained by each algorithm and the true tall posterior distribution p(0 xi..)) (for Ntrain E [103,3 3.103, 104, 3.10%)). Mean and std over 25 different parameters 0" ~~~~A0 ~ A(0).
37

---

Published in Transactions on Machine Learning Research (02/2026)
H.2
Maximum Mean Discrepancy
n=1
n=
n= 14
12 = 22
#=
GAUSS/-CI1p
SLCP
MMD
0.4
JAC-clip
LANGEVIN/-cLip
0.2
0.0
Nolka-Volterra MMD
1.0
MMD
SIR
0.2
0.0
10
10
10*
10
10a
10a
10
10
10'
Muain
Nusin
Narain
Nurin
Nusin
Figure 13: MMD as function of Nuran [10*,3.10" 10 3.10*] between the samples obtained by each algorithm and the true tall posterior distribution p(O xi..) (for E [1, 14, 22, 30)). Mean and std over 25 different parameters A(0).
Mirain 1000
Mirain 3300
Nurin 10000
Ninin 30000
0.6
GAUSS/ GAUSS/-Clin -clip
SLCP
MMD
0.4
JAC-clíp
LANGEVIN/-CL1p
0.2
0.0
Lotka-Voll Volterra MMD
1.0
0.5
0.0
0.8
SIR MMD
0.6
0.4
0.2
0.0
8
14
22
30 1
8
14
22
30 1
8
14 11
22
30
14
22
30
n
11
f2
Figure 14: MMD as a function of n E [1,8, 8, 14, 22, 30) between the samples obtained by each algorithm and
the true tall posterior distribution p(0 xi,.) (for Ntrai E [103, 3.103,, 10 (.1100)))). Mean and std over 25 different parameters 1(0).
38

---

Published in Transactions on Machine Learning Research (02/2026)
H.3
Classifier- Two-Sample Test
R=
71
8
R2
14
n= 22
n==30
1.0
SLCP
C2ST
GAUSS/-clip JAC-cl1p LANGEVIN/-Clip
Lotka-Vollerra
1.0
Lotka-Voltern
C25 225
0.8
0.6
1.0
C2ST
SIR
10a
10
10a
10*
10a
10a
10a
10*
10a
10*
Nurin
Nirain
N
Nuain
N Nirain
Figure 15: C2ST as a function of Nuain [10 3. 10 104, 3.100 between between the samples obtained by each algorithm and the true tall posterior distribution p(0 xi,n) (for over 25 diffcrent parameters 8* X(0).
[1,8, 14, 22, 30)). Mean and std
Nuain 1000
Muain 3000
Nirain 10000
Ntrain 30000
1.0
SLCP C2ST
0.8
GAUS5/-cl1p
JAC-clip
0.6
LANGEVIN/-C115
Lotka-Volle Voltera
1.0
C2ST
0.8
0.6
1.0
C2ST
SIR
0.8
0.6
30
14
22
22
30
n1
71
71
n
Figure 16: C2ST as a function of n E [1,8, 8, 14, 22, 30] between the samples obtained by each algorithm and the true tall posterior distribution p(0 =i,..) (for Nivan E [103, 3,.00, 10 1 3.10")) and the Dirac of the true parameters 8* used to simulate the observations zi,-- Mcan and std over 25 different paramoters 8 X(0).
39

---

Published in Transactions on Machine Learning Research (02/2026)
Results for additional tasks from the SBI benchmark
To complete our empirical study, we added results for additional examples from the SBI benchmark (Lueck- mann et al. b: Gaussian Linear, GMM, GMM (uniform) 8 B-GLM/ (raw) and Two Moons. These new results allows us to compare the performance of our proposal on other challenging situations, such as when scaling to highly structured (e.g. multimodal) posteriors and high-dimensional observation spaces. Note that these examples go a step further as compared to the experiments carried out by Geffner et al. (2023). including non-Gaussian priors Figures 17
18 and 19 respectively report the sW, MMD and C2ST function of Ntrain- They all show that our algorithm outperforms the Langevin sampler, with smaller dis- tance values everywhere but for the GMM (uniform) example. Interestingly, the certainly very discriminative C2ST metric, shows that all three algorithms fail to infer the tall posterior for the Two Moons and B-GLM/ (raw) examples.
8
14
22
30
0.8
GAUSS/-clip
Gaussian
JAC-c1íp
LANGEVIN/-C11P
0.2
0.0
0.8
GMM
0.6
0.4
GMM (uniform)
1.5
19
1.0
0.5
B-GLM
0.2
0.0
B-GLM (raw)
0.8
0.6
Two Moons
0.3
0.2
0.1
0.0
10
105
10a
10a
10a
10a
10a
109
10a
109
Nirain
Mirain
N Mizain
Nirain
N Mirain
Figure 17: Results for additional benchmark examples. Sliced Wasserstein (sW) distance as a function of Nirain [10", 3.10 10a, 3.10 between between the samples obtained by each algorithm and the true tall posterior distribution p(0) zT.,.) (for [1,8, 22, 20]). Mean and std over 25 different parameters 0* ~ X(0)
Saam CHM but with Uniform prior.
BBrnoulll GLM with summary statics high dimensional raw data. 1oNote that this was already the case for SLCP with Uniform prior.
40

---

Published in Transactions on Machine Learning Research (02/2026)
n=8
A1
==30
1.5
Gaussian
JAC-clip LANGEVIN/-CI1P
MMD
0.5
GMM
MMD
1.0
0.5
0.0
GMM (uniform)
1.5
B-GLM
MMD
B-GLM (raw)
MMD
1.0
0.5
Two Moons
MMD
0.0
10
10a
10
10a
10
10a
10a
10a
10a
Nirain
Nirai
Nurin
Nuran
N
Figure 18: Results for additional benchmark examples. MMD as a function of Ntrain E 3110, 104,3.105 between between the samples obtained by each algorithm and the true tall posterior distribution p(0 zt..) (for n E [1,8, 14, 30).. Mean and std over 25 different. parameters 0* X()).
41

---

Published in Transactions on Machine Learning Research (02/2026)
11
n ==
Gaussian
C2STT
0.8
1.0
GMM
C2ST
0.8
GMM (uniform)
1.0
C2ST
0.6
B-GLM M
C2ST
0.8
B-GLM (raw)
1.0
C2ST
0.8
Two Moons
1.0
C2ST
GAUSS/-CL1p
JAC-ci1o
LANGEVIW/-clip
0.6
10a
10a
10a
10a
10a
10a
10a
10a
10a
104
Nerain
Nizain
Nuran
Nirain
Nurain
Figure 19: Results for additional benchmark examples. C2ST as a function of Ntrain E [_ 3.30 between between the samples obtained by each algorithm and the true tall posterior distribution p(0 zt..) (for n E [1,8, 14, 30)). Mean and std over 25 different. parameters 0* A()).
104,3.10
42

---

Published in Transactions on Machine Learning Research (02/2026)
J. The Jansen and Rit Neural Mass Model (JRNMM)
J.1 The JRNMM as a system of stochastic differential equations
The JRNMM (Jansen & Rit 1995) serves as an illustrative example from computational neuroscience in Section 4.3 Specifically, we consider the implementation by Buckwar et al. (2019) of the stochastic version of the JRNMM presented in (Ableidinger et al., 2017). In our experiments this model simply presented as a black-box simulator with four input parameters represent the EEG measurement of brain signal s(I). We will now give insights into the inner workings of this simulator, which behind the generation of s(t), following the description provided in (Rodrigues) ct al. 2021). The stochastic JRNMM models the interaction between excitatory and inhibitory neuronal populations in the cortical column, described by a system of three coupled non-linear stochastic differential equations, that can be written as six-dimensional first-order stochastic differential system:
(C, p4, o, that outputs signal x(t) 10 10, which
xo()) = Xx(t)
x1()) Xa(t)
X2(1) Xo(1)
(50)
Xz(1) (Aa (ug + Sigm (Xi(1) - X2(1)) - 20Xz(1) - x2 (1)
x1(1) = (Aa (Mg + C2 Sigm (C)X,(1)) - (1xx(1) - -2xx(1))+ + +,,(((1)
Xs(t) (Bb (Ms C4 Sigm (C)Xo(t)) ( - -b2x2(t) +otWi(()
The observed EEG measurement then corresponds to xx(() 100/20(X1()) (x(())xx(()) where g is a gain factor expressed in decibels. According to Jansen & Rit 1995). most of the parameters in 50 are expected to be almost constant across different individuals and different experimental conditions, except for (Ci, C2, C2, C4), that represent connectivity, and μ4, T4, the statistical parameters of the input signal from neighboring cortical columns. Based on this assumption, [Buckwar al. (2019) propose a simplified implementation of the model defined by a reduced set of four parameters
(C, where p f44, σ4 and all connectivity parameters are related via C1 C, C2 0.8 C, Css = 0.25, C1 0.25. All other parameters in 50 are fixed at their "constant" value.
J.2 Inference for the JRNMM
This section details some high level choices to perform inference on the JRNMM that is used to illustrate proposal real-world example. They are exactly the same used in the experiments from the original HNPE paper (Rodrigues et al. 2021).
Prior distribution. To avoid any bias due to misspecification issues, inference is often performed with simple and rather uninformative prior. This ensures that the parameter space is sufficiently explored to provide reliable results. In this work, the prior distribution is chosen to be Uniform over the range of scientifically plausible parameter values (Rodrigues et al. (2021) Buckwar et al. (2019) Ableidinger et al.
2017 Jansen & Rit 1995
C' ~ 2510, U(50, 500), u(0,5000), ~ u=((-0,+2)
(51)
Summary Statistics. In traditional SBI methods, is standard to use summary statistics to improve the inference quality by reducing the dimensionality while describing sufficiently well the statistical features of the data. When using neural density estimators, it is possible to learn these summary statistics with specified embedding networks, such as LSTMS for time scries data (Rodrigues & Gramfort 2020). However, for comparison purposes and enhanced clarily of our ontribbutions we stick to the summary statistics traditionally used for neural time series data (Buckwar et al. 2019). They consist of the logarithm of the power spectral density (PSD) of each obscrved time scries. The PSD is evaluated in 33 frequency bins between zero and 64 Hz (half of the sampling rate). This leads to setting with parameters estimate given observations defined in 33-dimensional space.
43

---

Published in Transactions on Machine Learning Research (02/2026)
J.3 Validation results with C-C2ST
In this section, report results on the validity of the inferred posterior for the simplified version (3D) and the full (4D) Jansen and Rit Neural Mass Model (JRNMM).
In both cases, we evaluate the accuracy of cach posterior sampling algorithms (GAUSS, JAC and LANGEVIN) w.r.t. the true tall posterior, using the local Classfiier-Tww-samppp. Two-sample Test (L-C2ST) diagnostic proposed by Linhart et al. (2023) and implemented in the toolbox. E-C2ST is frequentist hypothesis test that evaluates the null hypothesis of sample equality between two conditional distributions by training classifier data from the respective joint distributions. The test statistic takes values in [0, 0.25] and is defined by the mean squared error (MSE) between the predicted class probabilities and one half, i.e. the chance level where the classifier is unable to distinguish between the two data classes. Higher MSE values indicate more distributional differences.
Figures 20a and 20b respectively show the E-C2ST results for the 3D and 4D JRNMM cases, obtained using the MLPCCassifii from scikit-learn, with default parameters from the abi implementation and trained using 10000 samples (,, from the joint distributions corresponding to the estimated and true tall data posteriors. We report the mean and std over different seeds used for the initialization of the classifier. On the left, we plot the E-C2ST statistics, computed on 10000 samples of the approximate posterior obtained by each sampling algorithms for a given observation set Fin ~ p( 0*). On the right, we display the associated p-values computed by comparing the obtained statistics to the mull hypothesis, estimated over 100 trials. In both cases, the GAUSS algorithm consistently yields test statistics close to 0 and p-values above the significance level set at a 0.05, which means that there is no evidence that the estimated posterior is not statistically consistent with the true tall posterior at vims This is not the case for JAC, whose un-clipped version yields MSE values outside the plot limits and LANGEVIN, for which the test is more easily rejected in the 3D case (large std on pvalues that overlap with the significance level) and very clearly rejected in the 4D case for > 14 (p-value below the significance level).
CCCSST statistics
p-values
1.0
0.2
GAUSS/-clip JAC-clip
0.5
0.1
LANGEVIN/-clip
significance level
0.0
0.0
1
8
14
22
30
8
14
22
30
n
n
(a) Inference on the JRNMM (fixed g 0).
CCCST statistics
p-values
1.0
0.2
GAUSS/-clip
JAC-clip
0.5
0.1
LANGEVIN/-Clig
significance level
0.0
0.0
8
14
22
30
8
14
22
30
(b) Inference on the full (4D) JRNMM.
Figure 20: Accuracy of the sampling algorithms w.r.t. the true tall posterior. We show (-C2ST statistics (left) and corresponding p-values (right) computed for samples obtained with GAUSS, JAC and LANGEVIN at zi. for [1,8, 3, 14, 30). Mean and std over seeds.
44

---

Published in Transactions on Machine Learning Research (02/2026)
K Limitations of Langevin sampling
Sensibility to the quality of the score model. Our results in Section 4.1 analyze the robustness of the different sampling algorithms and essentially show that the Langevin algorithm is very sensible to noisy score networks.
Step-size choice. We found that different step sizes can generate very different results for given score model. Figure 21 shows comparison between LANGEVIN the unadjusted Langevin algorithm (ULA) from (Geffner et al., 2023) used in all experiments from Section 4 - and tamed ULA, additional implementation with tamed step sizes from (Brosse et al., 2017] as means to stabilize ULA. We can see that LANGEVIN quickly diverges as increases. A possible explanation could be a learning rate that is too large for settings with big n (i.e. not enough steps are done). We can see that the stabilization tools from the tamed version yield a more stable ULA algorithm, but does not provide a satisfying solution either (lower sW, but still diverges). Fundamentally, there exists a setting where the Langevin algorithı will work (small enough learning rate, run for a long enough time), but this setting extremely dependent of the problem hand. This is precisely the strength of our algorithm compared to ULA: we do not need to sample several times for each marginal pt at each time step t. Note that unfortumately, the code for (Geffner et al.) 2023 is not available, so our results are based in best-elfort attempt to reproduce the proposed algorithm.
n
n
n =
n
2.0
LANGEVIN (clip) tamed ULA (clip)
1.5
SLCP SW
1.0
0.5
0.0
Lotka-Volf Volterra
0.04
SW
0.03
0.02
0.01
0.00
0.008
SIR
sW
0.006
0.004
0.002
0.000
10
10
10
10a
10
Nusin
Nuzin
Nuran
Nuvin
Nusin
Figure 21: Comparison between the LANGEVIN algorithm from Geffner et al. 2023 (used in all our experiments) and more stable tamed ULA version with tamed step size from (Brosse et al. (2017). The plots show the sliced Wasserstein (sW) w.r.t. the true tall posterior p(0 aim) as function of Nemin E (103, 3.103,, 10 3.10*) and for E {1, 8, 14, 20,
45

---

Published in Transactions on Machine Learning Research (02/2026)
L
Extensions: classifier-free guidance and partial factorization
The implementation of both of these extensions can be found in our Code repository: https://github. com/ Their performance was investigated on the three benchmarks from Lotka-Volterra, SLCP and SIR.
bullaliihharttiiffusssoooo--or
Section 4.2
L.1
Classifier-free guidance
It is possible to implicitly learn the prior score via the flassifier-free guidance (CFG) approach (Ho & Salimans )021), which essentially consists in randomly dropping the context variables when training the posterior score model (e.g. 20% of the time). This is useful in cases where the diffused prior score cannot be computed analytically. Figure [2 displays the sliced Wasserstein (sW) distance as a. function of Nirain and compares the results obtained for GAUSS with the learned vs. the analytical prior score (GAUSS (CFG) GAUSS). We also report the results obtained for the Maximum Mcan Discrepancy (MMD) and Classifier- Two- Swo-Sample Test (C2ST) accuracy in Figure 23b The results highly accurate for the Log-Normal priors of Lotka-Volterra and SIR, but less satisfying for the Uniform prior in SLCP. We think that this is caused by the discontinuities of the Uniform distribution. In summary it seems that, under some smoot assumptions, is indeed possible to learn the prior score via the classifier freeee guidance approach.
==8 8
n= 14
n= 22
== 30
2.0
SLCP
1.0
0.5
0.0
Lotko-Vol Voleraa
Volttrra
SW
0.03
0.02
0.01
0.00
GAUSS GAUSS (CFG)
SIR
0.006
0.004
0.002
0.000
10
10
10
10
10a
10a
10a
10
10
Nuwin
Nuain
Nuain
Nuui
Ninan
Figure 22: Results obtained for GAUSS with the learned vs. the analytical prior score (GAUSS (CFG) vs. GAUSS). sW w.r.t. the true tall posterior p(0) ii..) a function of Naran {103, 3.103, 10 3. 3.10*) and for {1, 8, 14, 22,30
46

---

Published in Transactions on Machine Learning Research (02/2026)
R 14
n = 22
n= 30
0.6
GAUSS
SLCP
MMD
0.1
GAUSS (CFG)
0.2
Lotka-Volterra
MMD
1.0
0.0
SIR
0.0
10
10
104
10a
104
10a
10a
1004
Nirain
Nirain
Nirain
Nirain
Nirain
(a) Maximum Mean Discrepancy (MMD).
n= 8
#=
14
n= 22
== 30
1.0
SLCP
C2ST
0.6
Kolka-Volterra
(2ST
0.8
0.6
SIR
C2ST
0.8
GAUSS
GAUSS (CFG)
10a
104
10a
104
10a
104
10a
10a
10a
10a
Nirain
Nirain
Nuaan
Nirain
Nuain
(b) Classifier-Twor Sample Test (C2ST) accuracy.
Figure 23: Results obtained for GAUSS with the learned vs. the analytical prior score (GAUSS (CFG)
vs. GAUSS). MMD and C2ST w.r.t. the true tall posterior p(0 | xi:n) as function of Nirain E {10 3.103, 104, 3.10 and for a E {1, 8, 14, 20.
47

---

Published in Transactions on Machine Learning Research (02/2026)
L.2 Partial factorization (PF-NPSE)
In the same way as in (Gellner et al. 2023 our proposed algorithm can naturally be extended to partially factorized version. Specifically, it consists in approximating the tall posterior by factorizing it over batches of context observations (instead of a single r). To do so, the score model is modified to take as input context sets with variable sizes (between 1 and mmax). The given sampling algorithm (e.g. GAUSS, LANGEVIN) is then modified to split the context observations "f:...,,"; into subsets of smaller size < 1mmax < n, before passing them the trained score model. This approach should allow for a good trade-off between the accumulation of approximation errors due to multiple evaluations of the score model (n/ times) and the increased simulation budget (xflma))- This approach should allow for a good trade-off' between the accumulation of approximation errors due to multiple evaluations of the score model (n/numax times) and the incrcased simulation budget ( (X
We investigated the performance of PF-NPSE on the SBI benchmarks (Lotka-Volterra, SIR and SLCP). For each of the three examples we trained PF-NPSE model targeting the score models for the law of 0 given Ul:mmax for 7max E (1,3,6 6, 30). Figure 24 displays the sliced Wasserstein (sW), Maximum Mean Discrepancy (MMD) and the Classifier-Twol Two-Sample Test (C2ST) accuracy for 30 observations a function of the Nerain for samples obtained with the partially factorized LANGEVIN and GAUSS samplers and for all maxx The extreme case Ramax = 1 corresponds to the original 'fully" factorized version of the samplers. Amax = 30 correspond to the her extreme case with no factorization, but maximum simulation budget. We can see that the optimal sW values lie in the middle of the spectrum (i.e. for mmax = 3, 6), which corresponds to what was concluded in (Geffner et al., 2023). Note that the performance of LANGEVIN drastically improved for mmax 1, while GAUSS all results are close. In any case, the results suggest that practitioner will gain in choosing (a small enough) mmaa > 1.
SLCP
SIR
2.0
0.05
0.010
1.5
0.04
0.008
SW
0.03
0.006
1.0
0.02
0.004
0.5
0.01
0.002
0.0
0.00
0.000
0.6
1.5
1.0
GAUSS, Paux =1
0.8
GAUSS, =
MMD
0.4
1.0
0.6
GAUSS, Maxs =
GAUSS, Mars= 30
0.2
0.5
0.4
LANGEVIN, Plexx
0.2
LANGEVIN, Maxx 3
LANGEVIN, Grsa
0.0
0.0
0.0
LANGEVIN, Rases 30
1.0
1.0
1.0
C2ST
0.8
0.8
0.8
0.6
0.6
0.6
10a
10
10a
10
103
104
Numin
Nuain
Nusia
Figure 24: Results obtained with the partially factorized LANGEVIN and GAUSS samplers to infer the tall
posterior conditioned on total number of observations n 30, for mmax 1, 6, 30. We report the sliced Wasserstein (sW), Maximum Mean Discrepancy (MMD) and the - Two-Sample Test (C2ST) accuracy w.r.t. the true tall posterior p(0 si.n) as a function of Ntrain E {103 ,,100, 3..10 The extreme case Rmmax = 1 corresponds to the original 'fully" factorized version of the samplers. 7mmax = n = 30 correspond to the other extreme case with no factorization, but maximum simulation budget. We can see that the optimal distance values lie in the middle of this spectrum (i.e. for 1a = 3, 6).
48

---

Published in Transactions on Machine Learning Research (02/2026)
M Comparison with the deterministic sampler from Geffner et al. (2023)
M.1 Theoretical comparison
The alternative sampler proposed in Geffner et al. (2023) (Appendix D) builds Markov chain for the tall posterior by composing the reverse kernels associated with each individual posterior. They follow the approach from Sock--ick et al. (2015) to approximately sample from product of distributions, cach associated to an independent diffusion process. Concretely, if
11 (0 0.) = N
(1( (1( (1 (1-
(52)
denotes the learned reverse kernel for zjj the proposed algorithm performs reverse step of the form
P N (Me Fim), I)
(33)
where 14 is the average over all individual means plus heuristic prior correction, and all observations share the same scalar variance of 1-2, of a single diffusion whose marginals match the tall posterior, but rather defines a surrogate Markov chain intended to approximate it.
[n-T- This transition is therefore not derived as the reverse process
Our approach instead samples backward from the single diffusion process whose marginals satisfy
Pe(Oe 21mn) x ((B0)² II p(o 13)) 4ejo( 0u) ddu
(54)
and then approximates its reverse process through a Tweedie approximation of each backward kernel
do0 (Co 0.) ~ N (14 E0,) with Ee,(O,) 1 at Vpr (0.),
(55)
az
where ()) (0.) = Vor (0, + (1 - a.) v. log (e x,)). This yields our approximate score
Σ Enj (0.) Vo, logpr 1z)) + (1 n) E7,) (0.) Vo, logp)(0))
Vo, logp:(0) xIm) R A(O.)-
(56)
with A(0) = Eド E (0) + (1-m) (0) and assumed constant covariances, as explained in Section 3.2
We therefore have the score of the marginal density of an actual known forward diffusion process and can thus run DDIM sampler along the correspouding reverse-time dynamics. Importantly, the backward covariance matrices in 90t: include important information about these dynamics (via. the score) and can be different for each individual posterior. In contrast, the sampler in Gelfner et al. (2233 uses fixed transition covariance I2,, enforcing isotropic contraction that does not reflect posterior geometry. We do acknowledge that more approximations come into play when we assume the covariances to be constant (and compute them following the GAUSS or JAC algorithms).. But they seem more subtle than just defining isotropic covariances as in the sampler from Geffner et al. (2023) (Appendix D) and adding prior correction term.
For example, take the case where each posterior and the prior are Gaussian (Appendix D). Our approximation recovers the exact tall-posterior precision (and score). The sampler proposed in Geffner et al. (0023 does not, unless all covariances are isotropic.
In short, both methods make Gaussian approximations when computing the backward transitions, but not at the same level. Ours is structurally closer to the true reverse diffusion, enabling DDIM sampling, correct posterior contraction, and exactness the Gaussian limit. Note that both methods can be equivalent if each individual posterior and the prior are Gaussian with isotropic covariance matrix (or the parameter space is one-dimensional).
49

---

Published in Transactions on Machine Learning Research (02/2026)
M.2 Empirical Comparison
We compare GAUSS, JAC, and LANGEVIN with the deterministic sampler from Gellner el al. (2023] (Appendix
D), referred to as _ET_ DET_GEF, on the toy models from Section 4.1 and the benchmark tasks from Section 4.2
The implementation for this new sampler and instructions to reproduce the following experiments can be found in our Code repository:
-or
Runtime. Table extends the speed-up comparison from Table for the Gaussian toy example. shows that DET DET_GEF is consistently the lastest sampler in our experiments, followed by GAUSS and JAC. For the remaining experiments, again consider the equivalent time setting with 400 and 1000 steps for JAC LANGEVIN and DET_GEF GAUSS respectively.
Speed up GAUSS
Speed up JAC
Speed up DET_GEF
2
0.39 + 0.04
0.45 + 0.00
0.22 + 0.00
0.37 + 0.01
0.45 0.00
0.22 + 0.00
0.37 + 0.01
0.46 + 0.00
0.22 = 0.00
10
0.37 + 0.01
0.45 + 0.00
0.22 = 0.00
16
0.37 : 0.01
0.49 +: 0.01
0.22 + 0.00
32
0.37 1 0.01
0.52 1 0.01
0.21 + 0.00
Table 8: Ratio of the runtime for GAUSS, JAC, and _ET_GEF different dimensions m. Averaged over the number of steps T E (50, 150, 400, 1000), different noise levels E {0, 10 10,1,01000000 2,102,,00 11100,,0011----, and the number of observations n E [L, 100]. Mean and std over 5 different seeds.
GEF w.r.t. LANGEVIN for the Gaussian example in
Accuracy and stability. Figure 25 shows the empirical results obtained for the Gaussian toy examples from Section 4.1 In the Gaussian example, GAUSS achieves the best performance across noise levels, as the Gaussian approximation exact in this case. JAC is also exact, but less stable when noise increases. DET_ remains close but shows a consistent gap, since relies on surrogate chain with fixed isotropic variance and is therefore not exact in this example where the covariance is not isotropic. As increases, DET_ _EF_ and LANGEVIN improve and become closer in performance, consistent with posterior concentration: the tall posterior becomes increasingly unimodal with small variance (Bernstein von Mises), making DET_GEF a better approximation and LANGEVIN easier to converge.
For GMM, LANGEVIN, JAC, and DET_GEF outperform GAUSS in low-noise regimes. In particular as grows, all methods become highly accurate, except GAUSS whose approximation gap remains visible. However, GAUSS remains the most stable method in the most noisy regime (e 0.1), where LANGEVIN and DET_ GEF degrade more strongly, and JAC diverges strongly, with sW values outside the plot limits.
We now move on to the bbi-henchmark tasks from Section 4.2 where the score is learned and imperfect, and compositional evaluation across observations causes error accumulation as n increases. Results are shown in Figure 26 In this regime, LANGEVIN and DET_GER diverge as Nusin grows: the learned score becomes more confident while retaining small systematic bias, which amplified by postcrior contraction. DET GEF is particularly sensitive since its approximation relies directly on the composed backward mean, which depends the learned score. Tn contrast, GAUSS more closely approximates the reverse diffusion dynamics of the tall posterior and remains significantly more stable, consistent with its robustness in the noisy GMM toy regime, JAC with clipping is promising but still biased.
Overall, these additional experiments confirm our main conclusions: LANGEVIN is the computationally most expensive and least stable method. DET_GEF provides useful fast deterministic (i.e. that does not rely on Langevin dynamics) baseline, but GAUSS, and stable variants of JAC, offer the most reliable sampling behavior in realistic learned-score tall-data. benchmarks.
Note that the results for GAUSS, JAC, and LANGEVIN may slightly differ from those reported in the main text, the experiments rerun from scratch and minor variations can arise from different, random seeds and GPU configurations.
50

---

Published in Transactions on Machine Learning Research (02/2026)
c=0.0
c=0.00
c=0.01
e= 0.1
1.00
Gaussian
0.50
DET 10000
0.25
GAUSS 1300;
440)
0.00
LANGEVIN -400
1.00
0.75
GMM
SW
0.50
0.25
0.00
'1
10
10
10
10a
Figure 25: Sliced Wasserstein (sW) distance as a function of n and for increasing noise levels e. Results are shown for both Gaussian toy examples with m = 10. Mean and std over different seeds.
11 ==
n=
n
n 30
1.5
SLCP
SW
0.5
Lotka-Vop Volterra
0.04
SW
0.02
0.00
DET. DET_ GEF
0.008
GAUSS -lp
SIR
SW
0.006
JAC-dip
0.004
LANGEVIN/.dI
0.002
0.000
Nazin
Nuain
Maxin
Nezin
Nazin
Figure 26: Sliced Wasserstein (sW) distance function of Nuaan and for increasing n, for the benchmark tasks SLCP, Lotka-Volterra, and SIR. Mean and std over 25 seeds.
Takeaway:
GAUSS is the most stable sampler, with the best results in highly noisy regimess or for learned scores (sbibm), even if its approximation appears to be less accurate than JAC DET GEF in the non- Gaussian (GMM) setting (for exact non-noisy scores).
JAC the most accurate method in the noie setting across both toy tasks, but remains unstable in practice. Clipping improves stability but introduces bias. A more principled stabilisation (beyond clipping) could yield further improvements and could be promising research direction.
DET GEF is fastest and competitive with our proposals, even outperforming GAUSS on GMM in the low-noise regime as n grows. However, still less stable than GAUSS as noise increases or in the learned score regime (sbibm), with similar behavior LANGEVIN.
51

---

Published in Transactions on Machine Learning Research (02/2026)
M.3 Clarifying Langevin behavior across settings
Noise-free toy regime. Figure 25 shows an interesting Irend: LANGEVIN improves n increases. This could be explained by the fact that, as n grows, the tall posterior becomes sharply concentrated and locally Gaussian (Bernstein-von Mises), making the sampling dynamics casier to handle with an appropriate stop size, whereas for small n the broader posterior would require more careful mixing (larger step sizes or more inner steps).
Noisy toy and learned-score regime. Figure 26 in contrast, shows that performance degrades with increasing (especially at large Miraia.. This is probably because small score approximation errors accumu- late across composed observations and are amplified by posterior contraction, leading to biased concentrated estimates and divergence. The fact that DET DET_GEE GEF exhibits similar instability indicates that this behavior is not merely a Langevin tuning issue, but rather an intrinsic sensitivity of samplers whose updates rely di- rectly on the composed learned scores, making them highly vulnerable to accumulated bias in concentrated regimes. In contrast, GAUSS and JAC approximate the tall reverse diffusion dynamics more explicitly through diffusion-coosistent backward-kernel constructions and associated covariance structure, which provides ad- ditional stabilization and yields more robust behavior under imperfect learned scores.
52


Reviews:
[{"review_id": "2H86ELbRdK", "reviewer": "Reviewer_p3EZ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper addresses the challenge SBI with \"tall data\", where the goal is to infer parameters given multiple independent observations. The authors comment on the limitations on the LAngevin sampling algorithm used by F-NPSE, and propose an alternative sampling algorithm that does not rely on Langevin dynamics (which at the same time involves selecting a few parameters, like number of steps per noise level, and step size). They derive this approximate sampling algorithm by approximating certain intractable terms in the \"true tall-data\" score with Gaussian distributions.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The authors provide an empirical evaluation showing that the method they propose, which does not rely on Langevin dynamics, yields good performance. While the paper has good evidence, I have some questions / concerns that I'd like to have addressed.", "audience_yn": "Yes", "audience_explanation": "I think people working in the field of SBI, specifically with problems involving large sets of observations, could find this work useful and interesting.", "requested_changes": "A few comments / questions / concerns.\n\n**Comparison with Deterministic Baselines**\n\nI think it would be interesting if the authors compare against the non-Langevin sampler from Geffner et al. (2023), presented in their appendix. That sampling algorithm also circumvents Langevin dynamics, and appears to do it in a way that's not completely foreign to the directions pursued in this work. It appears that both works use Gaussian assumptions to render the sampling tractable, though they appear to apply these approximations for different terms: the proposed method in this work approximates the backward kernel ($q_{0|t}$), whereas Geffner et al. approximate the backward transitions $q(t-1 | t)$. I think it would be valuable to compare these approaches, from both an empirical perspective (i.e. include it in the empirical evaluation), and a theoretical perspective, trying to shed light on how the different approximations used differ and whether they may be equivalent in some toy scenarios maybe?\n\n**Asymptotic Exactness and Structural Bias**\n\nA strong theoretical property of the Langevin sampler is its asymptotic exactness; with a step size tending to zero and a sufficient number of steps, it draws from the exact distribution defined by the score field. The method proposed here, however, relies on two distinct approximations: first, the approximation of the intractable backward kernels as Gaussians, and second, the approximation of the covariance matrices as constant (in the GAUSS algorithm) to avoid Jacobian computations. Consequently, regardless of the sampling budget or the precision of the score network, this method cannot theoretically recover the exact posterior structure for non-Gaussian distributions. While the authors correctly identify that the asymptotic behavior of Langevin is difficult to realize in practice (and I fully agree with this), the proposed method trades this for an irreducible structural bias. Am I missing something about the proposed method, regarding whether it can be made exact with sufficient sampling compute?\n\n**Analysis of Baseline Performance and Hyperparameter Tuning**\n\nIn Figure 3, the Langevin baseline performs worse as the amount of training data ($N_{train}$) increases. I find this counter-intuitive, as providing more data to the score network should result in a more accurate approximation of the gradients. Can the authors comment on why the performance gets worse with more training data?\n\n**Approximating the Intractable Correction Term**\n\nThe derivation of the tall data score hinges on the \"correction term\" $L_\\lambda$, which involves an intractable integral over the product of backward kernels. The authors currently address this by forcing the kernels to be Gaussian so the integral becomes analytically tractable. However, this imposes the structural bias mentioned previously. I would be interested to hear the authors' thoughts on whether this term could instead be approximated using a different method (e.g. Normalizing Flows, or a DeepSet-style network). If one could train a network to approximate this term directly (conditioned on the set of single-observation scores), it might be possible to avoid the Gaussian assumption entirely. Do you think this is possible? I understand that non Gaussian approximations will lead to issues when computing the resulting integrals, so this may not be possible. Just wondering if you had any thoughts on this.", "broader_impact_concerns": "No.", "additional_comments": "One more question. How much of these results do you think comes from poor tuning of the Langevin sampler? Having to tune it is, of course, a weakness of the method. That being said, if tuned properly, it should work better than the proposed methods, or is that not correct? I am mostly thinking about this from the perspective of the bias introduced by the approximations, which seems to be irreducible.", "extra_scores": {}, "extra_text": {}}, {"review_id": "ihziFqLr5l", "reviewer": "Reviewer_Tubi", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents an SBI method for tall data settings based on the scores of single-observation posteriors. This is made possible by an exact decomposition of the score function for multi-observation posteriors. The per-observation scores can be estimated purely from joint samples, which are available in the SBI setting. While the introduced decomposition of the score is not fully tractable, the authors present two different ways to approximate it. This gives (approximate) scores for a known diffusion process, allowing deterministic samplers such as DDIM to be applied.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The paper is well written overall, with a clear motivation by pointing out issues with F-NPSE. The exact decomposition of the tall data posterior score seems correct. Several approximations made to then estimate this seem reasonable, some questions on this below. The results show that this method cab work in principle and often better than F-NPSE.\n\n### Questions\n\n- The correction term in the tall data score, $L_\\lambda$, is approximated using the Tweedie framework. Can the authors elaborate on what kinds of assumptions this rests on and under which circumstances this is a good or bad approximation?\n- Why does LANGEVIN improve when conditioned on taller data? See e.g. Figure 2, Gaussian.\n- For problems in 4.1, does the method not learn the score but just computes it analytically except for $L_\\lambda$ and the difficulty is then to handle the noise?\n- Why not compare methods with the same time stepping? The \"equivalent time\" setting is interesting, but seeing methods with the same number of diffusion steps would also be relevant.\n- Why does Langevin performance worsen with increased $N_{train}$ (Figure 3,  LV and SIR)?\n- For GAUSS, is the time to estimate the $\\Sigma_{0,j}$'s included in Tables 1 and 2?\n- Can you better explain what you are doing with the clipping? The explanation in 4.2 is a bit limited.", "audience_yn": "Yes", "audience_explanation": "SBI and SBI in tall data settings is a relevant problem in many scientific domains. It is clearly not solved yet and this work makes an interesting contribution to the field.", "requested_changes": "I will wait for answers to my questions before I request changes.", "broader_impact_concerns": "This is fundamental research with no immediate impact concerns.", "additional_comments": "While I have a good understanding of SBI, I am not an expert in diffusion.", "extra_scores": {}, "extra_text": {}}, {"review_id": "mrboF7UV5y", "reviewer": "Reviewer_So1v", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Authors consider a problem of identifying parameters of the system $\\theta$ from multiple observations $x_1,\\dots, x_n$ in form of a posterior distribution $p(\\theta|x_1,\\dots, x_n)$. It is assumed that: (i)$p(x|\\theta)$ is intractable; (ii) it is possible/easy to generate $x$ from $\\theta$ by simulation. The approach proposed by authors is to:\n1. approximate $p(\\theta|x)$ by score matching\n2. approximate $p(\\theta|x_1,\\dots, x_n)$ using prior distribution $p(x)$, $p(\\theta|x_i), i=1,\\dots,n$ and data generated from $p(\\theta|x)$\n\nAs a result, sampling from approximated $p(\\theta|x_1,\\dots, x_n)$ can be done by backward diffusion process.\n\nAuthors demonstrate that their approach is cheaper and more accurate than Factorised Neural Posterior Score Estimation, a method that extends Neural Posterior Score Estimation to the tall data setting.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Authors made several claims:\n1. Introduction of a novel method that is similar to F-NPSE requires only learning $p(\\theta|x)$ but that works in a tall data setting\n2. Introduced method does not require Langevin dynamics\n3. Empirical results showing that new method is superior to F-NPSE (better accuracy and stability)\n4. Evaluation is done on multiple standard simulation-based inference benchmarks including challenging real-world problems.\n\nIn my view, the main claims are well-supported by the results of the authors.", "audience_yn": "Yes", "audience_explanation": "Given that the work is on the intersection of generative modelling, AI4Science and numerical methods, many readers of TMLR can be interested in the contribution.", "requested_changes": "The paper is well-written, with clearly stated motivation and detailed discussion of related work. I have only a few minor questions that appear below:\n1. The derivation of the method (Sections 3.1 and 3.2) is motivated by the assumption that the approximation of $p(\\theta|x_1,\\dots, x_n)$ should be suitable for sampling with a backward diffusion process. Is it possible to say something about the accuracy of this approximation? Possibly, given some convenient assumptions.\n2. Lemma 3.1 contains an assumption that a certain matrix is positive definite. In appendix A.1 authors provide some discussion when this is the case. This condition is possible to check numerically. Did the authors try that?\n3. Algorithm 2 requires an approximate covariance matrix. One may expect that with large number of parameters such an approximation will become less reliable. Can the authors please comment on that?\n4. In both Algorithm 1 and Algorithm 2 dense matrices are formed and stored. Does it mean that the methods are not applicable to parameter spaces of very large dimensions, e.g., $m = 10^{6}$?\n5. I suggest clarifying the meaning of $m$, $d$ in Section 4.2.\n6. In Appendix J. 2. there are references to undefined Section and equations. This appendix also mentions \"this thesis’ introduction\" which is likely a typo.\n7. Figure 5 is too small. Can the authors provide a larger version of this paper in the Appendix or in the main text of the paper?\n8. Similarly, Figure 10 in the appendix is too small. Is it possible to provide a larger version of it?\n9. In Section 4.3. authors refer to $g$ as \"gain\". I suggest removing that or to provide a more detailed discussion of parameters in the main text.", "broader_impact_concerns": "na", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)