Learning Diffusion Priors from Observations by Expectation Maximization
Authors: François Rozet
Gérôme Andry
Francois Lanusse
Gilles Louppe
Venue: NeurIPS 2024 poster (2024)
Keywords: diffusion
score-based
generative model
corrupted data
inverse problems
expectation-maximization
empirical bayes
Primary area: diffusion_based_models

Abstract: Diffusion models recently proved to be remarkable priors for Bayesian inverse problems. However, training these models typically requires access to large amounts of clean data, which could prove difficult in some settings. In this work, we present a novel method based on the expectation-maximization algorithm for training diffusion models from incomplete and noisy observations only. Unlike previous works, our method leads to proper diffusion models, which is crucial for downstream tasks. As part of our method, we propose and motivate an improved posterior sampling scheme for unconditional diffusion models. We present empirical evidence supporting the effectiveness of our method.

Learning Diffusion Priors from Observations by Expectation Maximization
François Rozet University of Liège
Gérôme Andry University of Liège gandrySuliege.. be
francois. rozetoulg....
François Lanusse Université Paris Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM francois lanusse@cnss.ff
Gilles Louppe University of Liège g louppe@uliege. be
Abstract
Diffusion models recently proved to be remarkable priors for Bayesian inverse prob- lems. However, training these models typically requires access large amounts of clean data, which could prove difficult in some settings. In this work, we present a novel method based on the expectation- expectation-maxxiizzttee maximizzation algorithm for training diffu- sion models from incomplete and noisy observations only. Unlike previous works, our method leads proper diffusion models, which is crucial for downstream tasks. As part of our method, we propose and motivate an improved posterior sam- pling scheme for unconditional diffusion models. We present empirical evidence supporting the effectiveness of our method.
Introduction
Many scientific applications can be formalized as Bayesian inference in latent variable models, where the target is the posterior distribution p(x y) x pyy x)) p(x:) given an observation y I8 Ad resulting from a forward process p(y z) and a prior distribution p(x) over the latent variable x E R Notable examples include gravitational lensing inversion [1-3], accelerated MRI [4-8], unfolding in particle physics [9, 10], and data assimilation [11-14]. In all of these examples, the observation alone is either too incomplete or too noisy to recover the latent æ. Additional knowledge the form of an informative prior p(z) is crucial for valuable inference.
Recently, diffusion models [15, 16] proved to be remarkable priors for Bayesian inference, demon- strating both quality and versatility [17-27]. However, to train a diffusion model for the latent æ, one would typically need large number of latent realizations, which by definition are or rarely accessible. This is notably the case in earth and space sciences where the systems of interest can only be probed superficially.
Empirical Bayes (EB) methods [28-31] offer a solution to the problem of prior specification in latent variable models when only observations y available. The objective of EB is to find the parameters of a prior model go(x) for which the evidence distribution go(y) = f p(y x:) go (x)) dæ is closest to the empirical distribution of observations p(y). Many EB methods have been proposed over the years, but they remain limited to low-dimensional settings [32-37] or simple parametric models [38, 39].
In this work, our goal is to diffusion models for the prior (o(x), they are best-in-class for modeling high-dimensional distributions and enable many downstream tasks, including Bayesian inference. This presents challenges for previous empirical Bayes methods which typically rely models for which the density (x) or samples x qo (x) are differentiable with respect to the
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

parameters 9. Instead, we propose an adaptation of the expectation canntttion-maimmzzttiiin [40-44] algorithm where we alternate between generating samples from the posterior qo (z) y) and training the prior go (3) on these samples. As part of our method, we propose an improved posterior sampling scheme for unconditional diffusion models, which we motivate theoretically and empirically.
2 Diffusion Models
The primary purpose of diffusion models (DMs) [15, 16], also known as score-based generative models [45, 46], is generate plausible data from a distribution p(z) of interest. Formally, adapting the continuous- Cantinuous-time formulation of Song et al. [46], samples E IR from p(x) are progressively perturbed through a diffusion process expressed as a stochastic differential equation (SDE)
dæ. = ft 2e dt + 9t dw;
(1)
where ft the drift coefficient, It E R the diffusion coefficient, Wt E V denotes standard Wiener process and z, E R. is the perturbed sample at time E [0, 1]. Because the SDE is linear with respect to zy, the perturbation kernel from æ to Te is Gaussian and takes the form
p(r, x) N (æt 04 Ex
where 02 and Σ, = I are derived from ft and gt [46-49]. Crucially, the forward SDE (1) has an associated family of reverse SDEs [46-49]
1 11 oe V., log dt + 1..t dw
da, = St 2t
where 77 2 is a parameter controlling stochasticity. other words, we can draw noise samples 21 ~ p(z1) x (0, ) and gradually remove the noise therein to obtain data samples xn p(xo) 2 p(x) by simulating Eq. (3) from t = to using an appropriate discretization scheme [16, 45, 46, 49-52]. In this work, adopt the variance exploding SDE [45] for which fi 0 and 0g 1.
In practice, the score function Vg,, log p(x,) in Eq. (3) unknown, but can be approximated by neural network trained via denoising score matching [53, 54]. Several equivalent parameterizations and objectives have been proposed for this task [16, 45, 46, 50-52]. this work, we adopt the denoiser parameterization dg (Tz and its objective [51]
arg @ (p(r))(()(rcrll)) >x || (zz, t) =112
for which the optimal denoiser is the mean B[r and] of p(1 wx). Importantly, B[x ard] is linked to the score function through Tweedie's ' formula [55-58]
E[z 21] = Iz + E, Ve, logp(zz).
(5)
which allows to use se(x;) = 23 (do (x xxe) as score estimate in Eq. (3).
3
Exppetation-Maxxiizzlli
The objective of the pnnectaaioonmaxiiizzttit (EM) algorithm [40-44] is to find the parameters 0 of a latent variable model go (z. y) that maximize the log-evidence log (v()) of an observation For
distribution of observations p(y), the objective is to maximize the expected log-evidence [43, 44] or,
equivalently, to minimize the Kullback-Leible (KL) divergence between p(y) and go(y). That is,
8" arg max Bp() log (v)
(6)
arg min KL (x(w)) go()))
(7)
0
The key idea behind the EM algorithm is that for any two sets of parameters 0, and 0,, we have qo_ (v) log
log qe, (x, v)
+ log go, (x y)
(8)
(x, y) Q0, (1) qo, (x,
qo. (x v)
log qa (25, y)
+ KT. (x 1)) go v))
(9)
Q0, (x., y)
2 Egg, (x1p)
log q0, (x, 8) - log qo, (x, v)]
(10)

---

