A Statistical Benchmark for Diffusion-Posterior-Sampling Algorithms
Authors: Martin Zach
Youssef Haouchat
Michael Unser
Venue: ICLR 2026 Poster (2026)
Keywords: Diffusion models
Bayesian inverse problems
statistical evaluation
Gibbs sampling
Primary area: datasets and benchmarks
TL;DR: We made an evaluation pipeline for diffusion posterior sampling algorithms for Bayesian linear inverse problems that relies on the construction of posteriors with known posteriors that we can efficiently sample from.

Abstract: We propose a statistical benchmark for diffusion-posterior-sampling (DPS) algorithms in linear inverse problems.
Our test signals are discretized Lévy processes whose posteriors admit efficient Gibbs methods.
These Gibbs methods provide gold-standard posterior samples for direct, distribution-level comparisons with DPS algorithms.
They can also sample the denoising posteriors in the reverse diffusion, which enables the arbitrary-precision Monte Carlo estimation of various objects that may be needed in the DPS algorithms, such as the expectation or the covariance of the denoising posteriors.
In turn, this can be used to isolate algorithmic errors from the errors due to learned components.
We instantiate the benchmark with the minimum-mean-squared-error optimality gap and posterior-coverage tests and evaluate popular algorithms on the inverse problems of denoising, deconvolution, imputation, and reconstruction from partial Fourier measurements.
We release the benchmark code at https://github.com/zacmar/dps-benchmark and invite the community to contribute and report results.

Published as a conference paper at ICLR 2026
A STATISTICAL BENCHMARK FOR DIFFUSION- POSTERIOR-SAMPLING ALGORITTMM
Martin Zach!? Youssef Haouchat', and Michael Unser!
'Biomedical Imaging Group, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland 2CIBM Center for Biomedical Imaging, Lausanne, Switzerland fmartin. zach, youssef. haouchat, michael. unser) (ense.)
ABSTRACT
We propose a statistical benchmark for diffusion- diffusion-posterion posterior-sampling (DPS) algo- rithms in linear inverse problems. Our test signals are discretized Lévy processes whose posteriors admit efficient Gibbs methods. These Gibbs methods provide gold-standard posterior samples for direct, distribution algorithms. They can also sample the denoising posteriors in the reverse diffusion, Monte Carlo estimation of various objects
comparisons with DPS
which enables the Ebitrary-precisiicn that may be needed in the DPS algorithms, such as the expectation or the covariance of the denoising posteriors. In turn, this can be used to isolate algorithmic errors from the errors due to learned components. We instantiate the benchmark with the cewiimun-mma---quureeeeee. evaluate popular algorithms the inverse problems of denoising, deconvolution, imputation, and reconstruction from partial Fourier measurements, We release the benchmark code ht tttps: / //gittull com/ za cmar /emar/dps--... waeaar/dds--echaar... benchmark and invite the community to contribute and report results.
optimality gap and posterior-coverage tests and
INTRODUCTION
Diffusion models are among the leading generative models in imaging (Rombach al., 2022), visual computing (Po et al., 2024), finance and time-series analysis (Huang et al., 2024; Rasul et al., 2021), I), de novo protein and drug design (Watson et al., 2023; Alakhdar al., 2024), natural language process- ing (Li et al., 2022), and other domains. Their ability to model complex distributions has motivated their use as priors in the Bayesian resolution of inverse problems. In fact, reconstruction methods that leverage diffusion models are competitive or state-of state-of-the-art for problems such as deconvolution (Ren al., 2023), magnetic resonance imaging and computed tomography reconstruction (Chung Ye, 2022; Liu et al., 2023), weather-artifact removal (Özdenizci & Legenstein, 2023), task-conditioned protein design (Bogensperger et al., 2025), audio bandwidth extension and dereverberation (Lemercier et al., 2024), and denoising of financial time-series (Wang & Ventre, 2024).
This empirical success has come in spite of a lack of a natural mechanism for the conditioning on measurements and active research explores how to incorporate the likelihood (Yismaw al., 2025; Erbach et al., 2025). Currently, conditioning strategies are evaluated in one two ways. (i) With respect to downstream applications: As an example, evaluations with respect perceptual metrics such as the structural similarity (Wang et al., 2004), the Fréchet inception distance (Heusel et al., 2017), or the learned perceptual image-patch similarity (Zhang et al., 2018) are common in the imaging sciences. As pointed out by Pierret & Galerne (2025b) and Cardoso et al. (2024), however, these metrics are ill-suited for the statistical evaluation of pesteriorsamplinn overly simplistic settings: A common fallback is to evaluate conditioning strategies in synthetic settings with (finite-com componeem))) the tail decreasing exponentially like the widest component. Consequently, they cannot reproduce extremes that are common in asset returns (Blattberg & Gonedes, 1974; Cont, 2001) or statistics of images (Wainwright& & Simoncelli, 1999). We illustrate signals with such power-law-like extremes later in Figure 3. Benchmarks built on such priors can therefore overstate posterior quality. A proper statistical evaluation in realistic settings is critical in high-stakee applications such as medical imaging, remote sensing, and finance, where decisions based on reconstructions and their associated uncertainties may have significant consequences.
algorithms. (ii) In
caaussian-mixtue priors. Such mixtures remain light-tailed with

---

Published as a conference paper at ICLR 2026
1.1 CONTRIBUTIONS
We propose such a statistical benchmark for linear inverse problems. Our test signals are discretized Lévy processes that admit efficient posterior- sampling algorithms. Indeed, they admit efficient Gibbs methods with exact conditionals that provide gold-standard posterior samples. Our framework supports general posterior-leve comparisons (e.g., (sliced) Wasserstein energy distances calibration via coverage and posterior predictive checks) by furnishing matched samples obtained from the DPS algorithms and the Gibbs methods.
(DPS) algorithms' for
The Gibbs methods also suited sample from the denoising posteriors in the reverse diffusion. This motivates our contribution of a new template for DPS algorithms, in which update steps utilize samples from the corresponding denoising posterior. These samples can be used for arbitrary-precisst Monte Carlo estimation of various objects that are needed in the update steps of the algorithms, such as the minimum-meane mneen-squared-enn of algorithmic errors from approximation errors due to learned components. We show how several popular DPS algorithms can be re-expressed within our template.
(MMSE) denoiser or its Jacobian, which enables the isolation
Finally, we instantiate the framework with the MMSE optimality gap and coverage checks across the inverse problems of denoising, deconvolution, imputation, and reconstruc- tion from partial Fourier measurements. We target the realistic scenario where a learned denoiser is used and check hyperparameter sensitivity by substituting the arbitrary-precision Monte Carlo counter- parts for the learned components. The benchmark code-which is another substantial contribution- available online. contains efficient implementations of sampling routines and a containerized run- time that allows novel algorithms to be easily benchmarked.
1.2 RELATED WORK
For unconditional sampling, many works derive theoretical bounds on distances between target distribution and the distribution obtained by (approximations of) the reverse stochastic differential equation (SDE) (see Section 2). For example, Gao et al. (2025) bound the Wasserstein-2 distance with respect to the discretization error of the SDE under the assumption that the target distribution is smooth and log-concave. This directly bounds the number of diffusii steps needed obtain a desired accuracy. Under absolute continuity of the target with respect to a Gaussian, Strasman et al. (2025) bound the Kullback-Leibler divergence with respect to properties of the noise schedule.
A common assumption that simplifies the analysis and facilitates the computation of various errors and bounds is that of Gaussian target. For example, Hurault et al. (2025) analyze the error incurred when using a finite number of prior samples for the estimation of the prior score and track its propagation through the iterations of the reverse-SDE solver. Pierret & Galerne (2025b) derive explicit solutions to the SDE and use them to derive bounds on the Wasserstein-2 distance to the distributions that are obtained via Malle--aruyyama discretizations.
For conditional sampling, Pierret & Galerne (2025a) derive expressions for the Wasserstein-2 distances between the conditional forward marginals and the distributions induced by specific likelihood approximations in the reverse SDE under the assumption of a Gaussian prior. Crafts & Villa (2025) systematically evaluate DPS algorithms numerically under the assumption of a comite-componet) Gaussian-mixe prior and provide reference objects to the DPS algorithms to ensure a fair evaluation. Cardoso et al. (2024) and Boys et al. (2024) also evaluate their algorithms on Gaussian-mixtut priors. These Gaussian-mixture mixture priors, however, cannot reproduce power-law- lawer-lawllike extremes and can overstate posterior quality.
Beyond diffusion-ssecifee theory, Thong et al. (2024) check the coverage of credible regions produced by different Bayesian recovery strategies and find that those that utilize diffusion models often under-report uncertainty. shorteoming of their approach is that they an empirical distribution of images as a surrogate for the prior distribution. Finally, Bohra et al. (2023) also used efficient Gibbs methods to obtain gold-standard posterior samples. Their main focus was to quantify the quality of neural MMSE estimators with different number of parameters. Our work extends this to posterior-level comparisons.
We use "DPS algorithms" as umbrella term for posterior-sampling methods with diffusion priors. The method due to Chung al. (2023), often called DPS in the literature, will he referred to later as C-DPS.
2

---

