Learning Diffusion Priors from Observations by Expectation Maximization

Authors: françois  rozet
gérôme  andry
françois  lanusse
gilles  louppe
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 14:28:57

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.

Main paper: 
 introduction: 
 Many scientific applications can be formalized as Bayesian inference in latent variable models, where the target is the posterior distribution p(x | y) ∝ p(y | x) p(x) given an observation y ∈ R M resulting from a forward process p(y | x) and a prior distribution p(x) over the latent variable x ∈ R N . Notable examples include gravitational lensing inversion [1] [2] [3] , accelerated MRI [4] [5] [6] [7] [8] , unfolding in particle physics [9, 10] , and data assimilation [11] [12] [13] [14] . In all of these examples, the observation y alone is either too incomplete or too noisy to recover the latent x. Additional knowledge in the form of an informative prior p(x) is crucial for valuable inference.Recently, diffusion models [15, 16] proved to be remarkable priors for Bayesian inference, demonstrating both quality and versatility [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] . However, to train a diffusion model for the latent x, one would typically need a large number of latent realizations, which by definition are not 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] [29] [30] [31] offer a solution to the problem of prior specification in latent variable models when only observations y are available. The objective of EB is to find the parameters θ of a prior model q θ (x) for which the evidence distribution q θ (y) = ´p(y | x) q θ (x) dx 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] [33] [34] [35] [36] [37] or simple parametric models [38, 39] .In this work, our goal is to use diffusion models for the prior q θ (x), as 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 on models for which the density q θ (x) or samples x ∼ q θ (x) are differentiable with respect to the parameters θ. Instead, we propose an adaptation of the expectation-maximization [40] [41] [42] [43] [44] algorithm where we alternate between generating samples from the posterior q θ (x | y) and training the prior q θ (x) 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. 