This inequality also holds in expectation over p(y). Therefore, starting from arbitrary parameters dp. the EM update
Dk+1 arg max (p(()) gu, (41a) log qe (zz, v) log qa, (2,z)) y)]
(11)
= arg max Bp Bu, (=ly) [log qo (z. v)]
(12)
leads sequence of parameters 0% for which the expected log-evidence By(w) [log qa, (v)] monotonically increasing and converges to a local optimum [42-44].
is
When the expectation in Eq. (12) is intractable, many have proposed to use Monte Carlo approxima- tions instead [59-66]. Previous approaches include Markov chain Monte Carlo (MCMC) sampling, importance sampling, rejection sampling and their variations [63-66]. A perhaps surprising advantage of Monte Carlo EM (MCEM) algorithms is that they may be able to overcome local optimum traps [60, 61]. We refer the reader to Ruth [66] for recent review of MCEM algorithms.
4 Methods
Although rarely mentioned in the literature, the expectation expectation-maximins -axxmizzation algorithm is a possible solution to the empirical Bayes problem. Indeed, both have the same objective: minimizing the KL between the empirical distribution of observations p(y) and the evidence go(y). In the empirical Bayes setting, the forward model ply 2)) is known and only the parameters of the prior qo (x) should be optimized. In this case, Eq. (12) becomes
0
arg max Ppy()) Eye, (ala) log go (x) + log x)]
(13)
arg max E Epy) ge, (arle)
log qe (z)]
(14)
arg min KL TT (z) 11 go()))
(15)
where Tx (z) = f go. (x y) p(y) dy. Intuitively, (x) and therefore 9, (r) assign more density latents which are consistent with observations y < p(y) than go, (zr). this work, we consider special case of the empirical Bayes problem where each observation has an associated measurement matrix A and the forward process takes a linear Gaussian form pvy xx A) = NN Ax, Eg). This formulation allows the forward process to be potentially different for each observation y. For example, if the position or environment of sensor changes, the measurement matrix may also change, which leads to an empirical distribution of pairs (y. A) ~ ply, A). As result, (x(2)) in Eq. (15) becomes TH (x) qo, (x y, A) p(y, A) dy.
4.1 Pipeline
Now that our goals and assumptions are set, we present our method to learn a diffusion model ga (4) for the latent from observations y by expectation expectation-maxiiizzie maximization. The idea is to decompose Eq. (15) into (i) generating dataset of i.idd. samples from T& (2(z) and (ii) training 40 (x) to reproduce the generated dataset. We summarize the pipeline in Algorithms 1, 2 and 3. provided in Appendix A due to space constraints.
Expectation To draw from (t),), we first sample pair (v, p(y, A) and then generate 90, (x y, A) from the posterior. Unlike previous MCEM algorithms that rely on expensive and hard tune sampling strategies [63-66], the use of diffusion model enables efficient and embarrassingly parallelizable posterior sampling [21-23) However, the quality of posterior samples critical for the EM algorithm converge properly [63-66] and, in this regard, we find previous posterior sampling methods [21-23, 25, 26] to be unsatisfactory. Therefore, we propose an improved posterior sampling scheme, named moment matching posterior sampling (MMPS), which we present and motivate in Section 4.2. We evaluate MMPS independently from the context of learning from observations in Appendix E.
Maximization We parameterize our diffusion model go by denoiser network de (are () and train it via denoising score matching [53, 54], as presented in Section 2. To accelerate the training, we start each iteration with the previous parameters 0..
3

---

Initialization An important part of our pipeline is the initial prior qu(z). Any initial prior leads to local optimum [42-44], but an informed initial prior can reduce the number of iterations until convergence. In our experiments, we take Gaussian distribution Max MMx; Ex) as initial prior and fit its mean and covariance by - you guessed it! expectation- ceneratioommaiimzzationn Fitting a Gaussian distribution by EM is very fast as the maximization step can be conducted in closed-form, especially for low-rank covariance approximations [67].
An alternative we do not explore in this work would be to use a pre-trained diffusion model as initial prior. Pre-trainee treitrain can be contucted on data we expect to be similar to the latents, such as computer simulations or even video games. The more similar, the faster the EM algorithm converges. However, the initial prior qu (z())) does not cover latents that are otherwise plausible under the observations, the following priors 90, (z) may not cover these latents either. A conservative initial prior is therefore preferable for scientific applications.
4.2 Moment Matching Posterior Sampling
To sample from the posterior distribution qo(x v) x go(x)) ply 1)) of our diffusion prior qa (1)), we have to estimate the posterior score Vz, log (o / y) and plug it into the reverse SDE (3). In this section, we propose and motivate an improved approximation for the posterior score. As this contribution not bound to the context of EM, we temporarily switch back to the notations of Section 2 where our prior is denoted p(x) instead of (xx)
Diffusion posterior sampling Thanks to Bayes' rule, the posterior score Va, log p(1e y) can be decomposed into two terms [17, 18, 21-25, 46]
22, Logg((z; y) = Va, log p(az) + Vn, logp(y 21).
(16)
As an estimate of the prior score V,, logp(x;)) is already available via the denoiser (I,, t), the remaining task to estimate the likelihood score Vz, logp(y st). Assuming a differentiable
measurement function A and Gaussian forward process ply x)) No A(z), Ey); Chung et al. [21] propose the approximation
ply xx) / p(y w)) p(a x.) de N (y A(E[(> xe)),
(17)
which allows to estimate the likelihood score Vz, log logp(y 22)) without training any other network than do t) Rs E[x zz].- The motivation behind Eq. (17) is that, when og is small, assuming that p(æ: | 22) is narrowly concentrated around its mean E[x 2x] is reasonable. However, this approximation is very poor when a, is not negligible. Consequently, DPS [21] is unstable, does not properly cover the support of the posterior p(r y) and often leads samples which inconsistent with the observation [22-25].
Moment matching To address these flaws, following studies [22-25] take the covariance V[r ar.] into account when estimating the likelihood score Vz, log p(y 21). Specifically, they consider the Gaussian approximation
q(x Ix)) N (a E[r xz). V[x xx])
(18)
which is closest to p(x ax) in Kullback-Leibler (KI.) divergence [68]. Then, assuming linear Gaussian forward process p(y x) = N2y Ar, E,), we obtain [68]
qlv xx)) = ply w) q x.)) dx = N (y AEE xel, + AV(x xx)AT)
(19)
which allows to estimate the likelihood score Va, log p(y xx2) as
Vz, q (x)) V.,, E[æ wr) TAT(E,) AV [ =TJATT- (y ABix xx))
(20)
under the assumption that the derivative of V[x [z]] with respect to x; is negligible [24, 25J. Even with simple heuristics for V(x wx., such as Σ, [20] or (E, + Σー [22, 23], this adaptation leads to significantly more stable sampling and better coverage of the posterior p(x y) than DPS [21]. However, we find that heuristics lead to overly dispersed posteriors g(z; y) x p(x:) g(y xx)
in the presence of local covariances L.e. covariances in the neighborhood of Z4.

---

A
Viaiie
10a
y
Az
Sx(p(m a/) || g u))
10a
41
p(a)
B
10
y = Ar
10-3
24
p(c)
C
10
y Az
St
10
109 10'
p(x)
Figure 1. Illustration of the posterior g( v)
Figure 2. Sinkhorn divergence [69] between the posteriors p(x, y) and g( y) for different heuristics of V [ xta] when the prior p(x) lies on 1-d manifolds embedded RR3. Lines and shades represent the 25-50-75 percentiles for 64 randomly generated manifolds [71] and mea- surement matrices A E IR 1 Using V| xa] instead of heuristics leads to orders of magni- tude more accurate posteriors g( y)
for the Gaussian approximation ((x: xr) when the prior p(z) lies on a manifold. Ellipses repre- sent 95 % credible regions of q(x zx). (A) With E, as heuristic for V(x 2z). any whose mean E[x zz] is close to the plane = Ar is consid- ered likely. (B) With V(z zz]. more regions are correctly pruned. (C) Ground-truth p(r, y) and p(æ) (2) for reference.
We illustrate this behavior in Figure 1 and measure its impact as the Sinkhorn divergence [69, 70] between the posteriors p(z, y) and q(2, y) when the prior p(x) lies randomly generated 1-dimensional manifolds [71] embedded in TR3 The prior p(x:) is modeled as a mixture of isotropic Gaussians centered around points of the manifold, which gives access to p(æ.). ), E[x: xea] and V[x aa)] analytically. The results, presented in Figure 2, indicate that using V[x st instead of heuristics leads orders of magnitude more accurate posteriors (xx, y). We expect this gap further increase with real high-dimen nigh-dimensional data as the latter often lies along low-dimensional manifolds and, therefore, presents strong local covariances.
Because the MCEM algorithm is sensitive to the accuracy of posterior samples [63-66], we choose to estimate V[x my] using Tweedie's covariance formula [55-58]
Vlz 21] E; Σn V log p( (I.) Σx Σ.
(21)
E Jzl
V do
(22)
Conjugate gradient method As noted by Finzi et al. [24]. explicitly computing and materializing the Jacobian V do N is intractable in high dimension. Furthermore, even if we had
access to V[x Ir), naively computing the inverse of the matrix Σ, + AV[r xx]AA T in Eq. (20)
would still be intractable. Fortunately, observe that the matrix E,y AV[x xx xx]ATT T is symmetric positive definite (SPD) and, therefore, compatible with the conjugate gradient (CG) method [72].
The CG method is an iterative algorithm to solve linear systems of the form Mr b where the SPD matrix M and the vector b are known. Importantly, the CG method only requires implicit access to M through an operator that performs the matrix-vector product Me given vector v. our case, the linear system to solve is
AE (1
x. = (Ey + AV x.)AA)
(23)
= =Ey + A AAE E[x x
(24)
vector-Jasobian product
Within automatic differentiation frameworks [73, 74], the vector-Jacobian product in the right-hand side can be cheaply evaluated. In practice, due to numerical errors and imperfect training, the Jacobian
5

---

zz] is not always perfectly SPD. Consequently, the CG method becomes unstable after a number of iterations and fails to reach an exact solution. Fortunately, we find that truncating the CG algorithm very few iterations (1 to 3) already leads to significant improvements over using heuristics for the covariance V/x zi)- Alternatively, the CG method can be replaced by other iterative algorithms that can solve non-symmetrie non-definite linear systems, such as GMRES [75] or BiCGSTAB [76], at the cost slower convergence.
Vo,, = E/r V
Results
We conduct three experiments to demonstrate the effectiveness of our method. We design the first experiment around low-dimensional latent variable x whose ground-truth distribution p(x) is known. In this setting, can use asymptotically exact sampling schemes such predictor-correctn sampling [23, 46] or twisted diffusion sampling [77] without worrying about their computational cost. This allows us validate our expectation maximization pipeline (see Algorithm 1) the limit of (almost) exact posterior sampling. The remaining experiments target two benchmarks from previous studies: corrupted CIFAR-10 10 and accelerated MRI. These tasks provide good understanding of how our method would perform in typical empirical Bayes applications with limited data and compute.
5.1 Low-dimensional manifold
In this experiment, the latent variable x I2 p(z) lies on a random I-dimensional manifold embedded in R represented in Figure 7. Each observation y E R² NN Axx Σy) is the result of a random linear projection of a latent x plus isotropic Gaussian noise (Sy = 10-47). The rows the measurement matrix 22x5 are drawn uniformly from the unit sphere S4, note that observing all push-forward distributions p(u x)) with
S of distribution p(z) in RN is sufficient to recover p(z) in theory [78, 79]. In practice, we generate a finite training set of 216 pairs (y. A).
We train DM (x())) parameterized by multi-layer perceptron do (az,; for K = 32 EM iterations. We apply Algorithm 3 to estimate the posterior score V., log (x; y: A), but rely on the predictor- corrector [23, 46] sampling scheme with large number (4096) of correction steps to sample from the posterior (x y, A). We provide additional details such as noise schedule, network architectures, and learning rate in Appendix C.
As expected, the model go, (x)) converges towards a stationary distribution whose marginals are close to the marginals of the ground- p(x), as illustrated in Figure 3. We attribute the remaining artifacts to finite data and inaccuracies in our sampling scheme.
gl (20, an)
0-
2-
Ue (11 .22)
2-
0-
2
(2(22, 12))
0
-2
2
-2
-2
-2 2
-2
0
-2
2
1:
1:
1:
32
Truth
Figure 3. Illustration of 2-d marginals of the model go, (x) along the EM iterations. The initial Gaussian prior go (Ir) leads to very dispersed first model ge, (I). The EM algorithm gradually prunes the density regions which are inconsistent with observations, until it reaches stationary distribution. The marginals of the final distribution are close to the marginals of the ground-truth distribution.

---

Method
FID
IS
0.20
11.70
7.97
27
AmbieniDifllsion [80]
0.40
18.85
7.45
0.60
28.88
6.88
24
0,25
5.88
8.83
Ours w/ Tweedie
0.50
6.76
8.75
13.18
8.14
Tweedle (23%)
Tweedie (50%)
Ours w/ (I + 2r 0.75
39.94
7.69
Tweedle (78%)
2
118.69
4.23
(1 + X, (75%)
M. (73%)
Table 1. Evaluation of final models trained on corrupted CIFAR-10. Our method outperforms AmbientDiffusion [80] similar corruption levels. Using heuristics for V/æ zy] instead of Tweedie's formula greatly decreases the sample quality.
25
30
Figure 4. FID of go, (z) along the EM iterations for the corrupted CIFAR- 10 experiment.
Figure 5. Example of samples from the model go, (x) along the EM iterations for the corrupted CIFAR-10 10 experiment with 0.75. We use the deterministic DDIM [50] sampling scheme for easier comparison. Generated images become gradually more detailed and less noisy.
5.2 Corrupted CIFAR-10
Following Daras et al. [80], we take the 000 training images of the CIFAR-IO [81] dataset as latent variables æ. A single observation y is generated for each image by randomly deleting pixels with probability p. The measurement matrix is therefore a binary diagonal matrix. We add negligible isotropic Gaussian noise (E, 10-10-0 for fair comparison with AmbientDiffusion [80] which cannot handle noisy observations.
For each corruption rate p E (0.25, 0.5,0.75), 0.75). we train a DM qo(x)) parameterized by U-Net [82] inspired network dg [z,, 1) for
32 EM iterations. We apply Algorithm 2 with 256 discretization steps and 77 = 1 to approximately sample from the posterior ge (x A.. We apply Algorithm 3 with several heuristics for V[z [x]] to compare their results against Tweedie's covariance formula. For the latter, truncate the conjugate gradient method in Algorithm 4 to single iteration.
For each model qo (2), we generate set of 50 000 images and evaluate its Inception score (IS) [83] and Fréchet Inception distance (FID) [84] against the uncorrupted training set of CIFAR-10. We report the results Table and Figures and 5. At 75 % of corruption, our method performs similarly to AmbientDiffusion [80] at only 20 % of corruption. On the contrary, using heuristics for V[x .r leads poor sample quality.
5.3 Accelerated MRI
Magnetic resonance imaging (MRI) is non-invasive medical imaging technique used in radiology to inspect the internal anatomy and physiology of the body. MRI measurements of object are obtained in the frequeney domain, also called k-space, using strong magnetic fields. However, measuring the

---

Mask
Zero-filld
Sample
Sample Ground-truth
R R=s
R= 16
==
Figure 6. Examples of posterior samples for accelerated MRI using diffusson prior trained from space observations only. Posterior samples detailed and present plausible variations, while remaining consistent with the observation. We provide the zero-filled inverse, where missing frequencies are set to zero, as baseline.
entire k-space can time-consuming and expensive. Accelerated MRI [4-8] consists in inferring the scanned object based on partial, possibly randomized and noisy, frequency measurements.
In this experiment, following Kawar al. [85], we consider the single-coil knee MRI scans from the fastMRI [7, 8] dataset. We treat each slice between the 10th and 40th of each scan as an independent latent variable æ, represented
320 320 gray-scale image. Scans min-max normalized such that pixel values range between -2 and 2. A single observation is generated for each slice by first applying the discrete Fourier transform and then a random horizontal frequency sub-sampling with acceleration factor R = [85, 86]. meaning that a proportion 1/R of all frequencies are observed on average. Eventually, obtain 24 853 k-space observations which we add isotropic Gaussian noise (Ey = 10-4T) 7) to match Kawar et al. [85].
Once again, train DM (I: parameterized by U-Net [82] inspired network dg t) for K = 16 EM iterations. We apply Algorithm 2 with T = 64 discretization steps and r 1 to approximately sample from the posterior qe (ar y. A) and truncate the conjugate gradient method in Algorithm 4 to 3 iterations. After training, we employ our final model gex (x) as a diffusion prior for accelerated MRI. Thanks our moment matching posterior sampling. scans at accelcration factors up to R 32, as shown in Figure 6. Our samples are noticeably more detailed than the ones of Kawar et [85]. We choose not to report the PSNR/SSIM of our samples as these metrics only make sense for regression tasks and unfairly penalize proper generative models [87, 88). We provide prior samples in Figure 13 and posterior samples for another kind of forward process in Figure 14.
able to infer plausible
Related Work
Empirical Bayes A number of previous studies have investigated the use of deep learning to solve the empirical Bayes problem. Louppe et al. [35] use adversarial training for learning prior

---

distribution that reproduces the empirical distribution of observations when pushed through non- differentiable black-box forward process. Dockhorn al. [33] use normalizing flows [89, 90] to estimate the prior density by variational inference when the forward process consists of additive noise. Vandegar et al. [36] also use normalizing flows but consider black-box forward processes for which the likelihood p(v x)) is intractable. They note that empirical Bayes is an ill-posed problem in that distinct prior distributions may result in the same distribution over observations. Vetter et al. [37] address this issue by targeting the prior distribution of maximum entropy while minimizing the sliced- Wasseesesn distance [78, 79] with the empirical distribution of observations. These methods rely on generative models go for which the density go(x) samples (0(x) are differentiable with respect to the parameters 6, which is not or hardly the case for diffusion models.
Closer to this work. Daras et al. [80] and Kawar et al. [85] attempt to train DMs from linear observations only. Daras et al. [80] consider noiseless observations of the form y = Ar and train a network do (Ara, A, t) to approximate R[[x| Arz] under the assumption that E[A 'A] is full-rank. The authors argue that TE Ar.]] can act as a surrogate for E(r ]r)]- Similarly, Kawar et al. [85] assume Gaussian observations y ~ Nyy Axx Mny) and train a network de (r,,, t) to approximate K[r Pxa] under the assumption that [[ is full-rank where P = A+A and A+ is the Moore-Penrose pseudo- inverse of A. The authors assume that dø(Px;; 1) can generalize to P I, even if the training data does not contain P = I. In both cases, the trained networks are not proper denoisers approximating E[r 22] and cannot reliably parameterize a standard diffusion model, which is problematic for downstream applications. Notably, in the case of Bayesian inference, they require custom posterior sampling schemes such as the one proposed by Aali et al. [91] for AmbientDiffusion [80] models. Conversely, in this work, we do not make modifications to the denoising score matching objective [53, 54], which guarantees a proper DM that is compatible with any posterior sampling scheme at every iteration. In addition, we find that our method leads to quantitatively and qualitatively better diffusion priors.
In a concurrent work, Daras et al. [92] propose an algorithm to train DMs from noisy (A = I and b, = o2/1) data by enforcing the "consistency" of the denoiser along diffusion paths. They prove that the mean E[x æx)] is the unique consistent denoiser. Interestingly, this training algorithm also relics on posterior samples, which are easy to obtain thanks to the white noise assumption.
Posterior sampling Recently, there has been much work on conditional generation using uncon- ditional diffusion models, most of which adopt the posterior score decomposition in Eq. (16). As covered in Section 4.2, Chung et al. [21] propose an analytical approximation for the likelihood score Vz, log logp(y rz) when the forward process p(y I) is Gaussian. Song et al. [22] and Rozet et al. [23] improve this approximation by taking the covariance V[æ zz] into account in the form of simple heuristics. We build upon this idea and replace heuristics with a proper estimate of the covariance V(x xy] based on Tweedie's 's covariance formula [55-58]. Finzi et al. [24] take the same approach, but materialize the matrix AV [x: 2: A which is intractable in high dimension. Boys ct al. [25] replace the covariance V[x xt] with a row-sum approximation diag(e T V|z æt]) where is the all-ones vector. This approximation is only valid when V(z 23] is diagonal, which limits its applicability. Instead, we take advantage of the conjugate gradient method [72] to avoid materializing AVEx x2. AT. A potential cheaper solution is to train sparse approximation of V/x wz). as proposed by Peng et al. [93], but this approach is less general and not immediately applicable to any diffusion model.
parallel line of work [94-96] extends the conditioning of diffusion models to arbitrary loss terms ((x,1) v) x - -log I). for which no reliable analytical approximation of the likelihood score Va, log p(y 2z) exists. Song et al. [94] rely on Monte Carlo approximations of the likelihood ply I x²) =f ply I) p(x 2;) dæ by sampling from Gaussian approximation of p(æ: | xx). Conversely, He et al. [96] use the mean E xz] as a point estimate for p(x xz).. but leverage a pre-trained encoder-decoder pair to project the updates of x, within its manifold. We note that our use of the covariance V[x ær]] similarly leads to updates tangent to the manifold 24-
Finally, Wu et al. [77] propose a particle-basee based posterior sampling scheme that is asymptotically exact in the limit of infinitely many partieles as long as the likelihood approximation g(y) xx.) here named the twisting function converges to p(y xx) as approaches O. Using TDS [77] as part of our
expectation-maxxinzaaooo pipeline could lead to better results and/or faster convergence, at the cost
of computational resources. In addition, the authors note that the efficiency of TDS [77] depends