Published as a conference paper at ICLR 2026
2 BACKGROUND
Bayesian Linear Inverse Problems
We seek to estimate a signal x E IRd from the measurements
Ax + n,
(I)
where the forward operator A E R models the noiseless linear-measurement acquisition and 12th is additive noise. In the Bayesian resolution of this problem (Stuart, 2010), the signals are modeled as a random variable, denoted X, with values in R and distribution px. referred to as the prior. Given any measurement y. the ultimate goal is to analyze the posterior Px It is related to the likelihood Xx and the prior px via Bayes' rule, which states that
(x) & P (x)xxxx)))xx))
In contrast to classical variational methods (Scherzer et al., 2008), the posterior distribution provides natural means to quantify uncertainty and can be summarized by various point estimators. We provide a precise description of point estimators that are relevant in this work in Appendix A.
For given signal x, the likelihood py x is fully specified by the distribution of the noise. common assumption on the noise is that it is vector of independent and identically distributed (i.id.) Gaussian random variables with mean zero and variance on In this case, the likelihood given by
Pr(x=x(y) x cxp(- 1 Ax max-yln yl 20
Thus, once the forward model and the noise distribution are specified, the remaining modeling choice is the prior. Diffusion models are good candidates due to their ability to encode complex distributions.
Diffusion Models Diffusion models were introduced by Song et al. (2021) by unifying the discrete approaches proposed by Song & Ermon (2019) and Ho et al. (2020) in continuous theory based on SDEs (Klenke, 2020, Chapters 25 and 26). We denote the (diffusion) SDE with drift coefficiens
f Rdd Rzo + Rd and diffusion coefficient g Rg0 + R as
dx.- f(X,, t) dt + dW,,
where W, is the standard Wiener process. In our setup, the initial condition X, is the random variable that describes the signal, thus, Xo = X. Under suitable choices for and g. the forward process admits limiting marginal X. as / + co. Sampling from pxo can then proceed by simulating the
SDE (4) in reverse with initial condition X By Anderson's 's theorem (Anderson, 1982), the reverse SDE that reproduces the forward marginals satisfies
dX, (I(X..t) (()( log (X.)) di + g())dW,,
(5)
where px, denotes the density of X, defined by the forward process, and dt is negative.
The primary challenge in this approach lies in the computation of the scores V log Px, for all 0. A fundamental relation known as Tweedie's 's formula connects the score with the MMSE denoiser: As we derive rigorously in Appendix B, for f(x, t) (- 2 ] and g = /3(t), we have that'
logxx,
V loggx, (x) = -x a()) X, x]).
(6)
where a(l) = exp(- 0o (s)dss) and o² o²) = (1)). This yields practical way compute Px, (x) through the resolution of the MMSE denoising problem of finding E[X, X, = x]. In standard applications where the goal is the generation of new signals, one typically tackles this by approximating the map (x, f) + X X, = x) with neural network that is learned in an offline step. In our framework, we can instead obtain rrbitrary-precision MMSE denoisers via Gibbs methods and thereby eliminate approximation errors from learned surrogate and isolate errors in DPS algorithms themselves.
V loggxx,
The implementation of the reverse SDE for generation requires its own time discretization, for instance with Euler MareraMaruyama techniques (Higham, 2001). In this work, we will base our backward processes on the alternative denoising diffusson probabilistic model (DDPM) backward process (starting from Gauss(0. (0, I)) X --11 (X, + 3,V log px. (X.) + BB,Z..
(7)
VI-B
3

---

Published as a conference paper at ICLR 2026
999
600
200
0
1L
50
S((1)
10
20
60
20
40
60
0
20
40
60
20
60
Figure 1: Unconditional reverse-diffusion trajectories obtained by DDPM using the arbitrary-precisn Monte Carlo denoiser. Rows: Increment distributions. Columns: Diffusion times. Line styles: Differ- ent random states.
that originates from the discrete-time Markov chain that was initially proposed by Sohl-Dickstein al. (2015) and revisited and popularized by Ho et al. (2020). We relate it to the MareraMa discretization of the reverse SDE through Taylor expansions in Appendix B.I.
Though we defer details on our signals and the Gibbs methods that we to obtain the cubitrary-precision MMSE denoiser to Section 3, demonstrate in Figure that our signals can be generated by coupling the unconditional back- ward process in (7) with this denoiser. We fur- ther motivate this arbitrary-precision denoiser Figure 2 by comparing histograms of signal increments produced by the learned denoiser and the arbitrary- preeaari-rreiioon denoiser for a St(1) increment target (notations summarized in Ap- pendix C.2). The increments of signals gener- ated with the precciar denoiser follow the target almost perfectly. Residual errors are due to the discretization the reverse diffusion and Monte Carlo error of the rrbitrary-precisiie denoiser.
I Leamed
0.3
Monte Carlo
Frequency
Target
0.2
0.1
-10
0
10
Increments
Figure 2: Histogram of increments of signals ob-
tained by DDPM with different denoisers.
Diffusion Posterior Sampling Our reverse-diffusion sampler be adapted to sample from posterior by replacing the prior score V log with the posterior score
V P..P..
V log + (x + logpy(x.=x
for some given measurement y. obtained by Bayes' theorem. Although the dependence between Y and Xn is known through (1) and the likelihood is explicitly modeled via (3), is generally challenging to relate Y and X, for any > 0. To overcome this, the conditioning the measurements is usually done in one of two ways. (i) A learned component models the conditional posterior score and also gets the measurements as input. This strategy (pursued by Liu et al. (2023); Özdenizei & Legenstein (2023); Bogensperger et al. (2025); Saharia et al. (2023)) is ad advantagcous when the measurement process unknown, difficult model, or prohibitively expensive to evaluate. However, its reconstructions typically degrade under shifts in measurement conditions since the learned components cannot adapt to the new measurement conditions. and the likelihood score is approximated. This strategy (pursued by Chung & Ye (2022); Xue et al. (2025) and reviewed by Lemercier et al. (2024)) advantageous when the measurement process is known, relatively inexpensive evaluate, and subject to change, but prior knowledge should be reused, which is frequently the case in, e.g., imaging or remote- menote-sensing applications. However, this requires approximations to the likelihood score V(x H log Pogxxx(=)()) for all t O.
(ii) The Bayesian separation that is described (8) is pursued
Our framework supports more general (possibly non- Gaussaassii) likelihoods, see Section 3. This is the veriance-preservng formulation (Song et al., 2021, Section 3.4) with standard normal limiting marginal, where 8 : R20 R2zo controls the speed of the contraction to and how much noise is injected.
4

---

Published as a conference paper at ICLR 2026
Gaussian
Studeni
Laplace
AAP. log-price
5.61
Pr([X] > A) P 1
10-1
5.59
10
09:00
10:00
11:00
12:00
13:00
14:00 15:00 16:00
Time day on 2025-11-21
Return
Incenitty
10-3
10-6
200
400
0.2
0.4
0.6
0.8
Difference
Figure 3: Examples of signals with heavy-tailed increment distributions. Top: Asset returns. Bottom:
Columns in natural images. Left: Signals. Right: Survival function of absolute increments (no marker: empirical; markers: best fit to empirical within distribution).
Our benchmark can evaluate either strategy, as well any other method that would claim to sample from a posterior distribution like in (2), Approach (i), however, relies on black-box learning of the conditional posterior score and its performance heavily depends on various implementation details. Thus, we primarily focus on approach (ii. which necessitates approximations score (and more general DPS algorithms with explicit conditioning, see our proposed generalization in Section 3). For those, our framework can supply arbitrary-precision Monte Carlo estimates of various objects to isolate and quantify the impact of these approximations.
the likelihood
PROPOSED FRAMEWORK
The prior distributions in our framework will be that of signals of length obtained by regularly spaced samples of processes with independent, stationary increments (Lévy processes, described in Appendix C). Let s be such process and let the unit-step increments be [u]x = (s(k) s(k 1)) for k = 1, 1,., d. Independence and stationarity imply that the distribution of [u]a, the increment distribution pu. does not depend on k. The increment vector is related to the signal x via u Dx. where D is a finite-differenee matrix with an initial condition that allows us to write x = D where 1 a lower-triangular matrix of ones. This gives a convenient way to synthesize signals once u is drawn. The independence of the increments implies that the density of the discrete signal is
(x) =
II pu (Dx]a).
We consider four increment distributions that are commonly used in sparse-process models: Gaussian; Laplace; Student-t; and Bernoulli-Laplace (spike-and-slab), Such increment distributions are sparse or heavy-tailed according the taxonomy due Unser & Tafti (2014) and are relevant in signal and image processing, finance, and other fields (Schoutens, 2003). We show instances signals with such heavy-tailed increment distributions in Figure 3. A precise definition of Lévy processes, the matrix D, the increment distributions and their notation along with a discussion about extensions to higher-dine digger-dimensional signals signals with more complicated graph structure given in Appendix C.
With the prior distribution specified in (9) and the assumption of Gaussian noise, the posterior associated to the inverse problem intrinsic in (1) is
Efficient Posterior Sampling
PXIY=y(X) x exp(expp Z07 |Ax y1²))xx(x exp( 20
Ax y 11)
II pu (Dx)):). (10) k=1
Unless pu is Gaussian (the simplified setting in Pierret & Galerne (2025b)). this posterior is not conjugate, so neither closed-form sampling nor direct evaluation of moments is available. Nevertheless, for the increment distributions used in this paper, the posterior distributions admit efficient Gibbs methods via standard ditent-variable augmentations. Motivation and more details about the Gibbs methods, such as the burn-in period and the number of samples S, are provided in Appendix D.
5

---

Published as a conference paper at ICLR 2026
The Gaussian, Laplace, and Student-t dis- tributioos admit latent representations as infnitt-component Gaussian mixtures, which makes them suitable for the Gaus- sian latent machine (GLM) that was cently introduced by Kuric et al. (2025).
Algorithm 1 GLM Gibbs method.
Require: Xo R., K E Raxd conditional la-
tent distributions {Piz)Ii k== and maps
{UA: [
1: for 8 = I...., , B do
It is generally applicable distributions
2: Draw [z]8 P= par. over k
11
3: Draw X, Gauss( E(x)
p(x) x II Φn (Kx))), k=1
(11)
4: return (x++s)s=1
18
where E guxd and all distributions Φ1 $2, On IR
have a latent representation
(xx(t) (p.(s)) (()) (e (z),
(12)
where the latent distribution f& and the latent maps Pk. of and 0u.a² is the density of a dieeniinnsional Gaussian distribution with mean μ and variance o². We can cast the posterior distribution in (10) into this framework by rewriting it as
depend on the distribution di-
122
m+d
PXIY=y (x)
II
givle o ([Ax]s)
II PU (Dx)a) k=1
II (Kx|&))
(13)
There, K = [A; [A: 0z Olyle.oz for = 1, Importantly, non-Gaussian likelihoods can be handled by some appropriate definition of the first m distributions.
..... ..... and dR py for l: m++ m ................................ -n+dd
The introduction of an appropriate redimennional random variable with nontrivial distribution (see the details in Kuric et al. (2025)) enables the efficient sampling from the conditionals: Sampling X Z z amounts to sampling a Gaussian with covariance and mean
Ez) (K (o) and (z() = E(z) Ez) (o)))
(14)
(wI(zz)),.... Men (lz),)). Sampling Z X = x amounts to sampling independent one-dimensional conditional latent distributions p -= P'zj.,.................... X="K),, that depend on the distributions Φ1...., On and are given in Table 3 in the appendix along with the corresponding latent distributions and latent maps. We summarize the GLM sampling in Algorithm 1.
respectively. where (z) diag( ( ((.))...... ([z),)) and no
of
For the Bernouli--aplace increment distribution, adapt the algorithm proposed by Bohra et al. (2023) that introduces two d-dimensional latent variables: a Bernoulli indicator ("on""'"of""); and haplace-distribute increment height. For self-contained exposition, we rigorously derive the resulting Gibbs method in Appendix D.1.
The Gibbs methods that just described are suitable for the generation of the gold-standard samples from the posterior that corresponds to
Seconds per Gibbs iteration (NVIDIA V100)
Baseline
the initial inverse problem intrinsic in (1) as well as the generation of samples from the de- noising posteriors in the DPS algorithms. In the latter case, the forward operator A is the iden- tity, the measurements are the noisy intermedi- ate reconstructions X,, and the noise variance og = of follows the schedule at timestep I.
Panaltel cuDA GIG sumpler
51.07
Purallel sumpling U
V W ww
-33.52
Woodbury-
sampling
W w.
11.7|
1.36s Speedup: 74.61
+ storah.compp compri le
When these Gibbs methods are embedded within the reverse-diffusion loop, efficient implementation is paramount to achieve ac- ceptable runtimes. This is most crucial for the Bernoulli Laplace increment distribution, where the sequential drawing of the binary support vector is nested within the Gibbs loop, which in turn may be nested within the deverse-diffusie loop. Accord- ingly, we tailored our implementation to modern, highly parallel compute units and optimized several components, including custom CUDA- and Tri ton-copiieed sampling routines and incremental updates based on the Woodbury-Sherman www-eeeaee-rrrriieeeeeeeee identities We achieved a cumulative speedup of 74.61 x over the baseline implementation (illustrated in Figure 4 with details in Appendix D.2).
Figure Runtime improvements of the Bernoulli- Laplace sampler.

---

Published as a conference paper at ICLR 2026
A Generalized DPS Template Widely used methods, such as diffusion plug-and-play (DPnP) (Xu & Chi, 2024), fall outside the pat- tern described in Section 2, where one approx-
Algorithm 2 Template for DPS algorithms.
Require: Initial point XT: y, A, A 1: for t= ....., ., do
D Diffusion process
imates likelihood score inside the reverse dif-
2:
Sample {xx 1 Xxx
3: Update X1-1
fusion. We therefore introduce a simple tem- plate that is natural in our setting and accom- modates a broader set of DPS algorithms. More precisely, we characterize the iteration rule of
S(x,, {*+ K- A,1,
4: return xalg Xo
> Posterior sample
DPS algorithms as a two-stage process: Given an iterate x, with associated noise variance of. the next iterate x 1 is computed by (i) drawing S samples denoted [xx] from the denoising posterior Px, O exp(-
1 20? xxll²) px, and (ii) the subsequent computation of an update step
S that may utilize the current iterate Xt, the samples {xx J the measurements y, the forward operator A, and, possibly, other algorithm-intere parameters such as a scalar that weights likeli- hood and prior terms or parameters that define the noise schedule. This template is summarized in Algorithm 2 and specialized instances for the update step S for a variety of popular algorithms are given in Appendix E.2. We have absorbed the (variance-preserving scaling into the step S since this template is not fundamentally limited to diffusion processes but supports any (also not monotonically decreasing) noise schedules. In addition. noise variances {ot]T= are usually derived from the internal parameters À that may include noise schedule.
Through this construction, DPS algorithms can use any statistic R of the samples {xx] in their update steps. Most methods use the mean R(x Xs) = = E X == jt, which is the Monte Carlo estimate of E X, = xz] An example a DPS algorithm that utilizes additional statisties Chung diffusion posterior sampling (C-DPS), which requires the Jacobian of x, E[Xo |X, x,]. As we show in Appendix 1. this Jacobian equals (up to the known variance- variance-preserving scaling) the conditional covariance of X, X, X,, unbiased estimator of which can be obtained through the statistic R(X,, xs) = I-1 I (*; - p An example of DPS algorithm that utilizes an alternative statistic is the DPnP algorithm that alternately samples from PX,(X,=x and data-proximal problem. There, (..,,...,, &g) = x, is used to obtain one sample from Px, X2 =x This statistic is frequently used in the asymptotically exact and the CSGM-type algorithms (using the taxonomy due to Daras al. (2024)). When only learned MMSE denoiser is available, obtaining this one sample requires a full reverse diffusion. In contrast, it requires only one iteration (and the burn-in period) with the Gibbs methods. Thus, these algorithms are typically faster when they are endowed with the Gibbs methods (see the runtimes in Tables 5 and 6), which enables easy benchmarking. However, CSGM-type algorithms typically do not aim at posterior sampling and we do not benchmark them here.
Since the denoising posteriors are always sub-Gaussian, the Monte Carlo estimation of any object enjoys favorable convergence. For instance, the computational complexity of estimating the covariance up to a desired precision the operator norm scales linearly with the dimensionality of the signal (Vershynin, 2018, Theorem 4.7.1).
Extensions A prerequisite for a quantitative evaluation of posterior-sampling algorithms is the availability of reasonably efficient samplers that can provide gold-standard samples. The development of such samplers for posteriors arising from nonlinear measurement models and non-Gaussian noise is challenging, and existing methods currently address only specific cases (e.g., Wang et al. (2017) study a nonlinear Nonlinear-Gaussian measurement model with a Laplace prior). Importantly, our framework modular: as more general- purpose samplers for these posteriors become available, they can be plugged into our benchmark directly. The denoising posteriors in the reverse diffusion do not change with the likelihood and can, therefore, always be efficiently sampled.
When going to higher dimensions, the primary challenge lies in the sampling of the high-dimensional Gaussian distributions required in the Gibbs methods. Luckily, the structure of the involved operators in our case is such that the Gaussians can be efficiently sampled with perturb-and-MAP approaches with matrix-free conjugate gradient implementations; we discuss this in more detail and show how the runtime of different samplers change with the dimensions in Appendix D.2. Sampling high- dimensional Gaussians is a well-studied problem and advances in that field can directly be used in framework.
7

---

Published as a conference paper at ICLR 2026
Our gold-standard posterior samples can be compared to samples obtained by any posterior-sampling algorithm. This includes classical Markov-chain Monte Carlo algorithms, algorithms that utilize flow-matching priors, and others. In this work, primarily focus DPS algorithms because our framework can supply arbitrary-precisonn Monte Carlo objects to them. We helieve that this fundamental principle can be extended to other algorithms, in particular those that utilize flow- matching priors. Such algorithms are frequently evaluated on toy examples based on Gaussian mixtures (e.g. by Pourya et al. (2025)), that are overly simplistic.
NUMERICAL EXPERIMENTS
We consider signals of dimension d 64 and four inverse problems that are frequently encountered in various estimation tasks throughout the natural sciences: denoising; deconvolution; imputation; and reconstruction from partial Fourier measurements. The dimension of the signal is large enough such that the corresponding operators can be sensibly defined, small enough such that the benchmark has acceptable runtimes. We provide experiments about the runtime with larger signals in Appendix D.2, details about the operators in Appendix F.1, and precise descriptions of the benchmarking pipeline (e.g., the number of training. validation, and test signals, and the number of iterations in the Gibbs methods) in Appendix F.2.
4.1 RECONSTRUCTION ALGORITHMS
Model-Based Methods
We consider the model-based methods
x (y, A) = arg min
- yyl2² + (2llAx-y11*
Dx|)))
(15)
xER2
and
x4 (y, arg min (111Ax- - D||Dx|111 xER
(16)
as baseline reconstruction algorithms. They coincide with the maximum- a --osstrrroo (MAP) estima- tors of Lévy processes associated with Gaussian and Laplace increment distributions, respectively.
Diffusion Posterior Sampling Algorithms We consider C-DPS (Chung et al., 2023), diffusion models for plug-and-play image restoration (DiffPIR), (Zhu et al., 2023) and DPnP (Xu & Chi, 2024). This selection demonstrates the applicability of the framework to algorithms that require denoising-posterior samples (DPnP), the MMSE denoiser (DiffPIR), and its Jacobian (C-DPS), which covers most of the existing DPS algorithms. For each DPS algorithm, we benchmark variant that uses learned components (learning details are provided Appendix F3) and a variant that uses Gibbs samples of the denoising posterior. For DPnP, this fully removes approximation errors. For the others, the learned components and the Monte Carlo estimates of those components have varying quality for different distributions and noise variances that we systematically investigate in Appendix F.4. We provide our main results, the MMSE optimality gap, for the learned variant and then investigate changes when substitute the Gibbs samples for the learned components.
The model-based methods and the DPS algorithms require the tuning of some hyperparameters. These were found by grid search on validation data independently for each algorithm, increment distribution, and forward operator. The precise setup for this grid search is given in Appendix F.5, The hyperparameters for the DPS algorithms were tuned to the learned denoiser. Parameters obtained with this procedure are later denoted with star in the superscript.
Gold-Standard Gibbs Methods The Gibbs methods are used to obtain gold-standard samples from the posterior. As described in Section 3, the Gibbs methods are parameter- and bias-free and efficient. Consequently, they are well-suited for our purpose. Chain lengths, diagnostics, and implementation details are given in Appendix F.2; reuse the same settings across operators and increment families.
4.2 RESULTS
For any measurement y. some DPS algorithm alg that depends on the parameters À will produce samples that we denote {x (y, A)) Namppes, We moreover denote x galg
MMSE (y, A)
8

---

Published as a conference paper at ICLR 2026
Table 1: MMSE optimality gap in decibel (mean + standard deviation; lower is better; 0 is a perfect reconstruction) of various estimation methods over the test set. Bold: best among DPS algorithms.
Gaues(0.0222 0.25)
Laplace[1]
BL(0.1.1)
St(1)
St:2)
Si(3)
C-DPS
0.12 +: 0.18
0.12 + 0.20
2.22 + 2,26
3.26 1.01
0.28 0.30
0.10 + 0.18
DinPIR
0.16 + 0:21
0.09 + 0.16
0.72 + 1.10
0.93 - 1.06
0.07
0.14
0.15 4 0:21
Denoising
DPaP
0.24 2 0,25
0.11 0.117
1,33 4: 2.12
1.19 4: 1,38
0.10 = 0.17
0.10 + 0.17
0.15 4: 0.21
0.06 0.112
3.41 4.338
0.38 4: 0.13
0.14 0.19
0.11
0.00
0.18 1.281
8.61 3.10
3.25
0.74
0.25
C-DPS
0.12 + 0.20
0.22
4.30 + 3.87
18.30 +: 6.288
0.46 1.40
0.17 + 0.53
DiAPIR
0.0T + 0.17
0.07 0.19
1.09
+: 2.22
10:45 6.10
0.09
0.57
0.08 + 0.26
Deconvolution
DPaP
0.10
0.13 0.22
1.71 2.49
7.84 5.66
0.35
0.14
L.65 d: 0.84
1.38 0.88
1.86 2: 3.14
L.87 4.01
1.10 =: 1.19
1.28 4: 0.94
0.00 + 0.01
0.07 + 0.23
6.11 + 4.40
21.50
4.46
1.44 + 2.85
0.36 + 1.00
C-DPS
111:::0.2
2.99 4:
0.50
0.57
DiAPIR
0.09 d: 0.23
0,08 1: 0.24
0.24 d: 1.14
0.88 :: 3.50
0.11 d: 0.62
0.08 1: 0.42
Imputation
DPaP
0.14 0.312
0.17 2 0.38
0.50 # 1.28
10.89
+: 5.92
0.25 =: 0.82
0.27 + 0158
1.74 + 1.12
1.77 + 1:35
1.25 + 2.78
13:32 5.32
1.37 2.56
1,53 1.58
0:00 0.01
0:01 + 0.05
1.10 + 1,88
0-42 +: 0.86
0.06 = 0.34
0.02 + 0.28
C-DPS
0.15 d: 0,36
0.26 4: 0.65
5.90 +: 4.41
4.29 d: 6.78
0.53 =: 0.83
0.35 à 0.77
DiNPIR
+ 0,29
0.08 +:
0.83 4: 1.44
3:19
+ 0.39
0.12 0.37
Fourier
DPaP
0.11 0.95
0.20 0.51
1.88 2.47
2.45
4.83
0.39
0.24
4
1.50 + 1.59
0.73 + 0.944
3.57 + 2.852
1.07 + 2.98
0.71 = 0.99
0.78 + 0.97
0,00 0.002
0.36 0.73
12.22 d: 1.53
9:17 8.54
2,66 3.57
1,08 1.79
N Neamples Naamllns X (y, A). For an estimation method xe (-) and data with corresponding
data-generating signal measure the MMSE optimality gap (in decibel) defined by
(y))
x
10 10
(17)
Gibbs MMSE (y) x|1
where x (y) = xxiva(y, A*) for model-based methods and x (y) sag
(y, X*) for DPS algorithms. gap of indicates perfect recovery of the gold- standard MMSE estimate and any positive values show the orders of magnitude of the error relative to the reference error. We found
that Nsamples = 50 provided good tradeoff between runtime and accuracy by benchmarking the gold-standard Gibbs method with that number of samples.
We report in Table the mean and standard deviation of the MMSE optimality gap over all signal- measurement pairs yx in the test set obtained by the model- methods and the DPS algorithms endowed with the learned denoiser. The Gaussian increment distribution validates the implementation: Since the MMSE and the MAP point estimates coincide, the model- &g estimator matches the Gibbs reference up to the error due to the finite parameter-gre grid resolution. When the posterior mean is smooth (e.g.. imputation and some deconvolution cases). & is the best model- choice and frequently outperforms the DPS algorithms. When the posterior mean is close to piecewise- constant (typical in denoising of signals with sparse increments), the & estimator is preferred. Among DPS algorithms, DiffPIR is typically the top performer and often exceeds &g and e baselines in deconvolution, imputation, and reconstruction from partial Fourier measurements. For spike-and-slab settings (Bernoulii-appaace) DPS algorithms substantially outperform the model-based baselines across operators. In deconvolution and reconstruction from partial Fourier measurements, DPS algorithms frequently match or surpass the best model-based estimator.
We now inspect the change in performance after we substitute higher-quality Monte Carlo com- ponents for the learned components. We do this without retuning of the hyperparameters, which allows us to see if the performance of the algorithms increases automatically with the quality of the denoiser. Here, we discuss general trends; exhaustive quantitative evaluation and precise quantification of the quality of the learned and Monte Carlo objects is given in Appendix G. For the same hyperparameters, the performance of DPnP increases significantly with the quality of the mmmossing-posterii samples. For example, the optimality gap decreases by 10.46 dB for imputation of signals with S1) increments, and significantly for other measurement operators for signals with Si(1) increments and (L.1, 1) increments. By contrast, C-DPS and DifiPIR can require retuning when the denoiser changes: Scores can deteriorate after one has substituted higher-quality Monte Carlo denoiser for the learned one, but a brief hand-tuning of the hyperparameters on the validation set improves them way beyond the learned denoiser. For instance, for DifiPIR and imputation of signals with St( increments. reusing the hyperparameters deteriorates the gap by 13.56 dB, whereas a brief hand hund-tuning decreased the optimality gap by almost 10 dB over what is reported in Table Qualita- tive examples of the MMSE estimates and the marginal variances obtained by the DPS algorithms and the gold standad Gibbs methods are shown in Figures 11 to 18 in the appendix.

---

Published as a conference paper at ICLR 2026
Prototypical samples and the corresponding MMSÉ estimate obtained from a DPS algorithm (here DiftPIR for deconvolution of signal with (L. 1) increments) are shown in Figure 5. (The full conditional reverse-diffusion trajectory, the data-generating signal, the measurements, and the MMSE estimated obtained with the gold- standard Gibbs methods are shown in Figure 19 in the appendix.) The figure highlights a key dis- tinction: Posterior samples often preserve high- frequency structure and reflect prior variability. whereas the MMSE point obtianeedbtained by averaging all samples is much smoother. This explains why DPS methods tend to score higher on perception-orienter metrics, while regressors that target the MMSE point estimate (through training with the mean squared error) excel on distortion metrics like the peak aggnal-o-noise-rre
& (y))
DiEPIR SXMME (V)
x x (y)
40
Figure 5: Conditional generation for deconvolution of a signal with (0.11,,1 1) increments with Diff- PIR. The shaded area indicates the variance.
(PSNR). Consistent with this distinction, Saharia al. (2023) fairly compare a sampling- -based method to an MMSE regressor and find the expected tradeoff: higher PSNR and structural similarity for the regressor; and better perceptual scores for the sampler. We therefore recommend to make the Bayesian target explicit point estimate versus sample quality-and and to use evaluation protocols that are aligned to that target. Our framework supports this by offering gold-standard posterior samples and urbitrary-precisionn Monte Carlo estimates.
In addition to the evaluation of the MMSE optimality gap we analyze the highess-posteriir-densii coverage of the algorithms. Specifically, for any measurement and any k 1,2...., Naumples we define* (y)
logpx Y=y ( (x P(k) (y, Xals:*) where P is the permutation that ensures that threshold at a E 1] as I 1 (y). We declare the gentragnnn signal x covered Px|yxXy= (x)) 2 4 Nammmepe (y) and define the coverage of a method as the fraction of
(y)) 2 (y) ... 2..2) 2N (y) and define the empirical highest-posterior-ee
pairs yx (y)) The coverage of a calibrated posterior-sampling method will be up to Monte Carlo error. A coverage result that is less than indicates that the samples concentrate too heavily around the mode; a coverage result that is greater than indicates that the samples are too spread out. We again discuss general trends here and present an exhaustive quantitative evaluation in Appendix G. The coverages obtained by the DPS algorithms are generally much smaller than a, which indicates that they are uncalibrated and is in line with what is reported by Thong et al. (2024). For C-DPS and DiffPIR, the reported coverage values are almost always 0 except for (0.1,, 1) and St(1) increments, where the coverages are usually (close to) 1 for C-DPS and inconsistent for DifIPIR. For almost all increment distributions and forward operators, DPnP reports coverage values that are closest to but typically smaller than a.
in the test set for which x is covered by the threshold I [aa.Namphee
5 CONCLUSION
We have introduced a statistical benchmark for diffusion posterior sampling algorithms for linear inverse problems. The framework constructs signals with a known distribution, simulates the mea- surement process, and subsequently generates samples from the posterior distribution that arises through the combination of the known prior and the known likelihood. Gold-standard samples from this distribution are obtained via efficient Gibbs methods. These samples are then compared to those obtained by the diffusion posterior sampling algorithms. In addition, the Gibbs methods can be Monte Carlo estimates of objects that needed in the reverse
used to obtain arbitrary-precisst stochastic differential equation, such as the minimum- mean-squared--rroo denoiser or its Jacobian. Consequently, the framework also enables the isolation and quantification of the error attributable to the likelihood approximations in the conditional reverse diffusion. We have provided numerical results for three common diffusion posterior sampling algorithms applied four common inverse problems. A consistent theme across all tested algorithms is that they are not calibrated, which demonstrates that research into algorithms that perform better in this respect remains crucial. We invite other researchers to benchmark their algorithms our open implementation.
'With some slight abuse of notation, log Py is the unnormalized ground-truth log-posterior (10). Since the additive constant is the same across all methods, this ranking is valid.
10

---

Published as a conference paper at ICLR 2026
Reproducibility Statement We release online repository with complete algorithm implemen- tations and step-by-step instructions to reproduce all results. A containerized runtime enables one- command setup and fully automated execution via the provided scripts. Each algorithm is specified at a level that supports independent re-implementation: The main text precisely details Gaussian latent- machine sampling: and the appendix presents the Bernoulli-Laplace Gibbs method in implementation- aligned notation, together with practical optimizations required for acceptable runtimes. The appendix also enumerates all experimental settings, including the numbers of training/valiaation/e signals, the samples-pe perplee--pe-daaee per-datum for each sampler, and the exact grid-search procedure used to select hyperpa- rameters.
Usage of Large Language Models We used large language models to adapt passages of already- written text for readability and conciseness.
The research leading to these results has received funding from the European Research Council under Grant ERC- 2020-AdG MunLeean-10110000 and by the Swiss National Science Foundation under Sinergia Grant CRSI15 198569, The authors thank Ludovic Reymond for his valuable contributions preliminary experiments.
Acknowledgments
REFERENCES
Amira Alakhdar, Barnabas Poczos, and Newell Washburn. Diffusion Models in De Novo Drug Design, Journal of Chemical Information and Modeling, 64(19):7238-7255 September 2024.
Brian D.O. Anderson. Reverse-time diffusion equation models. Stochastic Processes and Their Applications, 12(3)3313322,, 1982.
Robert Bassett and Julio Deride. Maximum a posteriori estimators limit of Bayes estimators. Mathematical Programming. 174(1-2):129-44, January 2018.
Robert C. Blattberg and Nicholas J. Gonedes. A Comparison of the Stable and Student Distributions as Statistical Models for Stock Prices. The Journal of Business, 47(2):244-280, 1974.
Lea Bogensperger, Dominik Narnhofer, Ahmed Allam, Konrad Schindler, and Michael Krauthammer. A variational perspective on generative protein fitness optimization. In Forty-Second International Conference on Machine Learning, 2025.
Pakshal Bohra, Pol del Aguila Pla, Jean-Françoii Giovannelli, and Michael Unser. A statistical framework to investigate the optimality of signal-reconstruction methods. IEEE Transactions on Signal Processing, 71:2043-2088 2055,, 2023.
Benjamin Boys, Mark Girolami, Jakiw Pidstrigach, Sebastian Reich, Alan Mosca, and Omer Deniz Akyildiz. Tweedie moment projected diffusions for inverse problems. Transactions Machine Learning Research, 2024.
Gabriel Cardoso, Yazid Janati el idrissi, Sylvain Le Corff, and Eric Moulines. Monte carlo guided denoising diffusion models for bayesian linear inverse problems. In The Twelfth International Conference Learning Representations, 2024.
George Casella and Edward I. George. Explaining the Gibbs sampler. The American Statistician, 46 (3):167-174, August 1992.
Hyungjin Chung and Jong Chul Ye. Score-based diffusion models for accelerated MRI. Medical Image Analysis, 80:1 00479, 2022.
Hyungjin Chung, Jeongsol Kim, Michael Thompson McCann, Mare Louis Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. In The Eleventh International Conference Learning Representations, 2023.
Christian Clason, Tapio Helin, Remo Kretschmann, and Petteri Piiroinen. Generalized modes in bayesian inverse problems. SIAM/ASA Journal on Uncertainty Quantification, (2):662-684, 2019.

---

Published as a conference paper at ICLR 2026
Rama Cont. Empirical properties of asset returns: stylized facts and statistical issues. Quantitative Finance, 1(2):223-2) 22032223,222 2001.
Evan Scope Crafts and Umberto Villa. Benchmarking diffusion annealing based Bayesian inverse problem solvers. IEEE Open Journal of Signal Processing, 6:975-991 2025.
Giannis Daras, Alexandros G. Dimakis, and Constantinos Daskalakis. Consistent diffusion meets Tweedie: Training exact ambient diffusion models with noisy data. In Proceedings of the Forty-First International Conference Machine Learning, ICML'22. JMI..orrg 2024.
Luc Devroye. Random variate generation for the generalized inverse Gaussian distribution. Statistics and Computing, 24(2):239-246 December 2012.
Julius Erbach, Dominik Narnhofer, Andreas Dombos, Bernt Schiele, Jan Eric Lenssen, and Konrad Schindler. Solving inverse problems with FLAIR, 2025. arXiv.
Xuefeng Gao, Hoang M. Nguyen, and Lingjiong Zhu. Wasserstein convergence guarantees for general class of score-based generative models. Journal of Machine Learning Research, 26(43): 1-54, 2025.
Crispin W. Gardiner. Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences. Springer, Berlin, Germany, 1990.
Andrew Gelman, John Carlin, Hal S. Stern, David Dunson, Aki Vehtari, and Donald B. Rubin. Bayesian Data Analysis. Chapman and Hall/CRC, November 2013. ISBN 97804291 13079.
Clément Gilavert, Saïd Moussaoui, and Jérôme Idier. Efficient gaussian sampling for solving large- scale inverse problems using meme. IEEE Transactions on Signal Processing, 63(1):70-80, 2015.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. GANs trained by two time-scale update rule converge to a local Nash equilibrium. In Proceedings of the Thirty-Fii First International Conference Neural Information Processing Systems, NIPS' 17, pp- 6629-66666 Red Hook, NY, USA, 2017. Curran Associates Inc. ISBN 5781510860964.
Desmond J. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review, 43(3):525-546, 2001.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusion Probabilistic Models. In Proceedings of the Thirty-Fourthy International Conference on Neural Information Processing Systems, NIPS '20, Red Hook, NY, USA, 2020. Curran Associates Inc. ISBN 9781713829546.
Hongbin Huang, Minghua Chen, and Xiao Qiao. Generative learning for financial time series with irregular and invariant patterns. In The Twelfth International Conference on Learning Representations, 2024.
Samuel Hurault, Matthieu Terris, Thomas Moreau, and Gabriel Peyré. From score matching to diffusion: A fine-grained error analysis in the Gaussian setting, 2025. arXiv.
Ajil Jalal, Marius Arvinte, Giannis Daras, Eric Price, Alexandros G Dimakis, and Jon Tamir. Robust compressed sensing MRI with deep generative priors. In Advances in Neural Information Processing Systems, volume 34, pp. 14938-14954, Curran Associates, Inc., 2021.
Achim Klenke. Probability Theory; A Comprehensive Course. Springer International Publishing, 2020. ISBN 9783030564025.
Muhamed Kuric, Martin Zach, Andreas Habring, Michael Unser, and Thomas Pock. The Gaussian latent machine: Efficient prior and posterior sampling for inverse problems, 2025. arXiv.
Hoang Trieu Vy Le, Marion Foare, Audrey Repetti, and Nelly Pustelnik. Embedding Blake Zisserman regularization in unfolded proximal neural networks for enhanced edge detection. IEEE Signal Processing Letters, 32:1271-1277 2025.
12

---

Published as a conference paper at ICLR 2026
Jean-Marie Lemercier, Julius Richter, Simon Welker, Eloi Moliner, Vesa Välimäki, and Timo Gerk- mann. Diffusion models for audio restoration: A review. IEEE Signal Processing Magazine, 41(6): 72 72-84, 2024.
Xiang Lisa Li, John Thickstun, Ishaan Gulrajani, Percy Liang, and Tatsunori Hashimoto. Diffusion- LM improves controllable text generation. Advances in Neural Information Processing Systems, 2022.
Jaming Liu, Rushil Anirudh, Jayaraman J. Thiagarajan, Stewart He, K. Aditya Mohan, Ulugbek S. Kamilov, and Hyojin Kim. DOLCE: A model-based probabilistic diffusion framework for limited- angle CT reconstruction. In 2023 IEEE/CVF International Conference on Computer Vision (ICCV), pp. 10464-10474, 2023.
Ozan Özdenizci and Robert Legenstein. Restoring vision in adverse weather conditions with patch- based denoising diffusion models. IEEE Transactions on Pattern Analysis and Machine Intelligence, pp. 1-12, 2023.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Vetito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic Differentiation in PyTorch, In NIPS-W, 2017.
Emile Pierret and Bruno Galerne. Exact evaluation of the accuracy of diffusion models for inverse problems with Gaussian data distributions, 2025a. arXiv.
Emile Pierret and Bruno Galerne. Diffusion models for Gaussian distributions: Exact solutions and Wasserstein errors. In Proceedings of the Forty-Second International Conference on Machine Learning, 2025b.
Rian Po, Wang Yifan, Vladislav Golyanik, Kfir Aberman, Jonathan T. Barron, Amit H. Bermano, Eric R. Chan, Tali Dekel, Aleksander Holynski. Angjoo Kanazawa, Karen Liu, Lingjie Liu, Ben Mildenhall, Matthias Nießner, Bjöm Ommer, Christian Theobalt, Peter Wonka, and Gordon Wet- zstein. State of the Art on Diffusion Models for Visual Computing. Computer Graphics Forum, 43 (2), April 2024.
Mehrsa Pourya, Bassam El Rawas, and Michael Unser. Flower: flow-matching solver for inverse problems, 2025.
Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf. Autoregressive denoising diffusion models for multivariate probabilistie time series forecasting. In International conference on Machine Learning, pp. 8857-8868. PMLR, 2021.
Mengwei Ren, Mauricio Delbracio, Hossein Talebi, Guido Gerig, and Peyman Milanfar. Multiscale structure guided diffusion for image deblurring. In 2023 IEEE/CVF International Conference Computer Vision (ICCV), pp. 10687-10699, Los Alamitos, CA, USA, October 2023. IEBE Computer Society.
Severi Rissanen, Markus Heinonen, and Arno Solin. Free hunch: Denoiser covariance estimation for diffuuson models without extra costs. In The Thirteenth International Conference Learning Representations, 2025.
Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Bjorn Ommer. High- resolution image synthesis with latent diffusion models. In 2022 IEEE/CVF Conference Computer Vision and Pattern Recognition (CVPR), pp. 10674-10685. IEEE, June 2022.
Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J. Flect, and Mohammad Norouzi. Image super-resolution via iterative refinement. IEEE Transactions Pattern Analysis and Machine Intelligence, 45(4):713 2023.
Ken-Iti Sato. Lévy Processes and Infinitely Divisible Distributions. Cambridge studies in advanced mathematies. Cambridge University Press, Cambridge, England, 1999.
Otmar Scherzer, Markus Grasmair, Harald Grossauer, Markus Haltmeier. and Frank Lenzen. Varia- tional Methods in Imaging. Applied mathematical sciences. Springer, New York, NY, 2009 edition, October 2008
13

---

Published as a conference paper at ICLR 2026
Wim Schoutens. Lévy Processes in Finance. Wiley Series in Probability and Statistics. John Wiley & Sons, Chichester, England, 2003.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In, Proceedings of the Thirty-Second International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp.
2256-2265, Lille, France, 7 2015. PMLR.
Jiaming Song, Chenlin Meng, and Stefano Ermon. arXiv:2010.02500 October 2020.
Denoising diffusion implicit models.
Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz. Pseudoinverse-guided diffusion models for inverse problems. International Conference on Learning Representations, 2023.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Advances in Neural Information Processing Systems, pp. 1 11895-11907, 2019.
Yang Song, Jasha
Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021.
Stanislas Strasman, Antonio Ocello, Claire Boyer, Sylvain Le Corff, and Vincent Lemaire. An analysis of the noise schedule for score-based generative models. Transactions on Machine Learning Research, 2025.
Andrew M. Stuart. Inverse Problems: A Bayesian Perspective. Acta Numerica, 19:451-559, 2010.
Yu Sun. Zihui Wu, Yifan Chen, Berthy T. Feng, and Katherine L. Bouman. Provable probabilistic imaging using score-based generative priors. IEEE Transactions on Computational Imaging, 10: 1305, 2024.
David Y. W. Thong, Charlesquin Kemajou Mbakam, and Marcelo Pereyra. Do Bayesian imaging methods report trustworthy probabilities?, 2024. arXiv.
Michael Unser and Pouya D. Tafti. An Introduction to Sparse Stochastic Processes. Cambridge University Press, 2014.
David van Dyk and Taeyoung Park. Partially Collapsed Gibbs Samplers. Journal of the American Statistical Association, 103(482):790-796 2008
Roman Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.
Martin J. Wainwright and Eero P. Simoncelli. Scale Mixtures of Gaussians and the Statistics of Natural Images. In Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999.
Zheng Wang, Johnathan M. Bardsley, Antti Solonen, Tiangang Cui, and Youssef M. Marzouk. Bayesian inverse problems with & priors: A randomize-thennopiimizz approach. SIAM Jour- nal on Scientific Computing, 39(5):S140SS166 140-6166, 2017.
Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli. Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 13(4): 600-612, 2004.
Zhuohan Wang and Carmine Ventre. A financial time series denoiser based on diffusion models. In Proceedings of the 5th ACM International Conference on Al in Finance, ICAIF '24, pp. 72-80, New York, NY, USA, 2024. Association for Computing Machinery. ISBN 9798400710810.
Joseph L.. Watson, David Juergens, Nathaniel R. Bennett, Brian L. Trippe, Jason Yim, Helen E. Eise- nach, Woody Ahern, Andrew Borst, Robert J. Ragotte, Lukas F, Milles, Basile I. M. Wicky, Nikita Hanikel, Samuel J. Pellock, Alexis Courbet, William Sheffler, Jue Wang, Preetham Venkatesh, Isaac Sappington, Susana Vázquez Torres, Anna Lauko, Valentin De Bortoli, Emile Mathieu, Sergey Ovchinnikov, Regina Barzilay, Tommi S. Jaakkola, Frank DiMaio, Minkyung Back, and David Baker. De novo design of protein structure and function with RFdiffusion. Nature, 620(7976): 1089-1100, July 2023.
14

---

Published as a conference paper at ICLR 2026
Xingyu Xu and Yuejie Chi. Provably robust score- score-based diffusion posterior sampling for plug- and-play image reconstruction. In The Thirty-Eigh Eighth Anmual Conference on Neural Information Processing Systems, 2024.
Duoduo Xue, Wenrui Dai, Ziyang Zheng, Xinyu Peng, Junni Zou, and Hongkai Xiong. Fourier phase retrieval with diffusion priors guided by failure detection. In 2025 IEEE International Symposium Circuits and Systems (ISCAS), pp. 1-5, 2025.
Nebiyou Yismaw, Ulugbek S. Kamilov, and M. Salman Asif. Covariance-correcte diffusion models for solving inverse problems. In 2025 IEEE Statistical Signal Processing Workshop (SSP), pp- 26-30, 2025.
Richard Zhang, Phillip Isola, Alexei A. Efros, Eli Shechtman, and Oliver Wang. The unreasonable effectiveness of deep features as a perceptual metric. 2018 IEEF/CVE Conference on Computer Vision and Pattern Recognition, pp. 586-595, 2018.
Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool. Denoising diffusion models for plug-and-play image restoration. In IEEE Conference on Computer Vision and Puttern Recognition Workshops (NTTRE),
A BAYES ESTIMATORS
A benefit of the Bayesian approach over classical variational methods (see, e.g., (Scherzer et al., 2008)) is that different point estimates arise from a fixed prior. For a given measurement y, these point estimates summarize the posterior distribution Px|Y=y with respect to given loss f: Rd x R - R via the optimization problem of finding the point (y(y) that minimizes the posterior risk:
(/ ((xxx) (XIxxxxxx))
< = arg min *ER
(18)
In this paper, the Bayes estimator with respect to the mean-squared error (MSE) Ł= -11 - 1
plays key role due to its close relation to the prior score in the reverse diffusion (see Section 2) and because we quantify the performance DPS algorithms via the MMSE optimality gap in Section 4. With this choice of (18) can be written as
xpx Y=y (x) dx E|X Y y], (f{ 1 -||x-x|1 xx| dxx ffe
SMMSE (y)
arg min &ERd
3d
(19)
which is the expectation of the posterior Yx
Another widely-used estimator arises through the choice
E(x, x) = - (x)
(20)
where
Xa(x)= 11 if
(21)
else,
which leads to the MAP estimator that seeks the mode of the posterior:
(L
-X(x)(x) (X(xxs)(xx) PXIY=y dx
8MAP(Y) = arg min
=
arg max PX Y Y-y (8).
(22)
SER
RERd
Rewriting (22) as
XMAP(Y) = arg min
A y|12
- -lggpx(())
(23)
20
R
reveals close relation to classical variational approaches after identifying the regularizer with -log
"This definition is informal but sufficient for the purposes of this paper. For continuous posteriors, the strict 0-1 loss Bayes' rule is ill- posed. A common formalization defines MAP as the limit of Bayes estimators under shrinking small-ball 0-1 losses; under additional regularity, this limit agrees with the posterior mode (Bassett & Deride, 2018; Clason et al., 2019). The MAP estimator may also not be unique.
15

---

Published as a conference paper at ICLR 2026
B TWEEDIE'S FORMULA
In the setting of Section 2, we now derive an equality that relates V Plogxx, to E[x, X, = I.i.e
the MMSE estimate of X. given that X, takes on certain value. Similar derivations can be found in, e.g., (Song et al., 2021; Chung et al., 2023; Daras et al., 2024), but we include it to underscore the relevance of the MMSE estimate in this paper and facilitate the understanding of its relation to various objects. Under the suriance-preseree
201 x and g(t) = v ((t) choice for f(x, t) = 2
of the drift and diffusion coefficient, the diffusion SDE (4) simplifies time-inhomogeneous Orastein-Uhlenbeed SDE (see Klenke (2020, Example 26.5))
dx;
8(0) X, dt v3(t)dW..
(24)
whose pathwise solution
X, = xt)Xx + L.
((s)) /B(t)d...
(())
(25)
a(s)
where X, is an appropriate initial condition and a(1) = exp(- - So ((s) ds). can be computed with standard techniques, see, e-g., (Gardiner, 1990, Section 4.4.4). In addition, since
(( (a(t) (8)
B(t) ds = l B(8) ex)) (-/
B(u) ds = 1 a2),,, 1-)2((t,
(26)
we can write that
X, a
ot)N
(27)
in distribution, where (t) = (1-a (1-a((tt) Consequently, the density of X, is given by the convolution of px. with a Gaussian with variance (1) and appropriate scaling by at),, which we write
/ Rd
px. (x)
( (())))pxx, (x) (2) di.
(28)
where xxx) = (2T) that
Σ exp - f| - 1el 11) Finally, after taking the gradient, we see
(x) =
(0, (())))xxxx(()) (8) (x) d&
Rd
o2(1) 1 (x (xx((())))00001(111))))) 90, ((t))))xxxx)))xxxx166 dx
(29)
o²()) xpx (x) aa1)) / 890,02 a(l)x)xx (x) (x)d
Rd
0(0)2 (xpx, (x) (x) X. xx].
Finally, after dividing by Px, (x) and since Vpx, (x)
= V log Px, (x), we find the celebrated Tweedie
Px, [x]
identity
v loggx, (x) = -o(t) 2 a(t)E[Xo X, =x).
(30)
B. A CONNECTION BETWEEN THE DISCRETIZED REVERSE SDE AND DDPM
To show the connection between the Muler-Maruyama discretization of the reverse-diffusion SDE and the DDPM backward process, start by deriving the latter from the respective forward process. DDPM has been introduced by Sohl-Dickstein et al. (2015) as a discrete-time Markov chain of length T with Gaussian transitions
Px, Xx_ = Gauss /1 - 3,x, B.T)
(31)
such that the transitions from X, to X, are also tractable as
X, Va,Xa + VI-ā,
(32)
where ai = =(1 (1- (1--3i) -.)) āg = -I'' a2 and Ze ~ Gauss (0, I). By definition,
X, v1 B,X, v B,2.1
(33)
16