diffusion models: 
 The primary purpose of diffusion models (DMs) [15, 16] , also known as score-based generative models [45, 46] , is to generate plausible data from a distribution p(x) of interest. Formally, adapting the continuous-time formulation of Song et al. [46] , samples x ∈ R N from p(x) are progressively perturbed through a diffusion process expressed as a stochastic differential equation (SDE)EQUATIONwhere f t ∈ R is the drift coefficient, g t ∈ R + is the diffusion coefficient, w t ∈ R N denotes a standard Wiener process and x t ∈ R N is the perturbed sample at time t ∈ [0, 1]. Because the SDE is linear with respect to x t , the perturbation kernel from x to x t is Gaussian and takes the formEQUATIONwhere α t and Σ t = σ 2 t I are derived from f t and g t [46] [47] [48] [49] . Crucially, the forward SDE (1) has an associated family of reverse SDEs [46] [47] [48] [49] EQUATIONwhere η ≥ 0 is a parameter controlling stochasticity. In other words, we can draw noise samples x 1 ∼ p(x 1 ) ≈ N (0, Σ 1 ) and gradually remove the noise therein to obtain data samples x 0 ∼ p(x 0 ) ≈ p(x) by simulating Eq. ( 3) from t = 1 to 0 using an appropriate discretization scheme [16, 45, 46, [49] [50] [51] [52] .In this work, we adopt the variance exploding SDE [45] for which f t = 0 and α t = 1.In practice, the score function ∇ xt log p(x t ) in Eq. ( 3) is unknown, but can be approximated by a neural network trained via denoising score matching [53, 54] . Several equivalent parameterizations and objectives have been proposed for this task [16, 45, 46, [50] [51] [52] . In this work, we adopt the denoiser parameterization d θ (x t , t) and its objective [51] arg min θ E p(x)p(t)p(xt|x) λ t ∥d θ (x t , t) -EQUATIONfor which the optimal denoiser is the mean E[x | x t ] of p(x | x t ). Importantly, E[x | x t ] is linked to the score function through Tweedie's formula [55] [56] [57] [58] ]EQUATIONwhich allows to use s θ (x t ) = Σ -1 t (d θ (x t , t)x t ) as a score estimate in Eq. (3). 
expectation-maximization: 
 The objective of the expectation-maximization (EM) algorithm [40] [41] [42] [43] [44] is to find the parameters θ of a latent variable model q θ (x, y) that maximize the log-evidence log q θ (y) of an observation y. For a distribution of observations p(y), the objective is to maximize the expected log-evidence [43, 44] or, equivalently, to minimize the Kullback-Leibler (KL) divergence between p(y) and q θ (y). That is,EQUATION= arg min θ KL p(y) ∥ q θ (y) .The key idea behind the EM algorithm is that for any two sets of parameters θ a and θ b , we haveEQUATION= E q θ b (x|y) log q θa (x, y) q θ b (x, y)EQUATION≥ E q θ b (x|y) log q θa (x, y)log q θ b (x, y) . (10) This inequality also holds in expectation over p(y). Therefore, starting from arbitrary parameters θ 0 , the EM updateEQUATIONleads to a sequence of parameters θ k for which the expected log-evidence E p(y) log q θ k (y) is monotonically increasing and converges to a local optimum [42] [43] [44] .When the expectation in Eq. ( 12) is intractable, many have proposed to use Monte Carlo approximations instead [59] [60] [61] [62] [63] [64] [65] [66] . Previous approaches include Markov chain Monte Carlo (MCMC) sampling, importance sampling, rejection sampling and their variations [63] [64] [65] [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 a recent review of MCEM algorithms. 
methods: 
 Although rarely mentioned in the literature, the expectation-maximization 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 q θ (y). In the empirical Bayes setting, the forward model p(y | x) is known and only the parameters of the prior q θ (x) should be optimized. In this case, Eq. ( 12) becomesEQUATION= arg maxEQUATION= arg minEQUATIONwhere π k (x) = ´qθ k (x | y) p(y) dy. Intuitively, π k (x) and therefore q θ k+1 (x) assign more density to latents x which are consistent with observations y ∼ p(y) than q θ k (x). In this work, we consider a special case of the empirical Bayes problem where each observation y has an associated measurement matrix A and the forward process takes a linear Gaussian form p(y | x, A) = N (y | Ax, Σ y ). This formulation allows the forward process to be potentially different for each observation y. 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) ∼ p(y, A). As a result, π k (x) in Eq. ( 15) becomes π k (x) = ´qθ k (x | y, A) p(y, A) dy.Now that our goals and assumptions are set, we present our method to learn a diffusion model q θ (x) for the latent x from observations y by expectation-maximization. The idea is to decompose Eq. ( 15) into (i) generating a dataset of i.i.d. samples from π k (x) and (ii) training q θ k+1 (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 π k (x), we first sample a pair (y, A) ∼ p(y, A) and then generate x ∼ q θ k (x | y, A) from the posterior. Unlike previous MCEM algorithms that rely on expensive and hard to tune sampling strategies [63] [64] [65] [66] , the use of a diffusion model enables efficient and embarrassingly parallelizable posterior sampling [21] [22] [23] . However, the quality of posterior samples is critical for the EM algorithm to converge properly [63] [64] [65] [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 q θ (x) by a denoiser network d θ (x t , t) 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 θ k .Initialization An important part of our pipeline is the initial prior q 0 (x). Any initial prior leads to a local optimum [42] [43] [44] , but an informed initial prior can reduce the number of iterations until convergence. In our experiments, we take a Gaussian distribution N (x | µ x , Σ x ) as initial prior and fit its mean and covariance by -you guessed it! -expectation-maximization. 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-training 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, if the initial prior q 0 (x) does not cover latents that are otherwise plausible under the observations, the following priors q θ k (x) may not cover these latents either. A conservative initial prior is therefore preferable for scientific applications.To sample from the posterior distribution q θ (x | y) ∝ q θ (x) p(y | x) of our diffusion prior q θ (x), we have to estimate the posterior score ∇ xt log q θ (x t | 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 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 θ (x).Diffusion posterior sampling Thanks to Bayes' rule, the posterior score ∇ xt log p(x t | y) can be decomposed into two terms [17, 18, [21] [22] [23] [24] [25] 46 ] 17) is that, when σ t is small, assuming that p(x | x t ) is narrowly concentrated around its mean E[x | x t ] is reasonable. However, this approximation is very poor when σ t is not negligible. Consequently, DPS [21] is unstable, does not properly cover the support of the posterior p(x | y) and often leads to samples x which are inconsistent with the observation y [22] [23] [24] [25] .EQUATIONMoment matching To address these flaws, following studies [22] [23] [24] [25] take the covariance V[x | x t ] into account when estimating the likelihood score ∇ xt log p(y | x t ). Specifically, they consider the Gaussian approximationEQUATIONwhich is closest to p(x | x t ) in Kullback-Leibler (KL) divergence [68] . Then, assuming a linear Gaussian forward process p(y | x) = N (y | Ax, Σ y ), we obtain [68] q(yEQUATIONwhich allows to estimate the likelihood score ∇ xt log p(y | x t ) asEQUATIONunder the assumption that the derivative of V[x | x t ] with respect to x t is negligible [24, 25] . Even with simple heuristics for V[x | x t ], such as Σ t [20] or (Σ -1 t + Σ -1 x ) -1 [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 q(x t | y) ∝ p(x t ) q(y | x t ) in the presence of local covariances -i.e. covariances in the neighborhood of x t . 10 -1 [69] between the posteriors p(x t | y) and q(x t | y) for different heuristics of V[x | x t ] when the prior p(x) lies on 1-d manifolds embedded in R 3 . Lines and shades represent the 25-50-75 percentiles for 64 randomly generated manifolds [71] and measurement matrices A ∈ R 1×3 . Using V[x | x t ] instead of heuristics leads to orders of magnitude more accurate posteriors q(x t | y).10 1 10 3 S λ (p(x t | y) q(x t | y)) V[x|xt] (Σ -1 x + Σ -1 t ) -1 Σt 0 Figure 2. Sinkhorn divergenceWe illustrate this behavior in Figure 1 and measure its impact as the Sinkhorn divergence [69, 70] between the posteriors p(x t | y) and q(x t | y) when the prior p(x) lies on randomly generated 1-dimensional manifolds [71] embedded in R 3 . The prior p(x) is modeled as a mixture of isotropic Gaussians centered around points of the manifold, which gives access to p(x t ), E[x | x t ] and V[x | x t ] analytically. The results, presented in Figure 2 , indicate that using V[x | x t ] instead of heuristics leads to orders of magnitude more accurate posteriors q(x t | y). We expect this gap to further increase with real high-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] [64] [65] [66] , we choose to estimate V[x | x t ] using Tweedie's covariance formula [55] [56] [57] [58] VEQUATIONConjugate gradient method As noted by Finzi et al. [24] , explicitly computing and materializing the Jacobian 20) would still be intractable. Fortunately, we observe that the matrix Σ y + AV[x | x t ]A ⊤ 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 M v = 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 M v given a vector v. In our case, the linear system to solve is∇ ⊤ xt d θ (x t , t) ∈ R N ×N is intractable in high dimension. Furthermore, even if we had access to V[x | x t ], naively computing the inverse of the matrix Σ y + AV[x | x t ]A ⊤ in Eq. (EQUATION)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∇ ⊤ xt d θ (x t , t) ≈ ∇ ⊤ xt E[x | x t ]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 to very few iterations (1 to 3) already leads to significant improvements over using heuristics for the covariance V[x | x t ]. Alternatively, the CG method can be replaced by other iterative algorithms that can solve non-symmetric non-definite linear systems, such as GMRES [75] or BiCGSTAB [76] , at the cost of slower convergence. 
results: 
 We conduct three experiments to demonstrate the effectiveness of our method. We design the first experiment around a low-dimensional latent variable x whose ground-truth distribution p(x) is known. In this setting, we can use asymptotically exact sampling schemes such as predictor-corrector sampling [23, 46] or twisted diffusion sampling [77] without worrying about their computational cost. This allows us to validate our expectation-maximization pipeline (see Algorithm 1) in the limit of (almost) exact posterior sampling. The remaining experiments target two benchmarks from previous studies: corrupted CIFAR-10 and accelerated MRI. These tasks provide a good understanding of how our method would perform in typical empirical Bayes applications with limited data and compute.In this experiment, the latent variable x ∈ R 5 ∼ p(x) lies on a random 1-dimensional manifold embedded in R 5 represented in Figure 7 . Each observation y ∈ R 2 ∼ N (y | Ax, Σ y ) is the result of a random linear projection of a latent x plus isotropic Gaussian noise (Σ y = 10 -4 I). The rows of the measurement matrix A ∈ R 2×5 are drawn uniformly from the unit sphere S 4 . We note that observing all push-forward distributions p(u ⊤ x) with u ∈ S N -1 of a distribution p(x) in R N is sufficient to recover p(x) in theory [78, 79] . In practice, we generate a finite training set of 2 16 pairs (y, A).We train a DM q θ (x) parameterized by a multi-layer perceptron d θ (x t , t) for K = 32 EM iterations. We apply Algorithm 3 to estimate the posterior score ∇ xt log q θ (x t | y, A), but rely on the predictorcorrector [23, 46] sampling scheme with a large number (4096) of correction steps to sample from the posterior q θ (x | y, A). We provide additional details such as noise schedule, network architectures, and learning rate in Appendix C.As expected, the model q θ k (x) converges towards a stationary distribution whose marginals are close to the marginals of the ground-truth p(x), as illustrated in Figure 3 . We attribute the remaining artifacts to finite data and inaccuracies in our sampling scheme. -2 0 2 q k (x 0 , x 1 ) -2 0 2 q k (x 1 , x 2 ) -2 0 2 k = 1 -2 0 2 q k (x 2 , x 3 ) -2 0 2 k = 2 -2 0 2 k = 4 -2 0 2 k = 8 -2 0 2 k = 16 -2 0 2 k = 32(I + Σ -1 t ) -1 (75%) Σt (75%)k = 1 k = 2 k = 4 k = 8 k = 16 k = 32Following Daras et al. [80] , we take the 50 000 training images of the CIFAR-10 [81] dataset as latent variables x. A single observation y is generated for each image x by randomly deleting pixels with probability ρ. The measurement matrix A is therefore a binary diagonal matrix. We add negligible isotropic Gaussian noise (Σ y = 10 -6 I) for fair comparison with AmbientDiffusion [80] which cannot handle noisy observations. For each corruption rate ρ ∈ {0.25, 0.5, 0.75}, we train a DM q θ (x) parameterized by a U-Net [82] inspired network d θ (x t , t) for K = 32 EM iterations. We apply Algorithm 2 with T = 256 discretization steps and η = 1 to approximately sample from the posterior q θ (x | y, A). We apply Algorithm 3 with several heuristics for V[x | x t ] to compare their results against Tweedie's covariance formula. For the latter, we truncate the conjugate gradient method in Algorithm 4 to a single iteration.For each model q θ k (x), we generate a 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 in Table 1 and Figures 4 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 | x t ] leads to poor sample quality.Magnetic resonance imaging (MRI) is a non-invasive medical imaging technique used in radiology to inspect the internal anatomy and physiology of the body. MRI measurements of an object are obtained in the frequency domain, also called k-space, using strong magnetic fields. However, measuring the In this experiment, following Kawar et 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 x, represented as a 320 × 320 gray-scale image. Scans are min-max normalized such that pixel values range between -2 and 2. A single observation y is generated for each slice x by first applying the discrete Fourier transform and then a random horizontal frequency sub-sampling with acceleration factor R = 6 [85, 86] , meaning that a proportion 1 /R of all frequencies are observed on average. Eventually, we obtain 24 853 k-space observations to which we add isotropic Gaussian noise (Σ y = 10 -4 I) to match Kawar et al. [85] .Once again, we train a DM q θ (x) parameterized by a U-Net [82] inspired network d θ (x t , t) for K = 16 EM iterations. We apply Algorithm 2 with T = 64 discretization steps and η = 1 to approximately sample from the posterior q θ (x | y, A) and truncate the conjugate gradient method in Algorithm 4 to 3 iterations. After training, we employ our final model q θ K (x) as a diffusion prior for accelerated MRI. Thanks to our moment matching posterior sampling, we are able to infer plausible scans at acceleration factors up to R = 32, as shown in Figure 6 . Our samples are noticeably more detailed than the ones of Kawar et al. [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 . 
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 a prior distribution that reproduces the empirical distribution of observations when pushed through a nondifferentiable black-box forward process. Dockhorn et 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(y | 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-Wasserstein distance [78, 79] with the empirical distribution of observations. These methods rely on generative models q θ (x) for which the density q θ (x) or samples x ∼ q θ (x) are differentiable with respect to the parameters θ, which is not or hardly the case for diffusion models.Closer to this work, Daras et al. [80] and Kawar et al. [85] [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 Σ y = σ 2 y I) data by enforcing the "consistency" of the denoiser along diffusion paths. They prove that the mean E[x | x t ] is the unique consistent denoiser. Interestingly, this training algorithm also relies 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 unconditional 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 ∇ xt log p(y | x t ) when the forward process p(y | x) is Gaussian. Song et al. [22] and Rozet et al. [23] improve this approximation by taking the covariance V[x | x t ] 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 | x t ] based on Tweedie's covariance formula [55] [56] [57] [58] . Finzi et al. [24] take the same approach, but materialize the matrix AV[x | x t ]A ⊤ which is intractable in high dimension. Boys et al. [25] replace the covariance V[x | x t ] with a row-sum approximation diag(e ⊤ V[x | x t ]) where e is the all-ones vector. This approximation is only valid when V[x | x t ] is diagonal, which limits its applicability. Instead, we take advantage of the conjugate gradient method [72] to avoid materializing AV[x | x t ]A ⊤ . A potential cheaper solution is to train a sparse approximation of V[x | x t ], as proposed by Peng et al. [93] , but this approach is less general and not immediately applicable to any diffusion model.A parallel line of work [94] [95] [96] extends the conditioning of diffusion models to arbitrary loss terms ℓ(x, y) ∝log p(y | x), for which no reliable analytical approximation of the likelihood score ∇ xt log p(y | x t ) exists. Song et al. [94] rely on Monte Carlo approximations of the likelihood p(y | x t ) = ´p(y | x) p(x | x t ) dx by sampling from a Gaussian approximation of p(x | x t ). Conversely, He et al. [96] use the mean E[x | x t ] as a point estimate for p(x | x t ), but leverage a pre-trained encoder-decoder pair to project the updates of x t within its manifold. We note that our use of the covariance V[x | x t ] similarly leads to updates tangent to the manifold of x t .Finally, Wu et al. [77] propose a particle-based posterior sampling scheme that is asymptotically exact in the limit of infinitely many particles as long as the likelihood approximation q(y | x t ) -here named the twisting function -converges to p(y | x) as t approaches 0. Using TDS [77] as part of our expectation-maximization 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. 
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] [29] [30] [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] [34] [35] [36] [37] . In addition, as mentionned in Section 6, these EB methods are not applicable to diffusion models. As such, we choose to benchmark our work against similar methods in the diffusion model literature [80, 85] , but stress that a proper comparison with previous empirical Bayes methods would be valuable for both communities. We also note that Monte Carlo EM [59] [60] [61] [62] [63] [64] [65] [66] can handle arbitrary forward processes p(y | x) as long as one can sample from the posterior q θ (x | 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.From a computational perspective, the iterative nature of our expectation-maximization method is a 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 q θ (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 a 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 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 θ (x) that is consistent with p(y), but generally different from p(x). Following the maximum entropy principle, as advocated by Vetter et al. [37] , is left to future work.Algorithm 1 Expectation-maximization pipeline 1 Choose an initial prior q 0 (x) 2 Initialize the parameters θ of the denoiser d θ (x t , t) 3 for k = 1 to K do 4 for i = 1 to S do 5 y i , A i ∼ p(y, A) 6 x i ∼ q k-1 (x | y i , A i ) # Posterior sampling 7 repeat 8 i ∼ U({1, . . . , S}) 9 t ∼ U(0, 1) 10 z ∼ N (0, I) 11 x t ← x i + σ t z 12 ℓ ← λ t d θ (x t , t) -x i 2 # Denoising score matching 13 θ ← GRADIENTDESCENT(θ, ∇ θ ℓ) 14 until convergence 15 θ k ← θ 16 q k := q θ k 17 return θ K Algorithm 2 DDIM-style posterior sampling 1 x 1 ∼ N (0, Σ 1 ) 2 for i = T to 1 do 3 s ← i-1 /T 4 t ← i /T 5 x ← x t + Σ t s θ (x t | y, A) # Estimate E[x | x t , y, A] 6 z ∼ N (0, I) 7 x s ← x + σ s 1 -η 1 - σ 2 s σ 2 t x t - x σ t + σ s η 1 - σ 2 s σ 2 t z 8 return x 0 Algorithm 3 Moment matching posterior score 1 function s θ (x t | y, A) # Estimate ∇ xt log q θ (x t | y, A) 2 x ← d θ (x t , t) 3 if Tweedie then 4 Σ x|xt ← Σ t ∇ xt d θ (x t , t) ⊤ 5 else 6 Σ x|xt ← Σ t or (I + Σ -1 t ) -1 or (Σ -1 x + Σ -1 t ) -1 7 u ← Σ y + AΣ x|xt A ⊤ -1 (y -Ax) # Solve with conjugate gradient method 8 s y|x ← ∇ xt d θ (x t , t) ⊤ A ⊤ u # Estimate ∇ xt log q θ (y | x t , A) 9 s x ← Σ -1 t (x -x t ) # Estimate ∇ xt log q θ (xEQUATIONWe provide proofs of Theorem 1 for completeness, even though it is a well known result [55] [56] [57] [58] .Proof.∇ xt log p(x t ) = 1 p(x t ) ∇ xt p(x t ) = 1 p(x t ) ˆ∇xt p(x, x t ) dx = 1 p(x t ) ˆp(x, x t ) ∇ xt log p(x, x t ) dx = ˆp(x | x t ) ∇ xt log p(x t | x) dx = ˆp(x | x t ) Σ -1 t (x -x t ) dx = Σ -1 t E[x | x t ] -Σ -1 t x t Proof. ∇ 2 xt log p(x t ) = ∇ xt ∇ ⊤ xt log p(x t ) = ∇ xt E[x | x t ] ⊤ Σ -1 t -Σ -1 t = ˆ∇xt p(x | x t ) x ⊤ dx Σ -1 t -Σ -1 t = ˆp(x | x t )∇ xt log p(x t | x) p(x t ) x ⊤ dx Σ -1 t -Σ -1 t = ˆp(x | x t ) Σ -1 t x -E[x | x t ] x ⊤ dx Σ -1 t -Σ -1 t = Σ -1 t E xx ⊤ | x t -E[x | x t ] E[x | x t ] ⊤ Σ -1 t -Σ -1 t = Σ -1 t V[x | x t ] Σ -1 t -Σ -1 tAll experiments are implemented within the JAX [73] automatic differentiation framework.The code for all experiments is made available at https://github.com/francois-rozet/ diffusion-priors.Diffusion models As mentioned in Section 2, in this work, we adopt the variance exploding SDE [45] and the denoiser parameterization [51] . Following Karras et al. [51] , we precondition our denoiser d θ (x t , t) asEQUATION)where h θ (x, log σ) is an arbitrary noise-conditional network. The scalar log σ is embedded as a vector using a sinusoidal positional encoding [100] . In our experiments, we use an exponential noise schedule σ t = exp (1t) log 10 We apply Algorithm 3 to estimate the posterior score ∇ xt log p(x t ) and truncate Algorithm 4 to 3 iterations. We rely on the predictor-corrector [23, 46] sampling scheme to sample from the posterior q θ (x | 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 1 correction step. At each EM iteration, we generate a single latent x for each pair (y, 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 λ = 1e -3.Corrupted CIFAR-10 The noise-conditional network h θ (x, log σ) is a U-Net [82] with residual blocks [104] , SiLU [101] activation functions and layer normalization [102] . Each residual block is modulated with respect to the noise σ t in the style of diffusion transformers [105] . A multi-head self-attention block [100] is inserted after each residual block at the last level of the U-Net. We train (I + Σ -1 t ) -1 is used for V[x | x t ].The generated images become gradually more detailed but some noise remains. 11 ), but do not present artifacts typical to unresolved frequencies (see Figure 12 ). (I + Σ -1 t ) -1 is used for V[x | x t ].The samples start to present vertical artifacts due to poor sampling. Figure 16 . Example of samples from the model q θ k (x) after k = 4 EM iterations for the accelerated MRI experiment when the heuristic(I + Σ -1 t ) -1 is used for V[x | x t ].The artifacts introduced by the poor sampling get amplified at each iteration, leading to a total collapse after few iterations.In this section, we evaluate the moment matching posterior sampling (MMPS) method presented in 4.2 independently from the context of learning from observations. The code for this section is made available at https://github.com/francois-rozet/mmps-benchmark.Tasks We consider four linear inverse problems on the 256 × 256 FFHQ [99] dataset. (i) For box inpainting, we mask out a randomly positioned 128 × 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 (σ y = 0.01). (iii) For motion deblur, we apply a randomly generated 61 × 61 motion blur kernel and add a medium amount of noise (σ y = 0.1). (iv) For super resolution, we apply a 4× bicubic downsampling and add a medium amount of noise (σ y = 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 DiffPIR [26] , ΠGDM [22] and TMPD [25] . We use the DDIM [50] sampler with η = 1 for all methods, which is equivalent to the DDPM [16] sampler. We fine-tune the hyperparameters of DPS (ζ ′ = 0.5) and DiffPIR (λ = 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 σ t , 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 images1 of the FFHQ [99] dataset. We generate a sample for each observation with all considered posterior sampling methods. All methods are executed with the same random seed. We compute three standard image reconstruction 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 As a side note, we 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 as the mean E[x | 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 simultaneously 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 1 and 3 iterations is much larger than between 3 and 5. MMPS is also remarkably stable with respect to the number of sampling steps in contrast to DPS [21] , DiffPIR [26] and ΠGDM [22] which are sensitive to the number of steps and choice of hyperparameters. Finally, MMPS requires fewer sampling steps to reach the same image quality as previous methods, which largely makes up for its slightly higher step cost. [22] and TMPD [25] .Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The code for all experiments is made available at https://github.com/ francois-rozet/diffusion-priors. Instructions to acquire the datasets are provided. Guidelines:• The answer NA means that paper does not include experiments requiring code. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Justification: There is no societal impact of the work performed.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a 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 or 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). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 to usage guidelines or 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 a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes] Justification: We acknowledge the use of the NYU fastMRI dataset [7, 8] in the acknowledgments section and follow its terms of use.Guidelines:• The answer NA means that the paper does not 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 a 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.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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 is not available online, the authors are encouraged to reach out to the asset's creators. 
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/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was 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. 
crowdsourcing 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 as 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 in the supplemental material is fine, but if the main contribution 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. 
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 an equivalent approval/review based on the requirements of your country or institution) were obtained?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. • Depending on 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.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to 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.Chung et al. [21] do not indicate which subset of 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. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: The paper has also several weaknesses.

1) 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.
2) 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.
3) 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.
4) 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.
5) 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.
6) 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.
7) 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.
8) 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.
questions: questions: See weaknesses above. I would be happy to further increase my score if my concerns are properly addressed.
limitations: limitations: The authors have adequately addressed the limitations of their method.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
2. The validity of the approximation needs to be shown theoretically. See the "Questions" section below.
questions: 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?
2. 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?
3. 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.
4. Are there any experiments on the improvement per iteration of the Conjugate Gradient method?
limitations: limitations: Yes, the authors adequately addressed the limitations.
paper_summary: paper_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.
strength_weakness: strength_weakness: The paper is poorly written, the method is not clear, derivations are incomplete. See my detailed review in the questions part.
questions: 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.