---

on how closely the twisting function approximates the exact likelihood. In this regard, our moment matching Gaussian approximation in Eq. (19) could be a good twisting candidate.
7 Discussion
To the best of our knowledge, we are the first to formalize the training of diffusion models from corrupted observations as an empirical Bayes [28-31] problem. In this work, we limit our analysis to linear Gaussian forward processes to take advantage of accurate and efficient high-dimensional posterior sampling schemes. This contrasts with typical empirical Bayes methods which target low-dimensional latent spaces and highly non-linear forward processes [33-37]. In addition, as mentionned in Section 6, these EB methods are not applicable to diffusion models. As such, we choose to benchmark stress that a proper comparison with previous empirical Bayes methods would be valuable for both communities. We also note that Monte Carlo EM [59-66] can handle arbitrary forward processes p(y w) as long as one can sample from the posterior qo y). Therefore, our method could be adapted to any kind of forward processes in the future. We believe that the works of Dhariwal et al. [97] and Ho et al. [98] on diffusion guidance are good avenues for adapting our method to non-linear, non-differentiable, or even black-box forward processes.
work against similar methods in the diffusion model literature [80, 85], but
From a computational perspective, the iterative nature of our pnnecctatio-mmxxmmiaae method is drawback compared to previous works [80, 85]. Notably, generating enough samples from the posterior can be expensive, although embarrassingly parallelizable. In addition, even though the architecture and training of the model (a(x) itself are simpler than in previous works [80, 85], the sampling step adds a significant amount of complexity, especially as the convergence of our method is sensitive to the quality of posterior samples. In fact, we find that previous posterior sampling methods [21-23, 25, 26] lead to disappointing results, which motivates us to develop a better one.
As such, moment matching posterior sampling (MMPS) is byproduct of our work. However, it is not bound to the context of learning from observations and is applicable to any linear inverse problem given a pre-trained diffusion prior. In Appendix E, we evaluate MMPS against previous posterior sampling methods for several linear inverse problems on the FFHQ [99] dataset. We find that MMPS consistently outperforms previous methods, both qualitatively and quantitatively. MMPS is remarkably stable and requires fewer sampling steps to generate qualitative samples, which largely makes up for its slightly higher step cost.
Finally, as mentioned in Section 6, empirical Bayes an ill-posed problem in that distinct prior distributions may result in the same distribution over observations. In other words, it is generally impossible to identify "the" ground tround-truth distribution p(x) given an empirical distribution of obser- vations ply). Instead, for sufficiently expressive diffusion model, our EM method will eventually converge to a prior go (z) that is consistent with p(y), but generally different from p(r). Following the maximum entropy principle, as advocated by Vetter et al. [37], is left to future work.
Acknowledgments and Disclosure of Funding
François Rozet and Gérôme Andry are research fellows of the F.R.S.-FNRS (Belgium) and acknowl- edge its financial support.
The present research benefited from computational resources made available on Lucia, the Tier- 1 supercomputer of the Walloon Region, infrastructure funded by the Walloon Region under the grant 101910247. The computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRR) under the grant n22.5020.1 11 and by the Walloon Region.
MRI data used in the preparation of this article were obtained from the NYU fastMRI Initiative database 17, 8). As such, NYU fastMRI investigators provided data but did not participate in analysis or writing of this report. A listing of NYU fastMRI investigators, subject to updates, can be found https: : //ff/fff /fasstrii.../// med. nyu. edu/. The primary goal of fastMRI to test whether machine learning can aid in the reconstruction of medical images.
10

---

References
[1] S. J. Warren and S. Dye. "Semilincar Gravitational Lens Inversion". The Astrophysical Journal (2003).
[2] Warren R. Morningstar et al. "Data-driven Reconstruction of Gravitationally Lensed Galaxies Using Recurrent Inference Machines". In The Astrophysical Journal (2019).
[3] Siddharth Hishra-Sharma and Ge Yang. "Strong Lensing Source Reconstruction Using Continuous Neural Fields". 2022..
[4] Shanshan Wang et al. "Accelerating magnetic resonance imaging via deep learning". In International Symposium on Biomedical Imaging. 2016.
[5] Kerstin Hammernik et al. "Learning a variational network for reconstruction of accelerated MRI data". In Magnetic Resonance in Medicine (2018).
[6] Yoseo Han al. "k-Space Deep Learning for Accelerated MRI". In Transactions on Medical Imaging (2020).
[7] Jure Zbontar al. "fastMRR: An Open Dataset and Benchmarks for Accelerated MRI". 2018. Florian Knoll et al. "fastMRI: Publicly Available Raw k-Space and DICOM Dataset of Knee Images for Accelerated MR Image Reconstruction Using Machine Learning". In Radiology: Artificial Intelligence (2020).
[8]
[9] G. Cowan. "A survey of unfolding methods for particle physics". In Conf. Proc. C (2002).
[10] Volker Blobel. "Unfolding Methods in Particle Physics". In PHYSTAT. CERN, 2011.
[11]
Saangois-Xaviees Le Dimet and Olivier Talagrand. "Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects". In Tellus A: Dynamic Meteorology and Oceanography (1986).
[12] Yannick Trémolet. " AAcounting for an imperfect model in 4D-Var". In Quarterly Journal of the Royal Meteorological Society (2006).
[13] Thomas M. Hamill. "Ensemble-baue atmospheric data assimilation". In Predictability of Weather and Climate. 2006.
[14] Alberto Carrassi al. "Data assimilation in the geosciences: An overview of methods, issues, and perspectives". In WIREs Climate Change (2018).
[15] Jascha Sohl-Dickstein et al. "Deep Unsupervised Learning using Nonequilibrium Thermo- dynamics". In Proceedings of the 32nd International Conference on Machine Learning. 2015.
[16] Jonathan Ho et al. "Denoising Diffusion Probabilistic Models". In Advances in Neural Information Processing Systems. 2020.
[17] Yang Song et al. "Solving Inverse Problems in Medical Imaging with Score- -osee- Generative Models", International Conference on Learning Representations. 2022.
[18] Bahjat Kawar et al. "SNIPS: Solving Noisy Inverse Problems Stochastically". In Advances in Neural Information Processing Systems. 2021.
[19] Bahjat Kawar
"Denoising Diffusion Restoration Models". In Advances in Neural
Information Processing Systems. 2022.
[20] Alexandre Adam et al. "Posterior samples of source galaxies in strong gravitational lenses with score-based priors". 2022.
[21]
Hyungjin Chung et al. "Diffusion Posterior Sampling for General Noisy Inverse Problems". In International Conference on Learning Representations. 2023.
[22]
Jiaming Song et al. "Pseudoinverse -uineed Diffusion Models for Inverse Problems". In International Conference on Learning Representations. 2023,
[23] François Rozet and Gilles Louppe. "Score-based Data Assimilation". In Advances in Neural Information Processing Systems. 2023.
[24] Mare Anton Finzi al. "User-defined Event Sampling and Uncertainty Quantification in Diffusion Models for Physical Dynamical Systems". In Proceedings of the 40th International Conference on Machine Learning. 2023.
[25]
Benjamin Boys et al. "Tweedie Moment Projected Diffusions For Inverse Problems". 22023.
[26] Y. Zhu et al. "Denoising Diffusion Models for Plug-and-Pia Image Restoration". In Confer- ence on Computer Vision and Pattern Recognition Workshops. 2023.

---

Noe Dia et al. "Bayesian Imaging for Radio Interferometry with Score-based Priors". 2023. [28] Herbert E. Robbins. "An Empirical Bayes Approach Statistics". In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability. 1956.
[27]
[29] George Casella. "An Introduction to Empirical Bayes Data Analysis". In The American Statistician (1985).
[30] Bradley P. Carlin and Thomas A. Louis. "Empirical Bayes: Past, Present and Future". In
Journal of the American Statistical Association (2000).
[31] Bradley Efron. "Two Modeling Strategies for Empirical Bayes Estimation". In Statistical Science (2014).
[32] G. D' Agostini. "A multidimensional unfolding method based on Bayes' theorem". In Nuclear Instruments and Methods in Physics Research (1995).
[33]
Tim Dockhorn al. "Density Deconvolution with Normalizing Flows". 2020.
[34] Anders Andreassen al. "OmniFold: A Method to Simultaneously Unfold All Observables". In Physical Review Letters (2020).
[35]
Gilles Louppe al. "Adversarial Variational Optimization of Von-Differentiable Simulators". In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics. 2019.
[36] Maxime Vandegar et al. "Neural Empirical Bayes: Source Distribution Estimation and its Applications to Simulation-Bast Inference". In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics. 2021.
[37] Julius Vetter et al. "Sourcerer: Sample-based Maximum Entropy Source Distribution Estima-
tion". 2024.
[38]
Bradley Efron. "Empirical Bayes deconvolution estimates". In Biometrika (2016).
[39]
Balasubramanian Narasimhan and Bradley Efron. "deconvolveR: A G-Modeling Program for Deconvolution and Empirical Bayes Estimation". In Journal of Statistical Software (2020).
[40] H. O. Hartley. "Maximum Likelihood Estimation from Incomplete Data". In Biometrics (1958).
[41]
A. Dempster et al. "Maximum Likelihood from Incomplete Data Via the EM Algorithm".
In Journal of the Royal Statistical Society (1977).
[42]
C. F. Jeff Wu. "On the Convergence Properties of the EM Algorithm". In The Annals of Statistics (1983).
[43]
Geoffrey McLachlan and Thriyambakam Krishnan. "The EM algorithm and extensions". John Wiley & Sons, 2007.
[44]
Sivaraman Balakrishnan al. "Statistical guarantees for the EM algorithm: From population to sample-based analysis". .. In The Annals of Statistics (2017).
[45]
Yang Song and Stefano Ermon. "Generative Modeling by Estimating Gradients of the Data Distribution". In Advances in Neural Information Processing Systems. 2019.
[46] Yang Song et al. "Score-Basse Based Generative Modeling through Stochastic Differential Equa- tions". In International Conference on Learning Representations. 2021.
[47] Brian D. O. Anderson. "Reverse-tine diffusion equation models". In Stochastie Processes and their Applications (1982).
[48] Simo Särkkä and Arno Solin. "Applied Stochastic Differential Equations". Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2019.
[49] Qinsheng Zhang and Yongxin Chen. "Fast Sampling of Diffusion Models with Exponential Integrator", In International Conference on Learning Representations: 2023.
[50] Jiaming Song al. "Denoising Diffusion Implicit Models". In International Conference Learning Representations. 2021.
[51] Tero Karras et al. "Elucidating the Design Space of Diffusion- Based Generative Models". In
Advances in Neural Information Processing Systems. 2022.
[52] Yaron Lipman et al. "Flow Matching for Generative Modeling". In International Conference Learming Representations. 2023.
[53] Aapo Hyvärinen. "Estimation Non-Normale Statistical Models by Score Matching". In Journal of Machine Learning Research (2005).
12