---

Published as a conference paper at ICLR 2026
and straightforward application of Tweedie's ' formula (6) gives that - X1]
Voa X + (1 a,) lx (X:)),
(34)
which leads to the DDPM backward transitions
X -í-Ø (X, + 8 logpx, (Xi)) v B.Z.
(35)
like they appear in (7).
Now, we recall the reverse-diifusion SDE which, under our choice of the drift and diffusion coefficient, given by
dx; (- [(t) X, 3(1)V logpx. (X:)) di + VB(() dw, A first-order step from 1 to (1 - 1) (dl = 1) of gives the Euler-Maruyama update
(36)
X (1 + ) 3. Px (X,) + VB,Z,.
(37)
where B. = (()) and Z, ~ Gauss (0, I).
The DDPM reverse process (35) can be related the the Fuler-Maruyama discretization of the reverse SDE (37) via Taylor expansions, since
B
1+
0(3))
(38)
B
2
and
B1
Br 0(3))
(39)
V1
3e
as Be + 0.
C LÉVY PROCESSES AND INCREMENT DISTRIBUTIONS
The prior distributions in our framework are those of signals obtained by regularly spaced samples of processes with independent, stationary increments (Lévy processes and their discrete-ting time counter- parts). We briefly recall the definition; see Unser & Tafti (2014); Sato (1999) for background and the link infinitely divisible laws.
Definition C.1 (Lévy process). A stochastic process (s(t) 2 0} is Lévy process if
1. (anchor at the origin) It holds that x0) = 0 almost surely:
2. (independent increments) for any N E N {0, 1 and 0 < < t2 ... < ty < xc, the increments (s(tt)=((t)))) s(tz)), (s(ta)=s(t2) ... (s(tx)- -s(x))1)) are mutually independent;
3. (stationary increments) for any given step h, the increment process un = {s(() - s(t - h)
> h) is stationary;
4. (stochastic continuity) for any
0 and 2 0,
lim PrpP (Is(t + h) s(1) > 5) = k=0
We form discrete and finite-length signals by sampling integer times and stacking the values
into x s(1),,((),....... ((2),....,(())) s(d)). Let the unit-step increments be [u] = (s(k) - s(k 1)) for k 1, 2, ... d. By independence and stationarity, the law of [u]x does not depend on & and we denote it pu. We define the fiiieeddifferenc matrix
0
-
1
1
0
D
1
(40)
-1
1
Forr choices, it always has density suitable reference measure.
17

---

Published as a conference paper at ICLR 2026
Table Univariate distributions used throughout this work. Parameters appear in the order they are specified in this table, e.g. Gauss(p o²).
Name
Distribution
Parameter(s)
Supp.
Notation
Gaussian
exp
/ R, o²
Rg0
R
Gauss
22ma²
o
Exponential
A Texp(
Ax)
À E Rs0
R20
Exp
Laplace
25 exp
[x|
b E Ry0
R
Laplace
u+4
Student-t
+
v E Ry0
R
St
V
Gamma
8 Mla ta) exp( -3x) (2)
a, B E R>a
R>0
Gamma
Gen. inv. Gaussian
p- exp
a2+bb/2
a, E E R
h>0
GTG
(Vp(vab) 18(x) - X) 26
Bernoulli-Lppace
exp
X [0, 1|. R
R
BL
fo .ー exp (-() dt for any x E IR The modified Bessel function of the second kind with parameter denoted by Ky.
Moreover, the gamma function is defined as ((r)
such that the increment vector satisfies
Dx.
(41)
Because s(0) = 0, the finiteddifference matrix D has an initial condition that makes it invertible and D-1 lower-triangular matrix of ones. This also implies that for all 1,2,...,d.
k
[x]
(u)an
(42)
which is a convenient way to synthesize signals once u is drawn. The combination of (41) with the independence of the increments implies that the density of the discrete signal
d
(x = II pu (Dx]x)-
(43)
k=1
C.I EXTENSIONS
The approach in this paper can extended two- or nigher-dimensionn signals on grids, such as images or videos, and even to more specialized structures like signals defined over trees or graphs. The structure the signal is effectively encoded through the choice of the matrix D. For instance, two-dimensional finite-difference matrix would result in signal vector that can be interpreted as woo-dimensional image. The main additional (computational) challenge is sampling during signal generation: Whenever D is not trivially reducible to dieedmensiooall operator, the model (43) will be overcomplete and, in general, no whitening transformation exists to decouple increments for independent sampling. The extension higher-diment digees-ainensional signals and the complications that arise in that context are rigorously treated in Kuric et al. (2025).
C.2 LATENT DISTRIBUTIONS AND NOTATION
Some of the distributions that we rely on in this work have multiple competing parametrizations. To avoid ambiguities, we provide precise definitions of the four increment distributions that we consider in this work: Gaussian; Laplace: Student-t; and Camoulli-Laplacee
spike--andsslaby We give in Table our notations of these and other distributions that we use in this work. We list in Table latent maps and conditional latent distributions that are needed for the GLM for the distributions in this work.
D GIBBS METHODS AND SAMPLING EFFICIENCY
Gibbs methods are Markov chain Monte Carlo (MCMC) methods to sample from joint distribution Px. Zn Zz....... Z. of (n+D blocks of variables that are advantageous when the direct sampling compu-
18

---

Published as a conference paper at ICLR 2026
Table 3: Latent variable representations and conditional distributions for common distributions.
Dist. Φx
Latent dist. fk
Latent maps
Cond. latent dist. P[z]a| X [Kx]
Gauss(ju, o²) Laplace(b)
8(0)
[=<(2) = (2) = x²
8(0)
Exp 26
(z( = 0, (z) = 2
GIG (+r Kx
k 2
St(v)
Gamma (2.
pu(2)) = O, of(2)
Gamma
p+1 2
y
Kx
Algorithm 3 Latent-variable Gibbs sampling of px ".....
Require: Burn-in period B e N, number of samples S E N, initial point (Xo, Zi.... Zn).
I: for k = 1....... B + do
2: X5 Px Z. 22..........
3: 21 X=x|xxxx
Latent blocks do not need to be stored
4:
5: return {xB++]Ā=1
tationally difficult but sampling from the conditional distributions PX|z Z, 22......... ZZ. X, Z2........ is easy. Gibbs methods cycle through the conditional distributions with repeated draws, which maintains the joint distribution invariant (Casella & George, 1992). The naming of the variables X, 21, Ze.... Zn is deliberately chosen to emphasize that we use latent-variable Gibbs methods that rely on auxiliary variables that are introduced solely to make the conditionals simple. The steps of general latent variable Gibbs sampler are shown in Algorithm 3, where the iteration counter in the sampling the latent variables is omitted since they need not be stored and previous iterations can immediately be overwritten.
Z..
Kuric et al. (2025) recently showed that such methods are significantly faster than other standard sampling routines that are commonly used in settings similar to the one in this paper. They report sampling efficiencies of close to 1, while alternatives, such as the Metropolis-adjusted Langevin algorithm, achieve sampling efficiencies? of around 1 x 10-3 3_ In addition, Gibbs methods require no stepsize or acceptance-rar tuning and introduce no discretization bias. These properties motivate our of Gibbs methods for the fast and robust posterior sampling throughout this work.
Like all MCMC methods, in practice Gibbs methods benefit from the discarding of some number of initial samples (the burn-in period) when the initial point is located in low-density regions. After the burn-in period, the quality of the Monte Carlo estimate of any object depends on the number of samples one uses in their estimation. We discuss our choice of the burn-in period and the number of samples for the various problems in Appendix F.2.
D.1 A GIBBS METHOO FOR BERNOULLI-LAPLACE INCREMENTS
Let be the Dirac distribution. Then, letting À be the Bernoulli parameter and b the scale parameter, we note that the Bernouli--applcee density
Po (u) Aō(u) (1 2) exp(-bbl
(44)
admits the representation
(u) = 11Σ PUiV= D. (ww) dw.
(45)
1=0
where
pv(u) = x2- 2)
(46)
for E {0. is a Bernoulli distribution.
5²
b²w
Pw (w(w)= =
2 exp
XRzo (w)
(47)
2
SSmpling efficiency refers to effective samples per iteration; an efficiency of means roughly 1/p iterations per "effective sample" (Gelman et al., 2013, Section 11.5).
19

---

Published as a conference paper at ICLR 2026
Algorithm 4 Bemoollii-Lplace Gibbs sampler.
Require: Initial increments Ug E IRd, initial support vector v E R 1: for = 1, 1,2,.... ... B + S do
2: Draw [w] Pw] wwwj-ju.-lll.V-W
parallel over
3: for k = ... do
4:
Draw [v]x ~ Bernoull((px(v))
5: Draw u,
6: return (D uB=
an exponential distribution, and
[ō(
if o 0,
PUlV=r,W=w(C0 Gauss(O, w)
(48)
if
The algorithm relies on the introduction of two latent vectors v, E Idd that satisfy
PUIV=y (u) II
Wo ([ul))
(49)
such that, as a result, the distribution conditioned on the measurements can be written as
d
PU.V PU. w) x exp(
1 Hu 20
y11
II
PulV= ([u]a)
1=1
(50)
6²
ba
w/
x II
1- M
2)(v
exp
2
2
1=1
1=1
where AD-1 Equations (48) and (50) imply that any sample from PUIV w. w.W.wwWW-.......... Y- takes the value zero at those indices where is zero, and values from multivariate Gaussian distribution with covariance C (oHHH* diag(w))
and mean o, "CH'', otherwise. Sampling W = u, V v, = amounts to the independent sampling of d one-dimensional distributions, which are Exp(2/b2) indices where is zero and GIG [ulß. 0.5) those indices k where is one. The conditional distribution of the binary support vector is
TV/ (v) & B(v, w) + expp(-y B(v, ' yy II (1 k=1
(51)
where B(v, = Hdiag(v ④ H The standard way to sample from this distribution is coordinate-wise Gibbs sampler that updates [v]s Bernoulli(ps with
pv(v) = (1 + exp( exp(-As(v)))-
(52)
where the log-odds increment
A&(v) = log 1-1 x -2 (log BB(v==1, w)l log B(( w)l)
(53)
B( - B(V w) -))
where Vk= (V1, (V1,..., Vk Vk-1, Vk+ ...,, va) is the difference between the log-posterior when the bit is on and when it is off. The resulting algorithm is summarized in Algorithm 4 and can be interpreted as (d + 2)-block (i.e,, dimension-dependdent) Gibbs method.
'This formulation is equivalent to the one presented by Bohra et al. (2023), who explicitly "slice" the matrices H and diag(w) with the indices where is We stick to this formulation since requires less notation and emphasizes that implementations need build variable-sized matrices, which crucial for efficient implementation on modem compute units that utilize highly parallelized computations.
9This is only interpretation hecause the density violates the classical positivity conditions that needed for Gibbs methods. is partiolly collapsed Gibbs method, see (Bohra al., 2023; Dyk & Park, 2008).
20

---

Published as a conference paper at ICLR 2026
Chalesky -C- Perturb and MAP
0.025
10a
Time seconds
cOOOE
10"
10-2
4.000
8000
4000
4000
8000
4000
8000
d
d
Figure 6: Runtimes needed to perform 20 Gibbs iterations on denoising posterior (Laplace 11) increment distribution, 10 parallel chains) depending on the dimensionality of the signal. Missing entries are due to excessive memory requirements.
D.2 PRACTICAL GIBBS IMPLEMENTATIONS
Sampling the Gaussians The sampling of X Z in the GLM and of U V, W, Y for the Bernoulii-Laplace case reduces to drawing from high-dimensional Gaussian, which a well- studied problem. For settings that necessitate matrix-fiee -fee implementation such those that commonly encountered imaging applications, Kuric et al. (2025) advocate a Perturb-and-Man sampler with preconditioned conjugate- -anjugate-gradient solvers. We report the runtime of the Gibbs method as a function of signal dimension for Taplace(1) increment distribution in Figure 6. A standard implementation based Cholesky factorization of the covariance matrix matrix-which requires explicitly instantiating the matrices A and D in memory-i is faster than the Perturb-and-MAN sampler with geguuatte-gradeen solver across a broad range of noise variances and dimensions. For our moderate- dimensional setting with d = 64, the Cholesky- -ased implementation is approximately an order of magnitude faster.
However, explicitly storing these matrices becomes infeasible at larger dimensions (in our setup, we ran out of memory at d = 8096), and the expected cubic scaling is apparent in the figure. In contrast, the Perturb-and-MA.
sampler (convergence criterion: squared residual norm below x 10-9), while slower than Cholesky small dimensions, exhibits substantially better scaling with signal dimension. In particular, it does not require materializing the operators: both the measurement operator A and the finiteddifference operator D can be implemented efficiently in a matrix-free manner. Moreover, the sublinear runtime observed in this experiment suggests that the corresponding linear systems well conditioned.
The sampling accuracy of Perturb-and-Man depends the termination criterion used by the optimization solver, and any finite stopping rule yields approximate samples. A principled refinement correction step remove bias, and tune the solver accuracy
is incorporate Seetropoiis-Hasting to optimize overall runtime; this strategy was proposed by Gilavert et al. (2015), to which we refer for details. Overall, these results indicate that the Gibbs method scales favorably to higher dimensions. Combined with the fact that the denoising posteriors are sub-Gaussian and with the relatively mild requirements for estimator accuracy this setting, this suggests that the proposed framework scales well as the dimension increases.
Sampling the Latent Variables The sampling of the latent variables necessitates the sampling of the -iddennnsional conditional latent distributions. All the conditional latent distributions that are relevant in this paper admit efficient samplers that are readily available in standard scientific computing packages or can be implemented with little effort. We reuse the CUDA implementation of the generalized inverse Gaussian sampler from Kuric et al. (2025) that implements the method proposed by Devroye (2012) and rely on PyTorch (Paszke et al., 2017) for all others. Wherever possible, latent updates are parallelized.
In the Gibbs methods for the CemoolliiLaplaan increments, the sequential drawing of the binary support vector V is embedded in the outer Gibbs loop which, in turn, may be embedded in the reverse- diferss-iffu loop. This makes it crucial to minimize the use of heavy algebra operations to achieve acceptable runtimes. Writing B(v, w) of + Hdiag(v H) we recognize that
21

---

Published as a conference paper at ICLR 2026
the flipping of the kth bit of v adds or removes rank-one term [w].H,,H.. where Hş is the kth column of H. Using the matrix matrix-determinan lemma and Woodbury-Sherman Sherman-Morrison, we update
log B( w) = log (V w)l + log(1 + [w]eTk)
(54)
and
y B( w) = o( w) Wwa (H, w)
(55)
where ты H B( w) HHg. Thus, efficient implementation factors B(v, w) once per latent state, obtains the needed scalars via triangular solves, and performs rank-one updates as bits
flip. We report our cumulative runtime improvement over naive implementation in Figure
E
DPS UPDATE STEPS
E.1 COVARIANCE IN C-DPS
C-DPS (Chung et al., 2023) uses the approximation of the likelihood
py(x,=x(y) PY|X, =Xx
(56)
When the noise in the inverse problem is Gaussian, the likelihood score
(xx log PY EX/ (y)) necessitates the computation of
V (x 2 AEEX X, x] - -
(57)
which is
J (x E X, ()) (AE
X, = y)
(58)
after application of the chain rule. The Jacobian J (x - E X, x)) is typically computed with automatic differentiation when (x, t) + E,Xo X, x] is approximated with a neural network. In our framework, use the connection with the covariance matrix Cov [Xo X, -]. Indeed, as also shown by Rissanen et al. (2025), if Xo and X, verify (32), then
T- Cov(Xo X, x] a. (I+
2 lxg (x))
(59)
This identity, combined with the derivative of (6), yields
V Or
J (x+ X, xl)
Cov|Xo X; =x]]
(60)
E.2 EXPLICIT UPDATE STEPS
We give the instantiations of the update step S a variety of DPS algorithms below. Each z, is d-dimensional random vector with i.i.d. standard Gaussian entries.
Score-ALD (Jalal et al., 2021) The input parameters of this algorithm composed of the follow-
ing: A noise schedule ( the noise level of the inverse problem on, and annealing parameters (n, 17-1 and (% 1T-1 The update step goes
1
X
5 a=1
(61)
8, ( x
1. A A Axi)) + 22 X X₁ + Me (s: + +ot
C-DPS (Chung et al., 2023) The input parameters are the the variance variance-preservee preserving scaling weight Gr in (32), the variance of the diffusion transitions 8, as in (35), and a scalar Ç that governs the
22

---

Published as a conference paper at ICLR 2026
likelihood-gee -giidance strength. The diffusion noise level that corresponds to the denoising posterior is denoted o, = (1 at)/ vā,, which is used to compute the samples {xx 18 The update step goes
11
Xxy
k=1
S
1
C
( (8x - p
8
(62)
1=
x/-1
Var(1-ar-1)
Xt + -Be-118 18 fa + 01hi,
1-6r
1-
x1-1
X(-1
TAp
1-o,
C A (Ap--pp)) yy.
Xt-1 8t-1/v 01-1-
DifIPIR (Zhu et al., 2023) The input parameters are similar to those of C-DPS. It also uses the noise level of the inverse problem on and additional balance hyperparameter 7, The update step goes
1 11
X/.
S
Pt
of
Xn arg min [ Ax - y112 + 2 1|x -m1l?,,
(63)
xE2#
== vT-ai 1 (x: vārXo);
8, V0 + v ar-ily zē vize),
X,-1 x Vo--i
HIGDM (Song et al., 2023) The input parameters are the noise schedule (o.) the data- dependent noise schedule (r:)r-> and the DDIM (Song et al., 2020) time-dependent coefficients (ct (1) T-11 and (g 1-= The update step goes
a
k=1
1
C =
-x) - -R) T.).
5 $=1
V C,
(64)
1 a².
= (y AW) T AAT + = AJ) AJ)",
ai (x: - v Va
_ Vo. C1 + ozê Varg
X--1 //-1/ vār-i.
DPnP (Xu & Chi, 2024) The diffusion noise level that corresponds the denoising posterior is denoted 1), which is used to compute the sample X1. This same t)2 defines the dielihood-guidanee strength. The update step goes
Xo = X1.
(65)
Xt-1 exp -
-yl
I
-xoll2)
247
23

---

Published as a conference paper at ICLR 2026
Diffision P...or
Sumjllnn
Denoising
Likelhood
PAglN Ny PXg Yay
Pkgy Aymay Ganseis,
Pag
Viis Pyrix.. (X)
[Ru]
Xo-1/2
Xe-1" Re-173
Siep 5 of DPS aguerilm
Figure 7: Illustration of the proposed template for DPS algorithms.
Annealed Plug-and-Paa Monte Carlo (PnP and RED variants) (Sun et al., 2024) The diffusion noise level that corresponds to the denoising posterior denoted Øs which are used compute the samples {x ) ==1 The parameter 1 denotes the likelihood guidance strength, and Ky is an annealing parameter.
The update step for the PnP variant goes
1
Xk:
$ (n xx) x/)/oe
(66)
X4-1 = - 7848 + 2
X4 X 7A (Ax y).
The update step for the RED variant goes
p 5 Σ 8k:
(67)
(p - ///
Xx - (A (AX -KKss + /2722.
The DPS template that is summarized in Algorithm is illustrated with one-dimensional exampee Figure 7.
F NUMERICAL EXPERIMENTS
F.1 FORWARD OPERATORS
We consider four forward operators A our experiments. The first operator the identity A I Rdxd This choice is motivated by the fundamental role that denoising algorithms currently play in many restoration algorithms and even in labeling problems such edge detection (Le et al.. 2025). The second operator A E ddd implements the convolution with kernel that consists of the 13 central samples of truncated Gaussian with variance that are normalized to unit sum. We adopt circular boundary conditions to enable fast computation of the proximal map that arises in the update step DiffPIR (see Appendix E.2) via the fast Fourier transform. Deconvolution a relevant problem with applications like microscopy astronomy. The third operator is sampling operator A mxxd that returns m < d entries of its argument unchanged. This operator is relevant many fields such image reconstruction and time-series forecasting. In particular, in forecasting prediction problem the operator would return the first m known entries, and the resolution of the inverse problem estimates the remaining (d m) entries. In experiments, each entry has independent chance of 40% % being kept. The fourth and last operator A MF Immxd where 2 (22(((22)+2))))) is the matrix representation of the "real" one-dimensional discrete Fourier transform with separated real and imaginary components, and M E x -1) is a sampling operator. Such operators are relevant in medical imaging astronomy. The sampling operator is constructed such that the lowest
24

---

Published as a conference paper at ICLR 2026
frequencies (the DC term included) are acquired, while the remaining frequencies independently have a 40 % chance of being kept.
For all operators, the noise variance of is chosen such that the median measurement signal-to-noise ratio (SNR) is around 25 dB. We set Ntrain = 1 x 10,, Nval = 1 x 10*, and Nuest = 1 10".
F.2 BENCHMARK IMPLEMENTATION
The benchmarking pipeline starts with the generation of Neest test signals denoted (xx==) Niett
I /== per increment distribution, each of which is independently synthesized by first drawing i.idd increments from the respective increment distribution and forming the signals via (42). It then proceeds to Muas synthesize the Ntest measurements (i.e. we use one noise instance per signal) denoted {yk=*t ]==l according to (1) and, for each of the measurements, computes the gold-standane posterior samples of the various inverse problems via the Gibbs methods described in Section 3. This stage is off-line (no reverse-diffusion loop) and trivially parallel across the measurements, which allows us to run long chains with burn-in periods of 1 x 10s iterations and obtain 2 x 10y draws from the posterior distribution. This far exceeds any values reported by Kuric et al. (2025) or Bohra et al. (2023) and results in precise MMSE estimates.
The datasct- stage also involves the generation of Ntrain training signals {x Extaa Nusin and Nval validation signals (mutually disjoint from the test signals) {xx M Nou. along with the corresponding validation measurements {yr ===1- The training signals are used for the learning of neural score function like those that used for the resolution of inverse problems when the prior is unknown or too expensive to evaluate. Training details are provided in Appendix F.3 The validation signals are used to monitor the performance of the neural score function on unseen signals during the training stage and to tune the regularization parameters for the model-based approaches as well as the parameters of the DPS algorithms, see Section 4.1
Jk=1
Noat
Unlike for the computation of the gold-standard MMSE estimate of the initial inverse problem, the denoising posteriors are sampled times per trajectory (we use T 1000). To ensure acceptable runtimes in this setting, we therefore pick the smallest burn-in period and sample count that still yield accurate estimates of the required statistics. We determine these settings with rigorous protocol that is detailed in Appendix F.4. Ultimately, this protocol resulted in the choice of a burn-in period of 100 iterations and sample count of 300.
F.3 LEARNING DETAILS
For learned-based denoisers, a noise-conditional neural network with UNet architecture (305 761 learnable parameters) is trained in an off- step on the Nursin training signals standard setup (Adam optimizer with learning rate 1 x 10 with exponential decay with factor 0.9999, 100 000 parameter updates, batch size 10 000). The noise schedule in C-DPS and DifPPIR is defined by the two endpoints So 1 10-4 and 3T 2 x 10-2 with linear equidistant samples in-between. The learned variant of DPnP is the "DDS-DDPM"' variant (Xu & Chi, 2024, Algorithms and 3) that contains inner denoising-samplinn loop. The arbitrary-preciste variant does not require inner loop at all (except for the burn-in period), which makes the arbitrary-precision variant the faster one for this case.
F.4 BURN-IN PERIOD AND THE NUMBER OF SAMPLES
As discussed in Appendix F.2, the burn-in period and the number of samples of the Gibbs samplers needs to be chosen appropriately to ensure an acceptable runtimes and a sufficiently small Monte Carlo error. We determine the burn-in period and the number of samples through the following protocol that is run in an off- stage prior to running the benchmark. We synthesize = Xo o,n where or is in the range defined by the noise schedule B. Xo is constructed via (42) for all four considered increment distributions, and n is some unknown but fixed vector of standard Gaussian noise. For each of the synthesized signals, we then launch C = 1000 parallel Gibbs chains on the corresponding denoising posterior and run those chains for Neufficient iterations, where Neufiiient is a sufficiently large natural number that guarantees that the chains are stationary for at least Navg (which is also relatively large) iterations and that, consequently, we compute precise estimates of
25

---

Published as a conference paper at ICLR 2026
d= 1.99
20.27
0
80.29
151.18
Gauss 0.25)
Laplare(1)
B1.00.11)
10a
Mean conn--uurreed
109
10
St(1)
St(2)
St)
Mean- squarsuurrrel
10a
10
1,000
500
1,000
500
1,000
Window longth
Window longth
Window length
Figure 9: Mean squared error between MMSE estimates and reference MMSE. Dashed lines: Learned neural MMSE estimate. Solid lines: Monte Carlo MMSE estimate in terms of the window length.
various statistics the posterior distribution from the iterates from the last Navg iterations across all C' chains.
To determine the burn-in period, we then pro- ceed to calculate statistic that we mon- itor throughout the iterations and that we can compare against the reference statistic. Specif- ically, denoting with X the random variable of the Gibbs sampler, we compute the empirical distribution of the increments at index 32 like (Xa3 X.2). The distribution of differences that is obtained by taking the last Navg iterations across all C chains is considered the reference distribution. Then, compute the Wasserstein- 1 distance of that distribution to the one obtained by taking the average across Navg iterations and all C' in a sliding-window starting from the first Gibbs iterations. This allows us gauge the burn-in period through a visual inspection of the Wasserstein- distance through the Gibbs iterations. In particular, we expect the Wasserstein-1 dis- tance to be large for a number of initial samples where the Gibbs sampler is not stationary and then to oscillate around small but nonzero value. The value will be nonzero due the finite sample size. The Wasserstein-1 distance between the reference distribution and the one obtained through the Gibbs iterations is shown in Figure 8 (for the exemplary case of a St(1) distribution and a selection noise variances). We observe that the empirical distribution of increments converges rapidly to the reference one. The Wasserstein-1 distance reaches the noise level after a single-digit number iterations, which is in line with the analysis provided by Kuric al. (2025). Based on these findings, chose the burn-in period as B 100 iterations for all our experiments, which is more than sufficient to reach stationarity and has acceptable runtime.
6.17
10a
Wassestein- distance
10a
100
Now wes
10
100
200
300
Figure 8: Wasserstein-1 1 distance of intermediate marginal distributions to that of the final sample.
26

---

Published as a conference paper at ICLR 2026
To determine the number of samples that needed for sufficiently accurate computation of various statistics that any DPS algorithm may utilize in their update steps, we compute precise estimation of the MMSE estimate of the denoising posterior by averaging the last Navg iterations across all C chains. Then, we pick one arbitrary chain and grow a window from iteration (Navg - 1) to the left, average the samples in that window, and compute the MSE from the MMSE estimates obtained in the one-chain window to the precise estimate obtained by averaging the C' chains and the last Navg iterations. We show this error in terms of the window length and the noise variance for all increment distributions in Figure 9. The quality of the learned denoiser and the Monte Carlo denoiser differ over the noise variances and the learned denoiser improves relative to the Monte Carlo denoiser as the noise variance vanishes. Our final choice of S 300 samples is motivated by the fact that the quality of the Monte Carlo denoiser, when averuged across all noise variances that appear in the reverse diffusion, is always strictly better than the learned denoiser. Since it is relevant for the discussion in Section 4.2, we highlight that for this choice the quality of the Monte Carlo denoiser is superior to the learned one across all noise variances for the St(1) and (0,01.1 1) increment distributions.
F.5 ALGORITHM PARAMETERS
The adjustable regularization parameter for est E {(2,4) satisfies
Nyal
A rst.* arg min AEA
3 xx
(y val A) - xx,
(68)
Nval
k=1
There, A is the loglinear grid A (Ai, (2,,),,-00 A2, ................................ ANmo with
An = 10a+(m
(b-a)
mmb
(69)
with a = (-5) and b = 5. Since the model-based methods are very fast, we can use the relatively high Nnb 1000.
The adjustable hyperparameters of the DPS methods were found by
Neal
A alg.* arg min
Nyal 2 aal
val A) My
xx 122
(70)
A
where the grid Oalg is method-ddependen. This tuning is tailored to the evaluation with respect to the MMSE optimality gap. Due to resource constraints, the parameters are tuned for the learned denoiser. We use Nsemples 10 for the grid search on the validation set. We define modest number of Naps = 40 grid-points and found the extreme points of the grid (i.e., the values of the parameters that clearly lead to worse results) by hand. For C-DPS and DifiPIR, we fix the diffusion schedule to standard choices (3p = 1 x 10-4, Br = 0.02). In addition to the diffusion schedule, C-DPS has one tunable parameter y that we tune on 40 loglinear grid points (n 1, 1...., Náps)
10 a+ 11) Naps
(b-a)
(71)
There, a = (-3) and b = 1. DifiPIR has two tunable parameters y and ñ, with y being typically considered uncritical. Thus, we split the 40 grid points into a two-dimensional grid DDIFPIR (0.3,0.7) OS, with 2 points for y and 20 points for ( given by O
,r Nap where
c 10 an (idp.
(72)
with a = (-4) and b 1. The DPnP algorithm only has the schedule {n } to tune. In this case, since DPnP is asymptotically correct, the schedule is practical vehicle that enables to trade off between speed and accuracy. Therefore, the schedule of this paper is similar to the one that was proposed by (Xu & Chi, 2024): We fix a small Manll = 0.15 and linearly decrease 17 from some Tinitial to Vinnal after K/5 initial iterations with Tinitial, like
Omitial
if n 1....,K/5...
7h
K/G R73 Mnitiil
(73)
Mtinal
if n K/5 + 1,..., K
Mostiial
27

---

Published as a conference paper at ICLR 2026
BL(0.1,1)
Laplace(1) -
St(2)
St(3)
St(1)
Gauss(O, 0.25)
Denoising
Debluriig
Imputattion
Fourier
2
1.5
0.5
10
1000
1020
6
100
1010
10
10f0
10
109
2
1.5
0.5
0
10
10f0
10
10f0-
10
1010
10
10
1.5
CDDPS
1.0
10
10a
10-3
10
10
10
10
10
10-3
10
10*
2
1.5
D:NPR
1
0.5
0
10-2 2
100
4
10
10
10
10
0.5
10
10a
10a
10"
10
10a 10a
10a
10
10a
10a
10a
Figure 10: Grid-search diagnostics (logarithm of the MSE over the validation data set) for the model- based methods and the DPS algorithms. Rows: Ezj bi; C-DPS; DiffPIR; DPnP. Columns: Denoising; deconvolution; imputation; reconstruction from partial Fourier measurements. For better visualization, each curve has had its minimum subtracted. To limit clutter, marks are spaced ten apart.
We treat Tinitial tunable parameter and search over DPuP {mi U2,.... :700), where
Tn (0 (t-))
(74)
There, (-1) and b 4. Like in the original publication, use the comparatively small K 40.
The MSE over the validation data depending on the value the adjustable regularization parameter of the ls and & estimators and the adjustable hyperparameters of C-DPS, DiffPIR, and DPnP is shown in Figure 10. Since the y parameter of DifiPIR is assumed to be uncritical, we only show the values of the MSE for various choices of Ç, where is set the value of the optimal (7, Ç) pair.
28

---

Published as a conference paper at ICLR 2026
Table 4: Change in MMSE optimality gap (mean = standard deviation) after substituting the learned denoiser with the arbitrary- arbitrary-precision precision denoiser. An asterisk indicates a significant changes according to a Wilcoxon signed-rank test (p- 0.05). Negative number with asterisk: MMSE estimates obtained with the arbitrary-precisee
denoiser are significantly better. Positive number with asterisk: MMSE estimates obtained with the learned denoiser are significantly better.
Camss(0.25)
Laplace(1)
BL
St(1)
St(2)
Se(3)
C-DPS
8:00 + 0.11
0.00 + 0.16
-0.46 1.16"
0.00 + 0.01
-0.01
0.14
Denoising
DiMPI
8.00 0.13
0.00 0.17
0.05 0.78"
0.41 0.80"
0.00 0.20
0.00
0.115
DPaP
0.04 4: 0.27
0.01 4: 0.22
0.55 + 1.31"
0.77 4: 1.31"
0.00 0.24
0.00 0.23
C.DPS
0.24
0.00 0.26
0.97
00.12 1.11"
0.43
Deconvoluticn
DRIPIR
-.001 4: D:23
0.00 + 0:23
0.01 + 1.12
13.56 4 9.90
-0.01 + 0.47
0.00:0..00 0.31
DPal
0:25
0:27"
1:20
3.80"
0.06 0.77
-0.02 0.34
C-DPS
0:00 1:00+:0.30
0.01 + 0.35
d.41 # 1.51"
3.41 1.90
-0.12 1.01"
-0.01 = 0.57
Imputation
DiMPIR
0.00
0.00 0.33
0.03 1.05
3.05'
0.03 0.71
0.00--0. 0..47
DPaP
0.00 4: 0.35
0.02 4: 0.38
0.02 +: 1.02
10.46 4 5.70
0.02 4: 0.67
0.01 0.48
CDDPS
-0.02 1: 0:43
-0.01 4: 0.49
0.80 1.45"
0.09 4 5.63"
00.03 0.79'
0,01 = 0.49
Fourier
DiIPIR
-0.01
1: 0.39
0.00 1: 0.40
0.12 4: 0.83"
0.64
00.03 2: 0.42"
-0.02 0.33
DPaP
-0.01 + 0:43
0.00 + 0:45
-0.33 + 1.13"
1.32 3,18
0.00 0.54
0.01 0.46
Table 5: Runtime of the benchmark with learned objects.
Gauss((0.25)
Laplace(L)
BL(0.1.1)
St(1)
St(2)
St(3)
C-DPS
00:04:52
00:04:52
00:02:56
00:44:52
00:4::52
00:04522
Denoising
DimPIR
00:01:59
00:01:58
08:01:12
00:01:58
00:01:59
00:01:59
DPul
0040233
00:04:58
00:0:11
00:59:33
00:06:3
00:44588
C-DPS
00:04:52
00:04:53
00:02:57
00:04:53
00:04:53
00:04:52
Decomvdurion
DiiPIR
00:11:99
00:01599
00:01112
00:01:59
00:01:59
00:01:59
00:13:54
00:46:39
00:05:48
00:53:30
00:28:24
C-DPS
00:4::3
00:04:53
08:02:59
00:04:53
00:04:53
00:04:53
Imputation
DHIPIR
00:01:59
00:01:59
00:01:13
00:01:59
00:01:59
00:01:59
DPuP
00:4:58
00:16:18
18:18:56
00:51:41
00:39:04
00:32:50
C-DPS
00:04:54
00:04:54
00:02:59
00:01:55
00:01:55
00:01:54
Fourier
DiiPIR
00:01:59
00:01:59
08:01:13
00:01:59
00:01:59
00:01:59
DPuP
00:06:13
00:13:53
08:04:42
00:51:41
00:23:39
00:16:18
G ADDITIONAL RESULTS
We provide in Table an exhaustive quantitative evaluation of the change in the optimality gap after we substitute the report for which cases the cubitrary-precisiin denoiser enjoys significantly better results than the learned denoiser according to Wilcoxon signed-rank test (p 0.05, Ntest pairs, two-sided test with the winner determined by the median of differences). We attribute a better performance of the learned denoiser to the fact that the algorithms are fine-tuned using the learned component or to the cases where the likelihood score approximation is compensated by the one of the learned component. Note that this table must be interpreted with the quality of the denoisers in mind. As we show in Figure 9, for our particular choice of 5 = 300 samples, the Monte Carlo denoiser is strictly better than the learned denoiser over all noise variances only for signals with (0, 1) and St( increment distributions.
Monte Carlo denoiser for the learned denoiser. We also
We show uncurated qualitative results of the MMSE estimate obtained by the DPS algorithms and the gold-standard Gibbs methods in Figures 11 to 18. The figures alternate between the rebitrary-precision denoiser and the learned denoiser and show the results for deconvolution, denoising, imputation, and reconstruction from partial Fourier samples, in that order. Each figure contains results for (L...,, 1), St(1), (2), and Laplace(1 (1) increment distributions.
The coverage results for a = 0.9 are presented in Table 7. The Gibbs row again validates the implementation; for all forward operators, they achieve coverages that are very close to 0.9. In contrast, the coverage values obtained by the DPS algorithms are generally much smaller than 0.9. For C-DPS and DifiPIR, the reported coverage values are almost always 0 except for (L, 1) and St(1) increments, where the coverages usually (close to) 1 for C-DPS and inconsistent for DifiPIR. For almost all increment distributions and forward operators, DPnP reports coverage values that are closest to, but typically smaller than, 0.9. Note that coverage of 1 can be considered the worst case even at a target of 0.9. For instance, it would be achieved by setting all samples to a constant vector with extremely large (i.e., "unlikely") entries.
29

---

Published as a conference paper at ICLR 2026
DiffPIR
DPuP
BL.0.1
-2
2
0
0.6
0.8
0.6
0.4
0.2
0.2
200
8((1)
200
200
-106
100
-10
1.5
1.6
1
0.5
Mmman Whn
Mm
Mm
20
(c)
-5
-10
-20
-1
15
10
10
Laplace(()
-5
-10
-15
1.5
0.5
20
40
GD
20
40
20
40
60
Figure 11: Qualitative results for deconvolution using the Monte Carlo. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
30

---

Published as a conference paper at ICLR 2026
C-DPS
3-
DiMPIR
DPnP
BL(0...1)
2
0.8
0.6
100
50
St
200
10
m m
20
Ss(2)
16
2
Tep]ace(() ()
-10
-15
0.5
20
40
GD
20
40
20
40
60
Figure 12: Qualitative results for deconvolution using the learned denoiser. Rows: increment distribu- tions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gtad-standard Gibbs methods are shown on top of the corresponding index- wise marginal variances. Columns: Different measurements.
31

---

Published as a conference paper at ICLR 2026
C-DPS
DifPIR
DPuP
BL.(0.1.1)
Gibhs
-10
0.15
0.1
0.4
10
0.2
0
20
40
St(1)
20
-100
200
300
10
St
-5
-10
-2
0.4
0.4
0.3
0.2
0.1
0.1
10
Laplace(()
mmm
-6-
-10
-10
0.15
0.15
0.15
0.1
0.1
u.1
Monaro M Ao
20
40
60
20
40
60
Figure 13: Qualitative results for denoising using the Monte Carlo denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
32

---

Published as a conference paper at ICLR 2026
C-DPS
DifPIR
DPuP
BL.(0-1.1)
Gibhs
-10
0.4
0.1
0.2
10
0
20
40
St(1)
20
-100
200
300
10
10
10
20
St
-0
10
-10
-20
0.3
0.3
0.2
10
Laplace(1)
Mmm
Mym
-5-
-10
-10
0.15
0.15
VAAM
0.1
0.1
5-
10
10-22
20
40 60
20
40
60
Figure 14: Qualitative results for denoising using the learned denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
33

---

Published as a conference paper at ICLR 2026
C-DPS
DifIIIR
DPuP
Gibhs
1.5
0.5
400
St(1)
-200
1000
M000 AM M A
100 110d Madsn
2
20
S(
-10
10
10
10
6
4
Mww.mma
Ao
10
Laplace (1)
-10
16
An A Am
A 20
Mma
20 49
60
40
60
Figure 15: Qualitative results for imputation using the Monte Carlo denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
34

---

Published as a conference paper at ICLR 2026
C-DPS
DiffPRR DPuP
131.(0111) 1,
Gihbs
0.5
200
Std
500
200
-1,000
309
108
1.5
2
0.5
0.5
4.
0
St(2)
-10
-10
15
15
AMA Wh W MAmam
5
10
Laplace(() (1)
-5-
5
-10
-1
4
M AAA.A.M.
A
2
APW W.W.A.
20
40
80
Figure 16: Qualitative results for imputation using the learned denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
35

---

Published as a conference paper at ICLR 2026
C-DPS
DifIIIR
DPuP
BL((.111))
Gibhs
0.4
0.2
0.2
20
40
St(1)
20
-100
-20
200
300
200
40
100
20
50
20
St
-0
-10
m
-20
-15
10
Laplace (1)
-10
1.5
0.5
20
49
60
20 40
80
Figure 17: Qualitative results for reconstruction from partial Fourier measurements using the Monte Carlo denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
36

---

Published as a conference paper at ICLR 2026
C-DPS
DifPIR
DPuP
BL((-.11))
Gibhs
0
-10 -10-2
0.2
0.1
20
40
St(1)
20
-100
-20
2000
300
200
40
100
20
50
S(
-6
-10
-20
16
Da hy W
10
Laplace (1)
-5
-10
1.6
0.5
0.6
20
40
60
20
40
80
Figure 18: Qualitative results for reconstruction from partial Fourier measurements using the learned denoiser. Rows: increment distributions. For each increment distribution, the MMSE estimates obtained by the different DPS algorithms and the gold-standard Gibbs methods are shown on top of the corresponding index-wise marginal variances. Columns: Different measurements.
37

---

Published as a conference paper at ICLR 2026
Table 6: Runtime of the benchmark with Monte Carlo objects.
Gorsss(0.25)
Laplace(A)
DI,(0.1.1)
St(1)
St(2)
St(3)
C.DPS
05:52:28
07:23:23
34:07:44
05:52:40
05:34:51
05:31:10
Denoising
DIMPRR
05:04:40
06:36:25
33:46:29
05:12:30
05:99:06
05:38:24
DPul
00:03:0
00:03:57
00:20:36
00:03:23
00:03:08
00:03:10
C.DPS
07:25:17
34:17:12
05:28:38
05:24:12
05:24:00
Decomvlution
DMIPIR
05:28:09
06:55: 34
34:16:17
05:31:29
05:32:32
05:22:39
DPuP
00:0300
00:03:99
00:21
00:03:13
00:03:21
00:03:21
C-DPS
05:19:07
07:15:41
34:29:37
05:53:44
05:27:M
05:26:05
Imgutation
DiAPIR
05:50:15
07:00:13 13
33:52:26
05:34:00
05:24:16
05:09:56
DPal
00:03:25
00:04:15 18
(8) 30:30:58
00:05:09
00:03-055
00103 22
C-DPS
05:49:49
07:09:51
34:30:13
05:49:44
05:49:26
05:49:07
Fourier
DHIPIR
05:13:06
06:38:32
34:31:38
05:17:58
06:14:52
05:15:14
DPuP
00:03:04
00:04:12
00:20:59
00:03:05
00:03:19
00:03:32
= 999
8 600
8 200
t= 0
2
My M000
2
2
0
DIPPER XXMSE (y))
Gilbb (xr)
x
Ax+n
-4
-4
40
20
60
20
60
20
40
Figure 19: Conditional generation for deconvolution of a signal with (L.1. increments with Prototypical sampling trajectory at times 999, 600, 200, O. Bottom: From left to right: MMSE estimate obtained by averaging all DiffPIR samples; gold-standard MMSE estimate obtained by the Gibbs method; the data-gene signal; the data.
DiffPIR. Top:
Table 7: Posterior coverage of various estimation methods at 0.9. MC: Monte Carlo.
Gaues(0, 0,25)
Laplace(1)
BL/0.1.1
St(L)
St(3)
St(3)
Leared
Leamed
Leamed
Leamed
Gibbs
0.90
0.91
0.91
0.89
0.99
01.9
Denoising
C-DPS
0.00
0.00
0.30
0.00
1.00
1.00
1.00
1.00
0.00
0.00
0.00
0.00)
DifIPIR
0.00
0.02
0.00
DPaP
0.58
0.67
0.11
0.11
1.00
0.41
0.53
0.08
0.09
0.00
0.09
0.10
Gibbs
0.89
0.50
0.50
0.91
0.91
0.91
Deconvolation
C-DPS
0.00
0.00
DiiPIR
0.00
0.00
0.00
0.00
1.00
1,00
0.97
0.92
0.00
0.00
0.00
0.00
DPnP
0.12
0.12
0.06
0.07
1.00
0.31
0.50
0.06
0.06
0.06
0.07
0.06
Gibbs
Imputation
C-DPS
0.00
0.00
0.00
0.00
1.00
1.00
0.94
0.78
0.15
0.15
0.00
0.00
DiAPIR
0.00
0.00
0.00
0.00
1.00
1.00
0.72
0.32
0.00
0.00
0.00
0.00
DPnP
0.28
0.31
0.09
1.00
0.41
0.56
0.07
0.14
0.15
0.12
0.15
Gibbs
0.50
0.50
0.91
0.52
0.91
Fourier
C-DPS
0.00
0.00
0.00
0.00
1.00
1.00
0.96
0.74
0.01
0.01
0.00
0.00
DifiPIR
0.00
0.00
0.00
1.00
1,00
0.92
0.65
0.00
0.01
0.00
0.19
0.19
0.08
0.06
1.00
0.32
0:50
0.06
0.07
0.07
0.07
38


Reviews:
[{"review_id": "qh8Nh3DeU4", "reviewer": "Reviewer_tkeZ", "summary": "- The authors introduce a benchmark suite for evaluating algorithms designed to solve linear inverse problems with diffusion model priors\n- The benchmark is built on a synthetic setup derived from discretized Lévy processes\n- It hence include setting of heavy-tailed/power-law–like distributions beyond the Gaussian case\n- The key motivation lies in the fact that Lévy processes possess explicit marginal distributions and can be targeted using Gibbs sampling\n- This property allows the benchmark to generate ground-truth posterior samples (from inverse problem and denoising posterior) for quantitative comparison across algorithms", "questions": "- I generally found the figures hard to understand and interpret, I'm referring namely to figure 1, it says that it shows reverse using the oracle denoiser, but it is not clear, similarly, for figure 3, it is hard to interpret namely to say wether the algorithms performs well or not \n- can the authors provides hints/explanation on the derivation of equation (13)\n- the authors claim that introduced framework can also assess the approach where the conditional components is learned (Line: 168-170), but it is not clear how it can be achieved given that in some tasks the likelihood is not known, e.g. tasks such deraining or dehazing, see for instance [1]; yet the benchmark is built on the ability to explicitly write the posteriors/marginals and target them using Gibbs sampling\n\n- A more broad question: did the authors think about how the benchmark can be extend to nonlinear inverse problems ?\n\n---\n\n.. [1] Wang, Hanting, et al. \"IRBridge: Solving Image Restoration Bridge with Pre-trained Generative Diffusion Models.\" arXiv preprint arXiv:2505.24406 (2025).", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written and accompanied with concise explanations in the appendix\n- The motivation of the paper is well articulated namely for principled benchmarking in diffusion-based inverse problem solvers\n- The proposed benchmark is a valuable contribution, as it extends evaluation on Gaussian setup to a broader family of distributions", "weaknesses": "**Overstated or misleading claims**\nThe repeated use of the term \"oracle\", e.g., Lines 56, 129, 277, 355 is misleading.\nThe samples used in the benchmark are produced via Gibbs sampling—an iterative procedure—hence they are approximate, not exact. The quality of these samples depends on choices such as burn-in time, which are hyperparameters of the framework.\nThis issue becomes more apparent when the benchmark is applied to algorithms requiring gradients of the denoiser (Line 257-263 and equation (60)): the paper substitutes the latter with a covariance estimator of $X_0 | X_t,$ and hence further deviating from the notion of an \"oracle\".\n\n\n**Template for posterior samplers**\nThe proposed benchmark template seems overly restrictive. By focusing on algorithms that use only the denoiser, it neglects methods that require the Jacobian of the denoiser.\nAlthough the paper connects this to the covariance $Cov(X_0 \\mid X_t)$, estimating this covariance is far more computationally demanding and less stable, and therefore it downgrade the claim that the benchmark offers \"oracle\" quantities with minimal approximation error.\n\n**Evaluation design**\n- The inclusion of learned denoisers in the evaluation is conceptually inconsistent with the paper’s stated goal of removing approximation errors (Section 1.1).\nIf the benchmark aims to isolate algorithmic performance, learned denoisers reintroduce training-dependent variability. While the authors justify this by citing robustness testing, the notion of robustness is loose and in practice requires hyperparameter tuning, which introduces additional confounding factors.\n- The experimental comparison is limited. Only 3 algorithms are evaluated, and these do not represent the diversity of available approaches, e.g., optimization-based, variational, or midpoint-guided methods; see the literature in [1] and [2]\n\n**Remarks and minor issues**\n\n- In background, rephrase the statement in Line 132 about DDPM, sampling in fact depend on several parameters and it is bold to say \" researchers typically use\"; I would argue that frequently DDIM sampling is used with $\\eta = 0$ (simulating the probability-flow ODE) for sharp samples with few diffusion steps\n- The used abbreviation **DPS** is already/actually the name of a well-known algorithm in diffusion models and inverse problems [4], hence the abbreviation might be misleading using it here to refer to something else may cause confusion.\n- The authors may also consider adding the following reference on inverse problems benchmarks [3]\n- Line 288: The statement that DiffPIR is an extension of C-DPS is incorrect. DiffPIR follows a distinct formulation based on quadratic half-splitting with an auxiliary variable and does not rely on the VJP of the denoiser.\n\n\n---\n\n.. [1] Daras, Giannis, et al. \"A survey on diffusion models for inverse problems.\" arXiv preprint arXiv:2410.00083 (2024).\n\n.. [2] Oliviero-Durmus, Alain, et al. \"Generative modelling meets Bayesian inference: a new paradigm for inverse problems.\" Philosophical Transactions A 383.2299 (2025): 20240334.\n\n.. [3] Zheng, Hongkai, et al. \"Inversebench: Benchmarking plug-and-play diffusion priors for inverse problems in physical sciences.\" arXiv preprint arXiv:2503.11043 (2025).\n\n.. [4] Chung, Hyungjin, et al. \"Diffusion posterior sampling for general noisy inverse problems.\" arXiv preprint arXiv:2209.14687 (2022).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "Ecx5TSO86m", "reviewer": "Reviewer_bQ8j", "summary": "Diffusion posterior sampling algorithms have become prominent methods for sampling from posterior diffusion with a denoising diffusion model prior. While many methods have been proposed in the recent years, most of the interesting benchmarks do not come with ground-truth posterior samples to which one can compare against. The aim of this paper is to close this gap by proposing a statistical benchmark that mimicks the behaviour of realistic data (power-law-like extremes as stated in the paper). To this end the authors consider the posterior associated to Lévy processes and use an efficient Gibbs sampler to obtain gold-standard posterior samples that serve as reference.", "questions": "I have a few suggestions and related works to be considered: \n- I think it would have been interesting to include samples from a conditional diffusion model, by either training the conditional denoiser or estimating the denoiser using Monte Carlo samples as is done for DPS methods. I believe that it could be relevant since it provides a lower bound on the performance that one hopes to achieve with DPS methods. \n- [1] considers an actual real world setting where gold standard samples can be obtained using MCMC. \n- The toy Gaussian mixture benchmark is introduced in [2, 3] \n\n[1] Cardoso, G.V. and Pereira, M., 2025. Predictive posterior sampling from non-stationnary Gaussian process priors via Diffusion models with application to climate data.  \n[2] Cardoso, G., Idrissi, Y.J.E., Corff, S.L. and Moulines, E., 2023. Monte Carlo guided diffusion for Bayesian linear inverse problems.  \n[3] Boys, B., Girolami, M., Pidstrigach, J., Reich, S., Mosca, A. and Akyildiz, O.D., 2023. Tweedie moment projected diffusions for inverse problems.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "- This paper tackles a fundamental problem in the evaluation of diffusion posterior samplers and proposes a very useful benchmark which in my opinion could be useful to the community and should be present in all the forthcoming papers. \n- The model is general enough to contain different instantiations such as Laplace and spike and slab and thus goes beyond the existing gaussian mixture toy examples. \n- The paper is rather well-written and quite pedagogical, I enjoyed reading it.", "weaknesses": "The only weakness I see is the structuring of the main paper. For example I think that some parts of the related works (such as the first two paragraphs) could be moved to the appendix as they are slightly relevant to the content of the paper. This space could be used to provide for example more background on the GLM framework, as one needs to go to the appendix to read more interesting details about it. \nI also think that Figure 1 and 2 are misplaced as at this stage of the paper the Lévy process is not introduced and we don't know yet what St(1) means.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "XUACQbtL5B", "reviewer": "Reviewer_TmEt", "summary": "This paper introduces a statistical benchmark for evaluating diffusion posterior sampling algorithms using discretized Lévy processes with tractable Gibbs posteriors as ground truth. While the framework enables rigorous distribution-level comparisons, the evaluation is severely limited to low-dimensional (d=64) linear inverse problems, raising serious concerns about scalability and practical relevance to realistic imaging applications.", "questions": "Does this benchmark can be used for amortized diffusion sampling methods, i.e., learning the full posterior?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Developing a benchmark for posterior sampling in high-dimensional problems is important.", "weaknesses": "• All experiments use d=64 signals with only linear operators. No evidence is provided that the framework scales to realistic dimensions (e.g., 256×256 images) or nonlinear problems, fundamentally limiting the practical applicability and making it unclear whether findings transfer to problems researchers actually solve.\n\n• The authors cite power-law phenomena in finance and images to motivate heavy-tailed priors, but never demonstrate that their 1D discretized Lévy processes meaningfully capture structure in realistic signals. The connection to actual image statistics remains unsubstantiated.\n\n• Table 4 shows learned denoisers often match or exceed oracle performance, undermining claims about isolating likelihood approximation errors. The paper doesn't establish whether likelihood errors dominate versus other sources (discretization, hyperparameter sensitivity), weakening the diagnostic utility argument.\n\n• DPS algorithms are tuned with learned denoisers but evaluated with oracle denoisers using the same hyperparameters (lines 276-278). This mismatch means oracle results may be suboptimal, contradicting claims about properly isolating algorithmic errors.\n\n• Claims of \"efficient implementations\" and \"acceptable runtimes\" (lines 231-234, 822-823) lack any quantitative evidence; no runtime comparisons, memory usage, or scalability analysis is provided to substantiate efficiency claims or assess practical feasibility at higher dimensions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "5qBVdgq9kk", "reviewer": "Reviewer_pM9c", "summary": "The authors propose a new benchmark for evaluating different posterior sampling algorithms using diffusion models (dubbed DPS in this paper; Different from DPS [1]), where the posterior samples can be computed analytically, so that the ground truth is given. \nPrevious *benchmarks* that admit analytical posterior samples were constrained to settings where the prior is a mixture of Gaussians, which largely differs from the natural data statistics. The prior distributions considered in this paper is much larger, and the authors propose methods to efficiently compute ground truth posterior distributions. Several widely established baselines are compared.\n\n**References**\n\n[1] Chung et al., \"Diffusion posterior sampling for general noisy inverse problems\", ICLR 2023", "questions": "Is there any reason to constrain the benchmark for *diffusion* posterior sampling algorithms?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. To the best of my knowledge, this is the first approach to go beyond mixture of gaussian priors when attempting to build a ground truth posterior distribution.\n\n2. The paper is well-written and easy to follow, with sufficient background given in the appendix.\n\n3. The method of acquiring the posterior distribution by extending Kuric et al. [1] is sound.\n\n\n**References**\n\n[1] Kuric et al., \"The Gaussian latent machine: Efficient prior and posterior sampling for inverse problems\", arxiv 2025", "weaknesses": "1. Being able to use different prior/posterior distributions as ground truth is, in and of itself, important. Nevertheless, the argument would be strengthened if the paper shows that the proposed distributions in this paper are closer to real-world statistics in some cases. Currently, only some references are given.\n\n2. The authors mention that the proposed framework can be extended to higher-dimensional settings, but there are complications. It would add much value if the authors were to include experiments with $d$ that match the typical image resolutions. Currently, it seems like the experiments are conducted with low dimensionality ($d$). What's the value of $d$ chosen here?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)