1) 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.

2) The authors put a prior on the matrix $A$. What is this prior?

3) 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?

4) 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.

5) 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.

6) 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?

7) 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?

I strongly suggest authors to clarify the pipeline clearly and precisely.
limitations: limitations: Yes
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 
2. 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. 
3. 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?
4. Does Equation (22) hold for any vector $v$ as the equation suggests?
5. 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).
6. Missing discussion and comparison with other baselines (e.g. [82]) highlighted in the related works Section 6.
7. 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. 

Missing Related Works

[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)

[2] [Improving Diffusion Models for Inverse Problems
Using Optimal Posterior Covariance](https://arxiv.org/pdf/2402.02149)
questions: questions: Please refer to the weaknesses above.
limitations: limitations: Yes, the authors have adequately addressed the limitations.
paper_summary: paper_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.
strength_weakness: strength_weakness: I see a few weaknesses; 
- *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. 

- *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. 


[1] Zhang, Guanhua, et al. "Towards coherent image inpainting using denoising diffusion implicit models." (2023).
questions: questions: N/A
limitations: limitations: see above

Comments: title: Rebuttal acknowledgement
comment: Here are preliminary quantitative results. We consider the same linear inverse problems as previously and add two more baselines (TMPD and DiffPIR). We consider 3 standard image quality metrics (LPIPS, PSNR and SSIM) which we average over 100 FFHQ observations. MMPS outperforms all baseline consistently. We note that PSNR and SSIM are biased towards smooth samples while LPIPS favors perceptually detailed samples.

|  |  |  |  | Box inpainting |  |  | Random inpainting |  |  | Motion deblur |  |  | Super resolution |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Method | Steps | LPIPS $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ |
| DiffPIR | 10 | 0.97 | 14.23 | 0.13 | 0.78 | 10.97 | 0.32 | 0.39 | 23.98 | 0.51 | 0.31 | 26.45 | 0.69 |
| DiffPIR | 100 | 0.88 | 15.03 | 0.20 | 0.68 | 10.26 | 0.25 | 0.17 | 25.17 | 0.71 | 0.19 | 26.30 | 0.73 |
| DiffPIR | 1000 | 0.90 | 14.65 | 0.18 | 0.74 | 9.51 | 0.21 | 0.17 | 24.67 | 0.68 | 0.19 | 25.67 | 0.68 |
| DPS | 10 | 0.64 | 10.41 | 0.34 | 0.58 | 12.68 | 0.43 | 0.75 | 8.63 | 0.27 | 0.58 | 12.01 | 0.41 |
| DPS | 100 | 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 | 1000 | 0.22 | 21.02 | 0.64 | 0.19 | 21.90 | 0.66 | 0.18 | 22.91 | 0.66 | 0.16 | 25.02 | 0.72 |
| PGDM | 10 | 0.48 | 15.32 | 0.57 | 0.52 | 14.44 | 0.47 | 0.46 | 18.61 | 0.54 | 0.37 | 20.47 | 0.67 |
| PGDM | 100 | 0.92 | 14.18 | 0.22 | 0.23 | 21.59 | 0.66 | 0.84 | 17.02 | 0.34 | 0.45 | 21.80 | 0.51 |
| PGDM | 1000 | 0.83 | 14.38 | 0.34 | **0.14** | 22.47 | 0.70 | 0.88 | 16.63 | 0.34 | 0.51 | 20.88 | 0.48 |
| TMPD | 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 | 100 | 0.27 | 19.86 | 0.64 | 0.58 | 10.73 | 0.31 | 0.17 | _26.22_ | 0.76 | 0.17 | 26.79 | 0.77 |
| TMPD | 1000 | 0.25 | 19.53 | 0.62 | 0.68 | 9.98 | 0.25 | 0.14 | 25.91 | 0.74 | 0.14 | 26.53 | 0.76 |
| MMPS (1) | 10 | 0.27 | 21.19 | _0.68_ | 0.26 | 22.41 | 0.69 | 0.33 | 22.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 (1) | 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 (5) | 10 | 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 |
Thank you for taking the time and effort to review our manuscript. Your constructive feedback is deeply appreciated. We will make sure to evaluate MMPS quantitatively and qualitatively against more posterior sampling methods in the camera-ready version.
Thank you for your response. I understand that when the denoiser is learned optimally, the covariance is also obtained exactly. My point is still that your approximation relies on taking the jacobian of the denoiser, which is not necessarily a good approximation of the true jacobian of the denoiser (a good parametric $f_\theta$ approximation of a function $f$ does not guarantee that $\nabla_x f_\theta$ is a good parametric approximation of $\nabla_x f$). It would be better to compare with the covariance approximation against the results obtained with the perfect Gaussian approximation. 

The paper has two contributions; the first one is the EM based algorithm for training a DM with incomplete data and the second contribution is the novel posterior sampler. While the benchmarks for the EM based algorithm are good in my opinion, the ones for the new posterior sampler are rather weak, even though the authors have compared with DPS and PGDM in the rebuttal. Still, this is not enough in my opinion since the comparisons provided are entirely qualitative and are not very convincing. More quantitative benchmarks are required, with comparisons against recent methods that use/do not use the Jacobian of the denoiser. Example: Diffpir [1] or DDNM [2]. 

I have increased my score 7 and I hope that the authors will do their best to strengthen the arguments about their posterior sampler. 

[1]Zhu, Yuanzhi, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool. "Denoising diffusion models for plug-and-play image restoration." In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 1219-1229. 2023.

[2] Wang, Yinhuai, Jiwen Yu, and Jian Zhang. "Zero-shot image restoration using denoising diffusion null-space model." arXiv preprint arXiv:2212.00490 (2022).
You raised two key weaknesses in your review and the authors addressed both in their rebuttal. Can you provide more detail about how their rebuttal affected your opinion about the paper? The authors asked what more you would have wanted to see to change your score; please respond to this.
The authors believe their rebuttal addresses most of your concerns; is that true? If not, please say what aspects of their rebuttal are insufficient.
The authors have provided a detailed response to your questions. How has it affected your opinion about the paper? The authors believe their rebuttal addresses most of your concerns; is that true?
Thank you again for taking the time and effort to review our manuscript. Your constructive feedback and recognition of our work are deeply appreciated. We will make sure to include the additional discussion and experiments in the camera-ready version.
Thank you for your answer.

> Regarding Figures 1 and 2 I do not see what is their point since you use the true covariance matrix. A more convincing should include the approximation obtained using the Jacobian of the denoiser and then showing that with this approximation one can get decent results, comparable to what you obtain using the true covariance matrix.

Thank you for clarifying your concern here. When the denoiser is trained optimally, that is when $d_\theta(x_t, t) = \mathbb{E}[x \mid x_t]$, there is no approximation in Eq. (21). Tweedie's formula using the Jacobian of the denoiser gives the covariance matrix exactly. Instead of training a denoiser for this toy problem, we assume an optimal denoiser, which ensures that the results are not biased by a choice of parameterization.

> This is only a minor weakness however. I leave my score unchanged.

Given that we have addressed your concerns, we are genuinely surprised by your decision to keep your rating unchanged (weak accept). We would greatly appreciate it if you could describe how to improve our submission.
Thank you for your response and for running the additional experiments. 

Regarding Figures 1 and 2 I do not see what is their point since you use the true covariance matrix. A more convincing should include the approximation obtained using the Jacobian of the denoiser and then showing that with this approximation one can get decent results, comparable to what you obtain using the true covariance matrix. The current plots tell us nothing about how the approximation< error impacts the final result. 

This is only a minor weakness however. I leave my score unchanged.
I would like to thank the authors for their efforts in the rebuttal. Indeed, most of my concerns are addressed. I highly encourage the authors to include this discussion and the additional experiments in the camera-ready version of their work.

Regarding W4, please also include corruption levels above 0.75 in your camera-ready version and experiments on other datasets beyond CIFAR-10. It is crucial to have a holistic evaluation of the method in the same setting as prior work so that we can make progress and set up a nice benchmarking environment for future works. I believe this is very important for the field and I highly encourage the authors to include these additional experiments.

Also regarding W4, for the MRI, I was not pointing to the Ambient Diffusion paper, but to the paper [Ambient Diffusion Posterior Sampling: Solving Inverse Problems with Diffusion Models Trained on Corrupted Data](https://arxiv.org/abs/2403.08728). This paper looks at exactly the same setting as yours and hence it should be cited and benchmarked.

Regarding W6, Consistent Diffusion Meets Tweedie provides a general algorithm to train diffusion models from noisy data. In the context of the paper, it seems that this method is evaluated for fine-tuning a pre-trained model to a different dataset, but in principle, it should be doable to use it to train from scratch.

I am increasing my score to 8. I hope that the authors will include the additional discussion and the experiments in their camera-ready and I am looking forward to reading this version.