---

[54] Pascal Vincent. "A Connection Between Score Matching and Denoising Autoencoders". In Neural Computation (2011).
[55] M. C. K. Tweedie. "Functions of a statistical variate with given means, with special reference to Laplacian distributions". In Mathematical Proceedings of the Cambridge Philosophical Society (1947).
[56] Bradley Efron. "Tweedie's Formula and Selection Bias". In Journal of the American Statistical Association (2011).
[57]
Kwanyoung Kim and Jong Chul Ye. "Noise2Score: Tweedie's Approach to Self-Supervised Image Denoising without Clean Images". In Advances in Neural Information Processing Systems. 2021.
[58]
Chenlin Meng et al. "Estimating High Order Gradients of the Data Distribution by Denoising". In Advances in Neural Information Processing Systems. 2021.
[59] Greg C. G. Wei and Martin A. Tanner. "A Monte Carlo Implementation of the EM Algorithm and the Poor Man's Data Augmentation Algorithms", In Journal of the American Statistical Association (1990).
[60] Gilles Celeux and Jean Dicbolt. "A stochastic approximation type EM algorithm for the mixture problem". In Stochastics and Stochastic Reports (1992).
[61] Bernard Delyon et al. "Convergence of a stochastic approximation version of the EM algo- rithm".. In The Annals of Statistics (1999).
[62] James G. Booth and James P. Hobert. "Maximizing Generalized Linear Mixed Model Likeli- hoods with Automated Monte Carlo EM Algorithm". In Journal of the Royal Statistical Society (1999).
[63] Richard A. Levine and George Casella. "Implementations of the Monte Carlo EM Algorithm". In Journal of Computational and Graphical Statistics (2001).
[64]
Brian S. Caffo et al. "Ascent-Based Monte Carlo Expectaton--Mxiiiiaatoori In Journal of the Royal Statistical Society (2005).
[65]
Wolfgang Jank. "The EM Algorithm, Its Randomized Implementation and Global Optimiza- tion". In Perspectives in Operations Research. 2006.
[66]
William Ruth. "A review of Monte Carlo-based versions of the EM algorithm". 2024. Michael E. Tipping and Christopher M. Bishop. "Mixtures of Probabilistic Principal Compo- nent Analyzers". In Neural Computation (1999).
[67]
[68]
Christopher M. Bishop. "Pattern Recognition and Machine Learning", Information Science
and Statistics. Springer, 2006.
[69] Lénaic Chizat et al. "Faster Wasserstein Distance Estimation with the Sinkhorn Divergence".
In Advances in Neural Information Processing Systems. 2020.
[70] Rémi Flamary et al. "POT: Python Optimal Transport". In Journal of Machine Learning Research (2021).
[71] Friedemann Zenke and Tim Vogels. "The Remarkable Robustness of Surrogate Gradient Learning for Instilling Complex Function in Spiking Neural Networks". In Neural Computa- tion (2021).
[72] Magnus R. Hestenes and Eduard Stiefel. "Methods of Conjugate Gradients for Solving Linear Systems". In Journal of Research of the National Bureau of Standards (1952).
[73] James Bradbury et al. "JAX: Composable transformations of Python + NumPy programs".
2018.
[74] Adam Paszke et al. "PyTorch: An Imperative Style, High-Performance Deep Learning Library". In Advances in Neural Information Processing Systems. 2019.
[75] Youcef Saad and Martin Schultz. "GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems". In Journal on Scientific and Statistical Computing (1986).
[76] H. A. Van der Vorst. "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems". In Journal on Scientific and Statistical Computing (1992).
[77] Luhuan Wu et al. "Practical and Asymptotically Exact Conditional Sampling in Diffussonn Models", In Advances in Neural Information Processing Systems. 2023.
13

---

[78] Nicolas Bonneel et al. "Sliced and Radon Wasserstein Barycenters of Measures". In Journal of Mathematical Imaging and Vision (2015).
[79]
Kimia Nadjahi et al. "Statistical and Topological Properties of Sliced Probability Diver- gences". In Advances in Neural Information Processing Systems. 2020.
[80]
Giannis Daras al. "Ambient Diffusion: Learning Clean Distributions from Corrupted Data". In Advances in Neural Information Processing Systems. 2023.
[81] Alex Krizhevsky, Geoffrey Hinton, al. "Learning Multiple Layers of Features from Tiny
Images". 2009.
[82]
Olaf Ronneberger et al. "U-Net:
Convolutional Networks for Biomedical Image Segmenta- tion". In Medical Image Computing and Computer-ssisted Intervention. 2015.
[83] Tim Salimans et al. "Improved Techniques for Training GANs". In Advances in Neural Information Processing Systems. 2016.
[84]
Martin Heusel et al. "GANs Trained by Two Time-Scale Update Rule Converge to Local Nash Equilibrium". In Advances in Neural Information Processing Systems. 2017. Bahjat Kawar et al. "GSURE-Basce Diffusion Model Training with Corrupted Data". In Transactions on Machine Learning Research (2024).
[85]
[86]
Ajil Jalal et al. "Robust Compressed Sensing MRI with Deep Generative Priors". In Advances in Neural Information Processing Systems. 2021.
[87] Yochai Blau and Tomer Michaeli. "The Verception-Distortin Tradeoff". In Conference on Computer Vision and Pattern Recognition. 2018.
[88] Mauricio Delbracio and Peyman Milanfar. "Inversion by Direct Iteration: An Alternative to Denoising Diffusion for Image Restoration". In Transactions Machine Learning Research (2023).
[89] E. G. Tabak and Cristina V. Turner. "A family of nonparametric density estimation algo- rithms". In Communications on Pure and Applied Mathematics (2013).
[90] Danilo Rezende and Shakir Mohamed. "Variational Inference with Normalizing Flows". In Proceedings of the 32nd International Conference on Machine Learning. 2015.
[91] Asad Aali et al. "Ambient Diffusion Posterior Sampling: Solving Inverse Problems with Diffusion Models trained Corrupted Data". 2024.
[92] Giannis Daras al. "Consistent Diffusion Meets Tweedie: Training Exact Ambient Diffusion Models with Noisy Data". 2024.
[93] Xinyu Peng et al. "Improving Diffusion Models for Inverse Problems Using Optimal Posterior Covariance". 2024.
[94] Jiaming Song et al. "Loss-Guided Diffusion Models for Plug-and- Controllable Genera- tion". In Proceedings of the 40th International Conference on Machine Learning. 2023.
[95] Arpit Bansal et al. "Universal Guidance for Diffusion Models". In International Conference Learning Representations. 2024.
[96] Yutong He al. "Manifold Preserving Guided Diffusion". In International Conference on Learning Representations. 2024.
[97] Prafulla Dhariwal and Alexander Quinn Nichol. "Diffusion Models Beat GANs on Image
Synthesis". In Advances in Neural Information Processing Systems. 2021.
[98]
Jonathan Ho and Tim Salimans. "Classifier Diffusion Guidance"". 2022.
[99]
Tero Karras al. "A Style-Basee Based Generator Architecture for Generative Adversarial Net- works". In Conference on Computer Vision and Pattern Recognition. 2019.
[100] Ashish Vaswani et al. "Attention is All you Need". In Advances in Neural Information
Processing Systems. 2017.
[101] Stefan Elfwing et al. "Sigmoid-weightee linear units for neural network function approxima- tion in reinforcement learning".. In Neural Networks (2018).
[102] Jimmy Lei Ba et al. "Layer Normalization". 2016.
[103] Diederik P. Kingma and Jimmy Ba. "Adam: A Method for Stochastic Optimization", In International Conference Learning Representations. 2015.
[104] Kaiming He et al. "Deep Residual Learning for Image Recognition". Conference on
Computer Vision and Pattern Recognition. 2016.
14

---

[105]
William Peebles and Saining Xie. "Scalable Diffusion Models with Transformers" In Inter- national Conference on Computer Vision. 2023.
[106]
Anton Obukhov et al. "High-fidelity performance metrics for generative models in PyTorch". 2020.
[107]
Wenzhe Shi et "Real-Time Single Image and Video Super-Resolution Using an Efficient Sub-Pixel Convolutional Neural Network". In Conference on Computer Vision and Pattern Recognition. 2016.
[108] Richard Zhang et al. "The Unreasonable Effectiveness of Deep Features as a Perceptual Metric". In Conference
Computer Vision and Puttern Recognition. 2018.
[109] Zhou Wang et al. "Image quality assessment: from error visibility to structural similarity". In Transactions Image Processing (2004).
15

---

A
Algorithms
Algorithm Expectation-maxxiizztion pipeline
Choose an initial prior (o(x))
Initialize the parameters the denoiser dæ t)
for k to K do
for = to S do
1/1 V/, A: p(y, A)
24 ~ Qx-1 Vis As)
Posterior sampling
repeat
7~ U(!,,....,..
t~1(011)
10
~~,(0,11)
11
21 I; 002
12
ft delldo(z1, (Z),, 2: 1|
# Denoising score matching
13
GRADIENTDESCENTT(R. Vol)
14
until convergence
15
Ok
16
Tk = 40
17
return 0K
Algorithm DDIM-style posterior sampling
T1 E0, for
= to do
f- i-1/1
- i/1
f + 34 Et sa(
y, A)
# Estimate B[x T4, y: A)
2 N (0, T)
0²
Ft
o² σi
I. - i - os
1 - 17
+ 0,
1
2
0
σt
return T0
Algorithm 3 Moment matching posterior score
function so(xe A) & + dg( t) Tweedie then
# Estimate Vz, log qa
y. A)
Enjor + Eg Vz, N2 (2)
else
Eeje, + Eg or (1 + E 1-1) (Ezl + +
+ EEy + AE,1 A A)
-1
(y Aa)
# Solve with conjugate gradient method
Sy f- Vz, ( t) Tuu
Estimate Vz, log q I A) Estimate Vg, log (xx))
82 + Σ, (i xz)
10
return 8 Syl
16

---

Algorithm 4 Conjugate gradient method
function CONUGATTEGRADEENTAA.... zo)
To + 6- Azo
Po < To for 0 N
if |r:ill < c then
return I;
r. T2 0g Pi Ap:
Ri+1 < I; + asps
T141 &-;; r; Api
1i+1
10
B,
T Ti
11
Pi41 < + Bipt
12
return EN
17

---

B Tweedie's formulae
Theorem 1. For any distribution p(x) and p(æ, | Mar, 2, E.), the first and second moments of the distribution p(z xx4) are linked to the score function Vx, log p(x)) through
2 Ø4 = 24 Σ Var logp(x;)
(25)
V/ xt. = Eg E log p(x)) Ex
(26)
We provide proofs of Theorem for completeness, even though it is well known result [55-58].
Proof.
1
Vz, p(z)
Van p(2t)
p(zs)
1
/ Va, p(x, dx
plaa))
/ p((,,0)) Vx, log plr, xr) dx
p(xx))
f p(z (x)) Van log logp(ar x) dx
/ p(z (x)) ا (xr - xz) dx
Ex Efx
[r - Σx x²
Proof.
22 logp(x;)
v, V. log p(xx)
nr Efx xil 1 -x11
( Ve, p( x1)) x do Dx - Ej²
(ʃ p(ar xe) Va, log
pau, x)) T dæ 5 - x
p(x))
(ʃ p(x xx) Er (x Elv xzl) x dx
Do1 - Eو¹
2T B [xT
[x - E|r
mrJ Bi1 xr]
V/2 srd] Ex1 x
18

---

C Experiment details
All experiments are implemented within the JAX [73] automatic differentiation framework. The code for all experiments is made available at https: / //gi tithhb..a.. comm/francoi--ooee// diffusion-pree prriors..
Diffusion models As mentioned in Section 2, in this work, we adopt the variance exploding SDE [45] and the denoiser parameterization [51]. Following Karras al. [51], precondition denoiser de (22, t) as
Gt Vo +1
St
do Ft, 1)
ho
log ot
(27)
0, +
Vo? +
where h0 log o)) is an arbitrary noise-conditional network. The scalar log is embedded vector using sinusoidal positional encoding [100]. our experiments, we use an exponential noise schedule
o2 = exp ((1 - t) log 10-3 + tlog 10") loss weights d = of + 1 and sample from Beta distribution Biaa = 3, B 3) during training.
(28)
Low-dimensional manifold The conditionl network (z, o a multi-layer perceptron with 3 hidden layers of 256 neurons and SiLU [101] activation functions. A layer normalization [102] function is inserted after each activation. The input of the network the concatenation of zt and the noise embedding vector. We train the network with Algorithm I for K 32 EM iterations. Each iteration consists 16 384 optimization steps of the Adam [103] optimizer. The optimizer and learning rate are re-initialized after each EM iteration. Other hyperparameters are provided in Table
Table Hyperparameters for the low-dimensional manifold experiment.
Architecture
MLP
Input shape
(5) (256, 256, 256)
Hidden features
Activation
SiLU
Normalization
LayerNorm
Optimizer
Adam 0.0
Weight decay
Scheduler
linear
Initial learning rate
1 x 10-3
Final learning rate
1 x 10-6
Gradient norm clipping
1.0 1024
Batch size
Steps per EM iteration
16384
EM iterations
We apply Algorithm 3 to estimate the posterior score Va, log p(x;) and truncate Algorithm to iterations. We rely the predictor-correcto [23, 46] sampling scheme sample from the posterior ( y, A). Following Rozet et al. [23], the predictor is a deterministic DDIM [50] step and the corrector is a Langevin Monte Carlo step. We perform 4096 prediction steps, each followed by correction step. At each EM iteration, we generate a single latent J for each pair (v. A).
We generate smooth random manifolds according to a procedure described by Zenke et al. [71]. We evaluate the Sinkhorn divergences using the POT [70] package with an entropic regularization factor A = le 3.
Corrupted CIFAR-10 The noise-conditional network (æ log o)) is a U-Net [82] with residual blocks [104], SiLU [101] activation functions and layer normalization [102]. Each residual block modulated with respect to the noise a; the style of diffusion transformers [105]. multi-head self-attention block [100] is inserted after each residual block at the last level of the U-Net. We train
19

---

the network with Algorithm 1 for K = 32 EM iterations, Each iteration consists of 256 epochs over the training set (50 000 images). To prevent overfitting, images are augmented with horizontal flips and hue shifts. The optimizer is re-initialized after each EM iteration. Other hyperparameters are provided in Table 3.
Table 3. Hyperparameters for the corrupted CIFAR-10 and accelerated MRI experiments.
Experiment
corrupted CIFAR-10
accelerated MRI
Architecture
U-Net (32, 32, 3) (5, 5, 5) (128, 256, 384)
U-Net
Input shape
(80, 80, 16)
Residual blocks per level
(3, 3. 3) (128, 256, 384, 512)
Channels per level
Attention heads per level
(0, 4, 0) 3 SiLU LayerNorm
(0, O. 0.
Kernel size
3 SiLU
Activation
Normalization
LayerNorm
Optimizer
Adam 0.0 2 x 10-4
Adam
Weight decay
0.0
Learning rate
10-4
Gradient norm cliping
1.0
1.0
EMA decay
0,9999
0.999
Dropout
0.1
0.1
Augmentation
h-flip, hue 256
h-flip, pad crop
Batch size
256
Epochs per EM iteration EM iterations
256
64
32
16
We apply Algorithm 2 with T = 256 discretization steps and 17 = 1 to sample from the posterior (x A). We apply Algorithm 3 with several heuristics for V[x | xx]] to compare their results against Tweedie's 's covariance formula. For the latter, we truncate the conjugate gradient method in Algorithm 4 single iteration. At each EM iteration, generate single latent x for each pair (y, A). Each FM iteration (including sampling and training) takes around 4 on A (40GB) GPUs
We evaluate the Inception score (IS) [83] and Fréchet Inception distance (FID) [84] of generated images using the torch-fidelity [106] package.
Accelerated MRI The noise-conditional network architecture the same as for the corrupted
CIFAR-10 10 experiment. The 320 320 x 1 tensor Ia is reshaped into 80 x 80 x 16 tensor using pixel shuffling [107] before entering the network. We train the network with Algorithm for = 16 EM iterations. Each iteration consists of 64 epochs over the training set (2 x 24 853 images). To prevent overfitting, images are augmented with horizontal flips and random crops. The optimizer is re-initialized after each EM iteration. Other hyperparameters
provided in Table 3.
We apply Algorithm with T 64 discretization steps and 17 to sample from the posterior go (x: y, A). We truncate the conjugate gradient method in Algorithm 4 to 3 iterations. At each EM iteration, we generate latents for each pair (y, A), which acts as data augmentation. Each EM iteration (including sampling and training) takes around 3 h on 4 A100 (40GB) GPUs.
20

---

D Additional figures
2
-2-
-2.5
0.0
2.5
-25
0.0
2.5
-25
0.0
2.5
-25
0.0
2.5
225
0.0
255
.70
22
33
Figure 7. 1-d and 2-d marginals of the ground-truth distribution p(x) used the low-dimensional manifold experiment. The distribution lies on random 1-dimensional manifold embedded R.
21

---

Figure 8. Example of samples from the posterior qe, (x y) along the EM iterations for the CIFAR-10
experiment. The generated images become gradually more detailed and less noisy.
Figure Example of samples from the posterior qo, (x y) along the EM iterations for the CIFAR-10 experiment when the heuristic (I + ) used for V(x x2).. The generated images become gradually more detailed but some noise remains.
Figure 10. Example of samples from the posterior qo, (ar y) along the EM iterations for the CIFAR-
10 experiment when the heuristic Σ, is used for V z 22). The generated images remain very noisy.
22

---

Figure 11, Example of scan slices from the fastMRI [7, 8] dataset.
Figure 12. Example of k-space sub-sampling observations with acceleration factor R - 6 for the accelerated MRI experiment. We represent each observation by its zero-filled inverse, where missing frequencies are set to zero before taking the inverse discrete Fourier transform.
23

---

Figure 13. Example of samples from the final model qe, (z) for the accelerated MRI experiment. The samples present varied and coherent global structures. Samples seem slightly less sharp than real scans (see Figure 11), but do not present artifacts typical to unresolved frequencies (see Figure 12).
Ground--rrnt trrnh
Obeevvation
Sumple
Sample
Figure 14. Examples of posterior samples using diffusion prior trained from k-space observations only. The forward process crops the latent r to centered 160 x 160 window. Moment matching posterior sampling is used sample from the posterior. Samples are consistent with the ground -ouund-truth where observed, but present plausible variations elsewhere.
24

---

Figure 15. Example of samples from the model qo, (x) after k 2 EM iterations for the accelerated
MRI experiment when the heuristic (I + E) is used for V(x 2z). The samples start to present vertical artifacts due to poor sampling.
Figure 16. Example of samples from the model qo, (x) after k 4 EM iterations for the accelerated MRI experiment when the heuristic (I + Σ, ">- is used for V[x 2z). The artifacts introduced by the poor sampling get amplified at each iteration, leading to total collapse after few iterations.
25

---

E Evaluation of MMPS
In this section, we evaluate the moment matching posterior sampling (MMPS) method presented in Section 4.2 independently from the context of learning from observations. The code for this section made available at https: https:///gtttin www/francois---oeeeeooooooo. wwwpa//oooossottttooooeeeossssse bencharrkbb.cca..........
Tasks We consider four linear inverse problems on the 256 x 256 FFHQ [99] dataset. (i) For box inpainting, we mask out randomly positioned 128 x 128 square of pixels and add a large amount of noise (Øy = 1). (ii) For random inpainting, we randomly delete pixels with 98 % probability and add a small amount of noise (0g = 0.01). (iii) For motion deblur, we apply a randomly generated 61 x 61 motion blur kernel and add medium amount of noise (a, y = 0.1). (iv) For super resolution, we apply a 4x bicubic downsampling and add a medium amount of noise (oy 0.1).
Methods For all inverse problems, we use the pre-trained diffusion model provided by Chung et al. [21] as diffusion prior. We adapt and extend the DPS [21] codebase to support MMPS as well as DitIPIR [26], LIGDM [22] and TMPD [25]. We use the DDIM [50] sampler with rj 1 for all methods, which is equivalent to the DDPM [16] sampler. We fine-tune the hyperparameters of DPS (C' = 0.5) and DifIPRR (A 8.0) to have the best results across the four tasks. With MMPS, we find that the Jacobian of the pre-trained model provided by Chung et al. [21] is strongly non-symmetric and non-definite for large oee which leads to unstable conjugate gradient (CG) [72] iterations. We therefore replace the CG solver with the GMRES [75] solver, which can solve non-symmetric non-definite linear systems.
Protocol We generate one observation per inverse problem for 100 images' of the FFHQ [99] dataset. We generate sample for each observation with all considered posterior sampling methods. All methods are executed with the same random seed. We compute three standard image reconstruc- tion metrics- LPIPS [108], PSNR and SSIM [109] for each sample and report their average in Table 4. We present generated samples for each inverse problem in Figures 17, 18 and 19.
As side note, emphasize that reconstruction metrics do not necessarily reflect the accuracy of the inferred posterior distribution, which we eventually care about. For example, PSNR and SSIM [109] favor smooth predictions such the mean E[z y] over actual samples from the posterior p(x y). Conversely, LPIPS [108] favors predictions which are perceptually similar to the reference, even if they are distorted. In general, it is impossible to simultancously optimize for all reconstruction metrics [87, 88].
Results MMPS consistently outperforms all baselines, both qualitatively and quantitatively. As expected, performing more solver iterations improves the sample quality, especially when the Gram matrix AA is strongly non-diagonal, which is the case for the motion deblur task. However, the improvement shows rapidly diminishing returns, as the difference between and 3 iterations much larger than between 3 and 5. MMPS is also remarkably stable with respect the number of sampling steps in contrast to DPS [21]. DifiPIR [26] and IIGDM [22] which are sensitive to the number of steps and choice of hyperparameters. Finally, MMPS requires fewer sampling steps reach the same image quality as previous methods, which largely makes up for its slightly higher step cost.
'Chung et al. [21] do not indicate which subset FFHQ [99] was used to train their model. Without further information, we choose to use the first 100 images for evaluation, which could lead to biased metrics if the diffusion prior was trained on them. However, since we use the same diffusion prior for all posterior sampling methods, the evaluation remains fair.
26

---

Table 4. Quantitative evaluation of MMPS with 1,3 and 5 solver iterations.
Box inpainting
Random inpainting
Motion deblur
Super resolation
Mothod
Steps
LPIPS
PSNR
SSIM
LPIPS
PSNR
SSIM
LPIPS
PSNR
SSIM
LPIPS
PSNR 1 SSIM 1
DIIPIR [26]
10
0.33
19.17
0,50
0.78
10.97
0.32
0.24
24.54
0.72
0.20
26.63
0.78
DiAPIR [26]
100
0.30
18.15
0.54
0.68
10.26
0.25
0.19
23.97
0.70
0.17
25.24
0.73
DifiPR [26]
1000
0.33
17.39
0.49
0.74
9.51
0.21
0.17
23.55
0.67
0.15
24.72
0.70
[21]
[21]
0.38
16.82
0.50
0.39
16.66
0.49
0.29
19.75
0.57
0.35
18.29
0.54
DPS [21]
1000
0.22
21.01
0.64
0.19
21.90
0.66
0.18
22.91
0.66
0.16
25.02
0.72
IIDM [22]
10
0.40
18.94
0.61
0.59
11.28
0.40
0.25
25.83
0.76
0.25
26.42
0.77
HIGD [22]
100
0.44
18.23
0.47
0.39
17.03
0.48
0.25
22.37
0.61
0.15
25.63
0.71
IGDM [22]
1000
0.81
14,80
0.31
0.14
22.32
0.69
1.06
13.12
0.21
0.64
18.41
0.29
TMPD [25]
10
0.36
19.90
0.64
0.59
11.08
0.40
0.27
25.28
0.74
0.26
26.07
0.76
TMPD [25]
0.27
19.86
0.64
0.58
10.73
0.31
0.17
26.22
0.76
0.17
26.79
0.77
[25]
25.91
0.76
MMPS (I)
0.27
21.19
0.68
0.26
22.41
0.69
0.33
22.12 12
0.66
0.24
26.94
0.78
MMPS (1)
100
0.20
21.19
0.67
0.18
22.18
0.69
0.20
23.92
0.71
0.15
27.32
0.79
MMPS (I)
1000
0.19
20.77
0.64
0.18
21.94
0.66
0.16
23.83
0.69
0.12
26.92
0.77
MMPS (3)
10
0.26
21.55
0.68
0.21
23.58
0.74
0.24
25.33
0.75
0.19
27.94
0.81
MMPS
0,20
21.29
0.67
0.15
22.76
0.71
0.15
26.16
0.76
0.13
27.18
0.78
21.01
0.15
22.45
0.12
(5)
0.23
21.73
0.69
0.20
23.72
0.75
0.20
26.70
0.78
0.18
28.02
0.81
MMPS (5)
100
0.20
21.30
0.67
0.15
22.82
0.72
0.13
26.70
0.77
0.13
27.12
0.78
MMPS (5)
1000
0.20
20.98
0.64
0.14
22.52
0.69
0.11
26.18
0.75
0.11
26.60
0.76
Table 5. Time and memory complexity of MMPS for the 4x super resolution task. Each solver iteration increases the time per step by around 16 ms. The maximum memory allocated by MMPS is about 10% % larger than DPS [21] and ITGDM [22].
Method
VJPs
Time [ms/stepp]
Memory [GB]
DiffPIR [26]
0
30.2
0.66
DPS [21]
1
40.5
2.29
ITGDM [22]
47.6
2,30
TMPD [25]
2
62.2
2.52
MMPS (1)
2
58.0
2.52
MMPS (3)
4
90.1
2.52
MMPS (5)
6
122.1
2.52
27

---

MMPS
MMPS
MMPS
MMMS
Reference
Oluervation
steps
steps
Tlo impareiin
Fand impareiing
Motion debdia
Suger remalaatio
Figure 17. Qualitative evaluation MMPS with and 5 solver iterations.
DPS 100 steque
THiPERR 10 xtrgps
DiAPIR 100 -trps
DEIPIR 1000 strp
Reference
Observation
10 strgw
1000-
Box ingaiiiin
ingainting
Metios debltx
Sugar moodeeiio
Figure 18, Qualitative evaluation of DPS [21] and DifiPIR [26].
28

---

HCCM 10
RIDDM 100 steps
HGDM 1000 steps
TMPD
TMPD 100 skeps
TMPD 1000 steps
Refenmec
10 stepe
Foxe inpaieiing
Paant inpareiing
Metion caldtar
Sigger remalaiion
Figure 19. Qualitative evaluation of IIGDM [22] and TMPD [25].
29

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We only claim to present new method, which we describe in Section 4. We compare our method against previous ones in Sections and 6.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitations in Section 7.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
30

---

Justification:
No new theoretical results are presented, but methods are motivated by estab-
lished literature.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The manuscript describes all methods and experiments. Algorithms are provided for the methods. The code for all experiments is made available at https: / ggithub.com// francois-rozet/diffue ppriors.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
31

---

The authors should provide scripts to reproduce all experimental results for the new reproducible, they Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes)
Justification: The code for all experiments made available at https://gii /ttpp://ggthub..... thub. com/ francois- stttaaarrrztt/ffdiiooooooo /ttdd/fiininiriioiocoo Instructions to acquire the datasets are provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (uttps://nips.. cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
/ ctttpp..cc/ubbiceoc htw/guidesss des CoodeSubmis ssi assonnolicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
proposed method and baselines. only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Experiment details provided in Appendix C. Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: Computing error bars for Table I and Figure 4 would require retraining every model several times for different datasets.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
32

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Experiment details are provided in Appendix C. We acknowledge the use of computer cluster in the acknowledgments section. Preliminary experiments are not reported in the manuscript.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: : /t/p/euur.pp./. http://ttrppc.ppppiiii////////// Cccc/uubbbi/Bthhttsssdeeeiiee
Answer: [Yes]
Justification:
We have reviewed and agree with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is no societal impact of the work performed.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
33

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We acknowledge the of the NYU fastMRI dataset [7, 8] in the acknow]- edgments section and follow its terms of use.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
34

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com /atasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
35

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
36


Reviews:
[{"review_id": "YRUBG6OaMG", "reviewer": "Reviewer_3yec", "summary": "The authors propose a new framework for training diffusion models from corrupted data, based on the Expectation-Maximization algorithm. Before this work, there were two approaches to the problem of learning from corrupted data: Ambient Diffusion and SURE. Ambient Diffusion was designed for linear inverse problems and SURE for denoising. The proposed method provides a different, unified methodology for training diffusion models from corrupted data. The authors show experimentally that their method outperforms the Ambient Diffusion baseline in the same corruption setting.", "questions": "See weaknesses above. I would be happy to further increase my score if my concerns are properly addressed.", "limitations": "The authors have adequately addressed the limitations of their method.", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1) The topic of learning diffusion models from corrupted data is very important and very interesting. The submission is timely and relevant as the interest in this research topic is growing.\n2) The authors propose a fresh idea in the space of learning diffusion models from corrupted data. The proposed idea has several nice properties: i) it provides a unified treatment for different corruptions, ii) it is simple to implement, and iii) it has strong experimental performance.\n3) The presentation of the paper is really good. The authors motivate the problem, present a new principled method to solve this, and show experimental results.\n4) A byproduct of this work is a new method for solving inverse problems with diffusion models. While still approximate, this new method seems to be more principled compared to previous approaches.\n\nOverall, this is a strong submission that offers new insights into the problem of learning from corrupted data.", "weaknesses": "The paper has also several weaknesses.\n\n1) The proposed framework introduces significant computational and engineering overhead. Even for local convergence, EM typically requires many steps. In this setting, each step is a new training of a diffusion model.\n2) The proposed algorithm is still approximate. The authors currently mention that their work is the first one that leads to \"proper\" denoisers. Yet, there is an approximation happening when the authors use diffusion models to perform posterior sampling. In fact, the authors had to develop a whole new method for solving inverse problems with diffusion models to get good results. The authors should emphasize this limitation more.\n3) There is a very large field of prior works in solving inverse problems with diffusion models. This paper proposes a new method to achieve that goal. This method should be independently tested and benchmarked, decoupled by the context of learning with corrupted data. If the method is truly superior to prior works, this would be a very interesting byproduct of this paper and it's worth understanding this. If not, there is still value in the proposed method for the context of learning from corrupted data, but there is more to understand regarding why it works so well in this particular setting.\n4) The authors currently provide comparisons with Ambient Diffusion only at the setting of $p=0.75$. It would be nice to see how their method scales for different corruption levels. It would also be helpful to provide comparisons in other datasets apart from CIFAR-10.\n5) Since this method can be in principle used to train diffusion models from various inverse problems, it would be nice to see results for learning from other corruption types, e.g. blurry data.\n6) Ambient Diffusion has been used to solve MRI inverse problems in the paper: \"Ambient Diffusion Posterior Sampling: Solving Inverse Problems with Diffusion Models trained on Corrupted Data\". It would be useful to provide comparisons with this work.\n7) The authors do not provide comparisons for the additive Gaussian noise case. A natural baseline would be the SURE method. I would also like to bring the paper: \"Consistent Diffusion Meets Tweedie\" to the attention of the authors.\n8) For some corruption types, it is impossible to ever reconstruct the true distribution from noisy observations. The authors do not discuss how the proposed algorithm would work in such cases.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your in-depth review and the legitimate concerns you have raised.\n\n* **W1** Indeed, this is one of the limitation of our method, which we mention in Section 3. However, we would like to mention that in our pipeline we start each training step with the previous parameters, which reduces the cost of each training step. Overall, the entire EM procedure for the corrupted CIFAR-10 experiment takes around 4 days for 32 EM iterations on 4 A100 GPUs (see Appendix C), which is similar to the training of AmbientDiffusion [77].\n\n* **W2** Our EM method is indeed sensitive to the quality of posterior samples, as mentioned in Sections 4.1, 4.2 and 7. We propose to add the following sentence at line 291.\n\n  > [...] sensitive to the quality of posterior samples. In fact, we find that previous posterior sampling methods [21, 22] lead to disappointing results, which motivates us to develop a better one.\n\n* **W3** We appreciate that you recognize MMPS as a valuable contribution in its own right. We follow your suggestion and benchmark MMPS against previous methods (DPS and $\\Pi$GDM). We invite you to consult the global rebuttal regarding these additional experiments.\n\n* **W4** We follow your suggestion and repeat the corrupted CIFAR-10 experiment at more corruption levels (0.25, 0.50 and 0.75). We invite you to consult the global rebuttal regarding these additional experiments.\n\n  We note that the other experiments of AmbientDiffusion [77] also use a diagonal measurement $A$ (a mask), which is why we conduct the accelerted MRI experiment, where the measurement is a more challenging undersampled Fourier transform.\n\n* **W5** Indeed, a strengh of our method is that it can handle a wide variety of measurement types, or even learning from several sources with different measurement types. Our three experiments present different measurement types: a linear projection $A \\in \\mathbb{R}^{5 \\times 2}$, a random binary masking, and an undersampled Fourier transform. We note that Gaussian blur, as you suggest, is equivalent to masking the high frequencies of an image, which is similar to the undersampled Fourier transform in the accelerated MRI experiment.\n\n* Concerning **W6** and **W7**, we thank you for bringing these concurrent works to our attention. We will discuss them in our related work section. Are we correct in understanding that in \"Consistent Diffusion Meets Tweedie\" the goal is to fine-tune a pre-trained diffusion model using data corrupted by isotropic Gaussian noise?\n\n* **W8** Indeed, we forgot to explain how our method would behave when $p(x)$ cannot be uniquely identified from the observations. We propose to replace lines 298-300 with the following paragraph\n\n  > Finally, as mentioned in Section 6, empirical Bayes is an ill-posed problem in that distinct prior distributions may result in the same distribution over observations. In other words, it is generally impossible to identify \"the\" ground-truth distribution $p(x)$ given an empirical distribution of observations $p(y)$. Instead, for a sufficiently expressive diffusion model, our EM method will eventually converge to a prior $q_\\theta(x)$ that is consistent with $p(y)$, but generally different from $p(x)$. In future work, we would like to follow the maximum entropy principle, as advocated by Vetter et al. [36], so as not to reject any possible hypothesis.\n\n  We emphasize that this identifiability issue is a limitation of the problem itself, and not of our method.\n\nWe believe that this rebuttal addresses most of your concerns and, therefore, kindly ask you to reconsider your score.", "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": "zicb5QdzPs", "reviewer": "Reviewer_XeV5", "summary": "This paper proposes a novel solution to a specific class of Empirical Bayes problems. The method is especially useful when the latent variable and the observations are closely related, such as in cases where the observations are incomplete or noisy latent variables. The major technical obstacle is modeling the posterior distribution for sampling and training, which is addressed through several approximations under Gaussianity.", "questions": "1. To what extent does Tweedie's formula violate the assumption that $\\mathbb{V}(x|x_t)$ is independent of $x_t$? Are there any bounds or asymptotic results?\n2. What is the role of $\\Sigma_y$? To what extent does it influence overall performance? What if $\\Sigma_y$ is anisotropic or non-diagonal? What if the Gaussian noise is not as small as it is in the experiments?\n3. How does the linearity assumption restrict performance? How are the priors over $A$ specified? A real-life example with explicit expressions of $A$ and $\\Sigma_y$ would be helpful.\n4. Are there any experiments on the improvement per iteration of the Conjugate Gradient method?", "limitations": "Yes, the authors adequately addressed the limitations.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. **Originality-Middle:** A novel combination of the Bayesian Inverse Problem and Diffusion Model under Gaussianity assumptions.\n2. **Quality-Middle:** Well-organized overall, but lacking sufficient details.\n3. **Clarity-Middle:** Clear intuitively and qualitatively, but requires more quantitative analysis.\n4. **Significance-Middle:** Impressive idea and experimental results, though applications are restricted.", "weaknesses": "1. There are not sufficient benchmark experiments, especially those involving other Empirical Bayes methods. A toy example with 3-4 benchmarks, focusing specifically on the trade-off between time consumption and model performance, would help justify the soundness of the model.\n2. The validity of the approximation needs to be shown theoretically. See the \"Questions\" section below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the pertinent questions you have asked.\n\n* **W1** Although we agree that a benchmark with previous empirical Bayes methods would be valuable for the statistical inference community, we don't think it would be relevant to justify our work. First, our method is based on the established EM algorithm, which has stood the test of time. Second, our goal is to train diffusion models from observations as they are best-in-class for modeling high-dimensional distributions and proved to be remarkable priors for Bayesian inference. However, as explained in Section 1 and 6, previous empirical Bayes methods are not applicable to diffusion models. These methods are also typically bound to low-dimensional latent spaces. For these reasons, it is challenging to design a benchmark between our work and these previous EB methods that would be both fair and informative. Instead, we choose to benchmark our work against similar methods in the diffusion model literature.\n\n* **Q1** This assumption is actually the same as the Gaussian approximation of Eq. (17). Indeed, assuming that $p(x \\mid x_t)$ is Gaussian, we have (Bishop [67])\n\n  $$ \\mathbb{E}[x \\mid x_t, y] = \\mathbb{E}[x \\mid x_t] + \\mathbb{V}[x \\mid x_t] A^\\top (\\Sigma_y + A \\mathbb{V}[x \\mid x_t] A^\\top)^{-1} (y - A \\mathbb{E}[x \\mid x_t]) $$\n\n  but Tweedie's formulae also gives\n\n  $$ \\mathbb{E}[x \\mid x_t, y] = x_t + \\Sigma_t \\nabla_{x_t} \\log p(x_t \\mid y) = \\mathbb{E}[x \\mid x_t] + \\Sigma_t \\nabla_{x_t} \\log p(y \\mid x_t) $$\n\n  Therefore, we have\n\n  $$ \\Sigma_t \\nabla_{x_t} \\log p(y \\mid x_t) = \\mathbb{V}[x \\mid x_t] A^\\top (\\Sigma_y + A \\mathbb{V}[x \\mid x_t] A^\\top)^{-1} (y - A \\mathbb{E}[x \\mid x_t]) $$\n\n  which is equivalent to Eq. (20) since $\\mathbb{V}[x \\mid x_t] = \\Sigma_t \\nabla_{x_t} \\mathbb{E}[x \\mid x_t]^\\top$.\n\n  Finzi et al. [24] provide a detailed analysis of the Gaussian approximation of Eq. (17) and prove (Theorem 2) that moments of order $n > 2$ converge to 0 at a rate of $\\sigma_t^n$.\n\n* **Q2** Our method does not make any assumptions on the covariance $\\Sigma_y$. It can be anisotropic and/or non-diagonal. In fact, it is always possible to rewrite the likelihood $\\mathcal{N}(y \\mid Ax, \\Sigma_y)$ as $\\mathcal{N}(\\Sigma_y^{-1/2} y \\mid \\Sigma_y^{-1/2} A x, I)$. However, you are right that the signal-to-noise ratio has an impact on our method. We expect larger noise levels to slow down the convergence of the EM algorithm, but lead to an equivalent stationary distribution, under the assumption of infinite data. The rationale is that it is possible to identify all the moments of $p(Ax)$ given $p(y)$ regardless of $\\Sigma_y$.\n\n* **Q3** The assumption of a linear Gaussian forward process $p(y \\mid x)$ does not \"restrict the performance\" of our method but limits its applicability. However, many real-life and scientific problems can be formalized with linear Gaussian forward processes. The accelerated MRI experiment is a good example. The measurement matrix $A$ and covariance $\\Sigma_y$ are defined by the specific instruments that are used to gather MRI scans. If the configuration or environment of the instruments changes, the measurement matrix $A$ and covariance $\\Sigma_y$ may also change. The prior $p(A)$ is therefore the empirical distribution of $A$ for the observations $y$ we have access to. To paraphrase, we do not specify $p(y, A)$, it is imposed upon us by the task at hand.\n\n  We propose to clarify this matter in the manuscript with the following sentence at line 82\n\n  > For example, if the position or environment of a sensor changes, the measurement matrix $A$ may also change, which leads to an empirical distribution of pairs $(y, A) \\sim p(y, A)$.\n\n* **Q4** There are no experiments on the improvement per iteration of the CG method currently, but we have conducted additional experiments to benchmark MMPS for this rebuttal and find that increasing the number of CG iterations improves image quality/sharpness, but only marginally and with rapidly diminishing returns. We invite you to consult the global rebuttal regarding these additional experiments.\n\nWe believe that this rebuttal addresses most of your concerns and, therefore, kindly ask you to reconsider your score.", "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": "nYnG0wod38", "reviewer": "Reviewer_s8XN", "summary": "This paper proposes a method to learn generative models from noisy and incomplete data. As opposed to general recent methodologies which require a clean unconditional dataset to solve inverse problems, the paper aims at providing a method to learn the prior from noisy data.", "questions": "The paper is on an interesting topic as I mentioned however many things are unclear and the algorithm/pipeline needs to be written in a much more clear way before this can be published.\n\n1) The algorithm does not appear at all in the main body of the paper. Please provide a clear and step-by-step demonstration of the training method in the main body of the paper.\n\n2) The authors put a prior on the matrix $A$. What is this prior?\n\n3) It seems that the main difference between the EM and this method is the parameterization of the prior (which is via a score network) described briefly in 4.1. The authors then briefly talk about how to estimate this parameter but then somehow the parameter disappears in the later sections. For example, in eq. (15), the authors say it is easy to estimate the prior score $\\nabla_x \\log p_t(x)$ but this would require unconditional (and clean!) samples. It is not clear how it is suggested that this quantity is easy to obtain?\n\n4) Section 4.2 describes the standard inverse problem solvers (or a very similar idea). However, as I said in the point above, this requires a pre-trained, clean sampler. This way of writing it is very unclear for readers.\n\n5) Algorithm 1 is the intended \"full pipeline\" yet the steps are not described clearly. I suggest authors to both move this into the main text as well as extend the discussion of the steps of this pipeline. Some questions are below.\n\n6) In Algorithm 1, is the posterior sampling done by a pretrained diffusion model? If so, this invalidates the main point of the paper. If not, then it means that $q_k$ using the score vector $d_{\\theta_k}$, is that a correct conclusion?\n\n7) In Algorithm 1, after sampling $x_i$ for $i = 1,\\ldots,S$, does one train the score network with this? In what sense is this an EM algorithm if score matching is done rather than maximum likelihood?\n\nI strongly suggest authors to clarify the pipeline clearly and precisely.", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper is on an interesting and timely topic; this is a much needed class of methods.", "weaknesses": "The paper is poorly written, the method is not clear, derivations are incomplete. See my detailed review in the questions part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely apologize for the inconvenience caused by the writing of our manuscript. We believe that most of your questions stem from the lack of clarity of Section 4.2, where we fail to state that MMPS is not bound to the EM context and can be applied to any diffusion prior. We propose to clarify this matter by replacing lines 108-112 in Section 4.2 with the following\n\n> To sample from the posterior distribution $q_\\theta(x \\mid y) \\propto q_\\theta(x) \\, p(y \\mid x)$ of our diffusion prior $q_\\theta(x)$, we have to estimate the posterior score $\\nabla_{x_t} \\log q_\\theta(x_t \\mid y)$ and plug it into the reverse SDE. In this section, we propose and motivate a new approximation for the posterior score. As this contribution is not bound to the context of EM, we temporarily switch back to the notations of Section 2 where our prior is denoted $p(x)$ instead of $q_\\theta(x)$.\n>\n> **Diffusion posterior sampling** Thanks to Bayes' rule, the posterior score $\\nabla_{x_t} \\log p(x_t \\mid y)$ can be decomposed into two terms [17, 18, 21-25, 53]\n> \n>  $$ \\nabla_{x_t} \\log p(x_t \\mid y) = \\nabla_{x_t} \\log p(x_t) + \\nabla_{x_t} \\log p(y \\mid x_t) \\, . $$\n>\n> As an estimate of the prior score $\\nabla_{x_t} \\log p(x_t)$ is already available via the denoiser $d_\\theta(x_t, t)$, the remaining task is to estimate the likelihood score $\\nabla_{x_t} \\log p(y \\mid x_t)$.\n\nWe now answer your questions, in light of this clarification.\n\n* **Q1** We describe the pipeline in plain text in Section 4.1 and as an algorithm in Appendix A.\n\n* **Q2** The measurement matrix $A$ is defined by the specific instruments that are used to gather the observations. If the configuration or environment of the instruments changes, the measurement matrix $A$ may also change. The prior $p(A)$ is therefore the empirical distribution of $A$ for the observations $y$ we have access to. To paraphrase, we do not specify $p(y, A)$, it is imposed upon us by the task at hand.\n\n  We propose to clarify this matter in the manuscript with the following sentence at line 82\n\n  > For example, if the position or environment of a sensor changes, the measurement matrix $A$ may also change, which leads to an empirical distribution of pairs $(y, A) \\sim p(y, A)$.\n\n* **Q3** You are right. The main difference is the parameterization of the prior with a diffusion model. At each step we use the current diffusion prior $q_\\theta(x)$ to generate samples from the posterior(s) $q_\\theta(x \\mid y)$ via the proposed MMPS method.\n\n* **Q4** Indeed, MMPS can be used to solve any linear inverse problem with a diffusion prior. However, this prior does not need to be the distribution of clean data. It can be any prior, including our intermediate priors $q_{\\theta_k}(x)$.\n\n* **Q5** We describe the pipeline in plain text in Section 4.1. We unfortunately do not have the space to move Algorithm 1 into the main text.\n\n* **Q6** We do not use a pre-trained diffusion model. As you correctly concluded, $q_k(x) := q_{\\theta_k}(x)$ is the current diffusion prior parameterized by the current denoiser $d_{\\theta_k}(x_t, t)$.\n\n* **Q7** We show in Section 4 that in the context of EB, the EM step is equivalent to minimizing the KL bewteen $\\pi_k(x)$ and $q_{\\theta_{k+1}}(x)$. For diffusion models, this KL minimization is generally conducted by denoising score matching or equivalent formulations (e.g. ELBO) using samples from $\\pi_k(x)$.\n\nWe believe that this rebuttal addresses most of your concerns and, therefore, kindly ask you to reconsider your score.", "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": "SGGwnTdHz3", "reviewer": "Reviewer_mKff", "summary": "This paper focuses on training diffusion models using incomplete or noisy data only, which is obtained through a linear measurement operator $A$. While prior works assume full rank of $E[A^TA]$ or $E[A^+A]$ and change the denoising score matching objective while training diffusion models, this paper uses a moment matching posterior sampling approach that does not modify the training objective. Yet, the authors demonstrate superior results both quantitatively and qualitatively on (i) a toy dataset (ii) corrupted CIFAR10 and (iii) accelerated MRI.", "questions": "Please refer to the weaknesses above.", "limitations": "Yes, the authors have adequately addressed the limitations.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. One major strength of the proposed method is that unlike prior works it does not modify the denoising score matching objective, which guarantees a proper diffusion model at every iteration. \n2. The experimental results in the toy setting clearly depicts the advantages of using higher order moments in pruning inconsistent regions.\n3. The paper is nicely written and the main contributions are clearly stated.", "weaknesses": "1. The core idea of moment matching posterior sampling (Section 4.2) has previously appeared in prior works STSL [1] and TMPD [25], where its benefits are clearly demonstrated in large-scale applications. \n2. The approximation used in Equation (21) computes gradient of the first order score which is already an approximation. This results in high time and memory complexity, which is precisely the reason why prior works [1,2,25, 65] seek alternatives. \n3. What is the typical range of $\\sigma_t$ in Fig. 2? It seems like the results are comparable in low noise regime. How does this observation translate into variance preserving SDE, which is the most commonly used form of SDEs for posterior sampling?\n4. Does Equation (22) hold for any vector $v$ as the equation suggests?\n5. Section 5: Taking gradients of the score becomes an issue especially in high dimensional setting. It is important to understand the time and memory complexity of the proposed algorithm in commonly used benchmarks such as FFHQ or ImageNet ( 256x256 or 512x512).\n6. Missing discussion and comparison with other baselines (e.g. [82]) highlighted in the related works Section 6.\n7. Theorem 1 in Appendix B is a well known result [63]. It'd be better to cite the original work and provide the key steps only for completeness. Or, the authors should highlight the key distinctions from the prior work. \n\nMissing Related Works\n\n[1] [Beyond First-Order Tweedie: Solving Inverse Problems using Latent Diffusion](https://openaccess.thecvf.com/content/CVPR2024/papers/Rout_Beyond_First-Order_Tweedie_Solving_Inverse_Problems_using_Latent_Diffusion_CVPR_2024_paper.pdf)\n\n[2] [Improving Diffusion Models for Inverse Problems\nUsing Optimal Posterior Covariance](https://arxiv.org/pdf/2402.02149)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the legitimate concerns you have raised.\n\n* **W1** Indeed, we are not the first to use the covariance $\\mathbb{V}[x \\mid x_t]$ to improve the approximation of the likelihood score. We believe that Finzi et al. [24] were the first to propose it, shortly followed by Boys et al. [25]. As explained in section 6, the difference of our approach resides in the use of the entire covariance matrix, without the need to materialize it (which is intractable) thanks to the conjugate gradient method. Instead, Boys et al. [25] use a row-sum approximation of the covariance's diagonal, which voids some of the benefits of using the covariance as illustrated in Figure 1 an 2. In addition, when the covariance is not diagonal, the row-sum approximation may result in null or negative values, which leads to total failures (NaNs) in our experiments.\n\n  We thank you for bringing STSL [1] to our attention. Are we correct in understanding that STSL uses the trace of the covariance/Hessian to guide the sampling, but does not justify it with a Gaussian approximation of $p(x \\mid x_t)$? If so, we are not sure how STSL is similar to the idea of TMPD [25] and/or MMPS.\n\n  In the manuscript we refer a few times to MMPS as a \"new\" posterior sampling scheme. We propose to replace these occurences with the word \"improved\". Would this be satisfactory?\n\n* **W2** Materializing the Jacobian $\\nabla_{x_t}^\\top d_\\theta(x_t, t)$ would indeed result in high time and memory complexity, as mentioned at lines 142-144. However, we never do. Instead, we leverage the CG method to solve the linear system in Eq. (22), which only requires a cheap vector-Jacobian product per iteration.\n\n* **W3** In our experiments, $\\sigma_t$ ranges between $10^{-3}$ and $10^2$. This figure can be easily translated to the VP SDE with the relation $\\bar{\\alpha}_t = \\frac{1}{1 + \\sigma_t^2}$. The divergence should also be scaled accordingly. We note that using the covariance $\\mathbb{V}[x | x_t]$ is orders of magnitude better than heuristics, especially in low noise regimes.\n\n* **W4** Yes, this equation simply rewrites $v = (\\Sigma_y + A \\mathbb{V}[x \\mid x_t] A^\\top)^{-1} (y - A \\mathbb{E}[x \\mid x_t])$ and has a solution as long as $\\Sigma_y + A \\mathbb{V}[x \\mid x_t] A^\\top$ is invertible, which is the case if $\\Sigma_y$ is SPD.\n\n* **W5** We follow your suggestion and benchmark MMPS against previous methods. We invite you to consult the global rebuttal regarding these additional experiments.\n\n* **W6** This is not true. Section 5.3 is dedicated to a comparison with GSURE-Diffusion [82] on the accelerated MRI experiment. We note that Kawar et al. [82] do not provide the identifiers of the scans they use for evaluation in their manuscript or code.\n\n* **W7** Indeed we provide these proofs only for completeness. We propose to add the following sentence in Appendix B for clarity\n\n  > We provide proofs of Theorem 1 for completeness, even though it is a well known result [62-65].\n\nWe believe that this rebuttal addresses most of your concerns and, therefore, kindly ask you to reconsider your score.", "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": "eBwLm7yk2L", "reviewer": "Reviewer_tnYw", "summary": "This paper tackles the problem of learning a diffusion model from incomplete and noisy observations. The problem is modelled as follows; the distribution of the noised and incomplete data is assumed to be $p(y) = \\int p(y|x) q(x) dx$. The authors introduce a parametric version of it $p^\\theta(y) = \\int p(y|x) q^\\theta(x) dx$ and then seek to learn $p^\\theta$ by minimizing the KL divergence $KL(p||p^\\theta)$. As $p^\\theta$ is a latent variable model, a natural way to learn it is via the Expectation-Maximization algorithm. The resulting algorithm is iterative and at each step, a diffusion model learned. Also, within each step, the E-step of EM requires sampling from the posterior diffusion $q^\\theta(x|y)$ which is intractable in practice. The authors resort to approximate inference and rely on recent advances in posterior sampling of Diffusion models. More specifically, to sample from $q^\\theta(x|y)$ it is required to be able to estimate the conditional score within the Diffusion. This conditional score involves computing the gradient log of the conditional distribution $p(y|x_t) = \\int p(y|x_0) q^\\theta _{0|t}(x_0 | x_t) dx_0$, where the time here referes to the Diffusion steps. The authors propose to replace $q^\\theta _{0|t}(x_0 | x_t)$ with its moment projection, which closely follows the work of Boys et al. 2023, and then integrating $p(y|x_0)$ against it. This step however requires an expensive matrix inversion, which the authors propose to approximate via conjugate gradient methods.", "questions": "N/A", "limitations": "see above", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper is well written and presents a natural, original and interesting methodology. The originality stems from the use of the EM algorithm which allows (1) learning at each step a Diffusion model until convergence, (2) using the said Diffusion model for posterior sampling to estimate the expectation. Learning the model until convergence allows using the recent posterior sampling algorithms, which have shown great promise. The experimental setting is also nice and quite diversified; the toy example is welcomed!", "weaknesses": "I see a few weaknesses; \n- *experiments*: for the posterior sampling algorithm, the authors compare to a narrow set of methods (actually only one method but with different covariance matrices). For example, why is there no comparison with DPS, DDRM or even CoPaint [1], which has been shown to exhibit very good performance (and can easily be extended to noisy inverse problems)? Also, are Figure 1 and 2 obtained using the true covariance matrix? Isn't this a misleading comparison? These figures should be compiled using the proposed method and not the ground-truth covariance. Overall, I think that there isn't enough evidence that the proposed posterior sampler is superior to what is proposed in the literature, given its large memory and time complexity. Especially since the current method requires drawing batches of samples during training, which may hardly fit in a modest GPU due to the computation of the Jacobian. \n\n- *identifiability*: doesn't the proposed framework suffers from an identifiability problem, meaning that the learned model may not actually learn clean data but something else? For example, let's assume that $A$ is the half mask matrix, $x$ is a distribution $q$ over images and $p(y|x,A)$ is a dirac delta at Ax. Then the optimal model can be either $q$ or the law of $A^\\dagger X$ with $X \\sim q$, where $A^\\dagger$ is the pseudo-inverse of $A$. Meaning that the model may as well learn images with missing half. Of course in practice the architecture has inductive biases that help avoid this but surely a UNet cannot simply guess the other half of the images while never seeing these during training. I believe that this is a significant drawback that should be addressed. \n\n\n[1] Zhang, Guanhua, et al. \"Towards coherent image inpainting using denoising diffusion implicit models.\" (2023).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the legitimiate concerns you have raised.\n\n* **W1** (Experiments) We follow your suggestion and benchmark MMPS against previous methods (DPS and $\\Pi$GDM). We invite you to consult the global rebuttal regarding these additional experiments.\n\n  Concerning Figures 1 and 2, we use Eq. (20) to compute the covariance $\\mathbb{V}[x \\mid x_t]$, which is tractable in this toy experiment. We do not see how this would be misleading.\n\n  We also note that we never materialize the Jacobian in the experiments of Section 5. Instead, we leverage the CG method to solve the linear system in Eq. (22), which only requires a cheap vector-Jacobian product per iteration.\n\n* **W2** (Identifiability) Indeed, we forgot to explain how our method would behave when $p(x)$ cannot be uniquely identified from the observations. We propose to replace lines 298-300 with the following paragraph\n\n  > Finally, as mentioned in Section 6, empirical Bayes is an ill-posed problem in that distinct prior distributions may result in the same distribution over observations. In other words, it is generally impossible to identify \"the\" ground-truth distribution $p(x)$ given an empirical distribution of observations $p(y)$. Instead, for a sufficiently expressive diffusion model, our EM method will eventually converge to a prior $q_\\theta(x)$ that is consistent with $p(y)$, but generally different from $p(x)$. In future work, we would like to follow the maximum entropy principle, as advocated by Vetter et al. [36], so as not to reject any possible hypothesis.\n\n  We emphasize that this identifiability issue is a limitation of the problem itself, and not of our method.\n\nWe believe that this rebuttal addresses most of your concerns and, therefore, kindly ask you to reconsider your score.", "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:
We would like to thank the reviewers for the quality and pertinence of their reviews. We are glad that all reviewers found the topic of our work interesting and timely.

Reviewers **3yec**, **XeV5**, **mKff** and **tnYw** found the method sound and well presented. Theirs concerns mainly regard the extent of the experiments. Notably, reviewers **3yec**, **mKff** and **tnYw** rightfully comment that the proposed Moment Matching Posterior Sampling (MMPS) method should be benchmarked independently from the context of learning from observations. We propose to address these concerns with additional experiments which we describe below.

Reviewer **s8XN** found the presentation of the method and its pipeline unclear, which prevented them from judging the contribution and results. We thank the reviewer for this opportunity to clarify our manuscript. We describe the relevant changes to the manuscript in reviewer **s8XN**'s rebuttal.

We would also like to emphasize that our work includes two contributions: a novel method to learn diffusion models from observations **and** an improved posterior sampling scheme. Many articles focusing on only one of these contributions have been published in major venues.

**Additional experiments**

* We **repeat the corrupted CIFAR-10 experiment at more corruption levels** (0.25, 0.50 and 0.75). We take the opportunity to refine the exponential moving average (EMA) decay rate of our training step to further improve the results of our method, which we present in Table 1 in the attached PDF. As expected, reducing the corruption level leads to even better final diffusion models. The EM algorithm also converges faster, which is expected as each observation $y$ conveys more information about its latent $x$. We propose to present and discuss these results in the main text.

* We **benchmark MMPS against SOTA posterior sampling methods**, namely DPS [21] and $\Pi$GDM [22], for several linear inverse problems on the FFHQ dataset. For a fair comparison, we adapt the official code published by Chung et al. [21] and use the provided pre-trained diffusion model as diffusion prior. We present premilinary results in the attached PDF. We consider 4 linear inverse problems:
    - Box inpainting with high noise ($\sigma_y = 1$)
    - Random inpainting with low noise ($\sigma_y = 10^{-2}$)
    - Motion deblur with moderate noise ($\sigma_y = 10^{-1}$)
    - Super resolution (4x) with moderate noise ($\sigma_y = 10^{-1}$)

    We find that MMPS requires very few sampling steps to generate qualitative samples and remains remarkably stable for challenging inverse problems (non-diagonal measurement and/or high noise). Conversely, DPS requires many steps to converge and $\Pi$GDM fails for moderate-to-high noise levels. We also find that increasing the number of CG iterations improves image quality/sharpness, but only marginally and with rapidly diminishing returns. Finally, our analysis (see Table 2) shows that an MMPS step is moderately slower (+16ms per CG iteration) than a DPS or $\Pi$GDM step, while only using 10% more memory. This hopefully addresses the concerns of high time and memory complexities raised by reviewers **tnYw** and **mKff**. We propose to present and discuss an extended version of these results in a new appendix section. We will also rewrite lines 291-297 in the discussion to better reflect these results.

Meta-review: This paper considers a class of latent diffusion models where a latent variable $x$ is associated with every observation $y$, and is modeled as coming from a diffusion model $x \sim q_\theta(x)$ with parameters $\theta$. The paper proposes a Monte Carlo EM algorithm to learn $\theta$ wherein the E-step consists of sampling $x$ from the posterior $q_\theta(x \mid y)$ and the M-step consists of updating $\theta$ by denoising score matching. Due to the sensitivity of this procedure the on posterior samples in the E-step, the paper also proposes an improved method (denoted MMPS) for posterior sampling from diffusion models. 

Most reviewers were excited by the paper, finding the methodological contributions timely, the results impressive, and writing to be of high quality.

There was initial concern among some reviewers about the computational  and memory overhead of the proposed procedure; however the authors responded convincingly to these concerns in their rebuttal. For example, reviewers were concerned that the EM algorithm essentially involved fitting a diffusion model from scratch every M-step. However, the authors clarified in the discussion period that the diffusion model is initialized at its last iterate, and further provided supplemental experimental evidence that indicated that the wall-clock time of the proposed method was on the same order as the baseline methods. 

Reviewers also noted that the proposed posterior sampling method was interesting in its own right and deserved separate recognition and independent benchmarking. In response, the authors provided a suite of experimental results that compared MMPS to two previous methods, showing favorable results.

One reviewer also suggested additional comparisons to AmbientDiffusion with different hyperparameter values, which the authors supplied in the response, again showing favorable results.

There were a range of other questions from reviewers (e.g., regarding identifiability, relation to related work, and the impact of various approximations) to which the authors replied thoroughly and convincingly.
Decision: Accept (poster)