Posterior sampling via Langevin dynamics based on generative priors

## Abstract

Abstract Posterior sampling in high-dimensional spaces using generative models holds significant promise for various applications, including but not limited to inverse problems and guided generation tasks. Despite many recent developments, generating diverse posterior samples remains a challenge, as existing methods require restarting the entire generative process for each new sample, making the procedure computationally expensive. In this work, we propose efficient posterior sampling by simulating Langevin dynamics in the noise space of a pre-trained generative model. By exploiting the mapping between the noise and data spaces which can be provided by distilled flows or consistency models, our method enables seamless exploration of the posterior without the need to re-run the full sampling chain, drastically reducing computational overhead. Theoretically, we prove a guarantee for the proposed noise-space Langevin dynamics to approximate the posterior, assuming that the generative model sufficiently approximates the prior distribution. Our framework is experimentally validated on image restoration tasks involving noisy linear and nonlinear forward operators applied to LSUN-Bedroom (256 x 256) and ImageNet (64 x 64) datasets. The results demonstrate that our approach generates high-fidelity samples with enhanced semantic diversity even under a limited number of function evaluations, offering superior efficiency and performance compared to existing diffusion-based posterior sampling techniques.

## 1 Introduction

Generative models that approximate complex data priors have been leveraged for a range of guided generation tasks in recent years [, ]. Early works focused on conditional synthesis using Generative Adversarial Networks (GANs) [, , , , ]. However, diffusion models have recently surpassed GANs as the state of the art in generative modeling [, ], demonstrating superior performance in guided generation tasks [, , ]. Posterior sampling, as a guided generation framework, has garnered significant interest [, , ], particularly for providing candidate solutions to noisy inverse problems.

Solving noisy inverse problems involves reconstructing an unknown signal $x$ from noisy measurements $y$, where the forward model is characterized by the measurement likelihood $p(y\,|\,x)$. The objective is to sample from the posterior distribution $p(x\,|\,y)\propto p(y\,|\,x)p(x)$. Such posteriors are often intractable in practical applications due to the complexity of the prior distribution $p(x)$. However, learned generative models that approximate complex data priors can enable approximate sampling from the posterior $p(x\,|\,y)$. Early approaches for solving inverse problems using diffusion models to approximate $p(x\,|\,y)$ relied on problem-specific architectures [, , ] and required training dedicated generative models for each task [, ]. In contrast, methods that utilize pre-trained diffusion models as priors for posterior sampling offer greater flexibility and are training-free [, , , , ], with recent extensions targeting nonlinear inverse problems [, , , ].

In the context of inverse problems, existing methods can be broadly categorized based on whether they yield a point estimate or multiple estimates. Existing approaches for posterior sampling focus primarily on providing point estimate solutions [, , , ], lacking the ability to generate a diverse set of posterior samples efficiently. For instance, a prominent method, Diffusion Posterior Sampling (DPS) [], predominantly produces point estimates for both linear and non-linear inverse problems. DPS leverages the prior $p(x)$ from a diffusion model and employs multiple denoising steps to transform isotropic Gaussian noise into a desired image, guided by the observations $y$. Generating posterior samples using this method requires re-running the entire sampling process using unique instantiations of Gaussian noise, which is computationally prohibitive and inefficient. Therefore, an algorithm that efficiently accumulates samples from the posterior is desirable.

Figure: Figure 1: (Left) : A schematic representation of posterior sampling via Langevin dynamics in our proposed framework. The sampling process begins with an initial sample $x_{1}^{(0)}$ from the noise space and maps to data space as $x_{0}^{(0)}$ using a deterministic mapper $\Phi$ and progressively updates the noise space input to obtain diverse posterior samples. (Right): Posterior samples generated by our method and DPS-DM. Our approach exhibits higher perceptual diversity, capturing variations in high-level features such as lighting, window style, and wall patterns. Uncertain semantic features are highlighted by red boxes, while persistent properties are shown by green boxes.
Refer to caption: extracted/5896580/figure_1.png

Figure: Figure 2: Reconstruction time comparison between DPS-DM and our method for varying numbers of posterior samples. DPS-DM scales poorly with the number of samples, while our method maintains a nearly constant time, demonstrating significantly lower computational cost. The corresponding Number of Function Evaluations (NFEs) (including NFEs for the warmup stage, refer to Section 5) values per image are annotated.
Refer to caption: extracted/5896580/figure_2_plot.png

In this work, we propose an efficient framework for posterior sampling by modeling it as an exploration process in the noise space of a pre-trained generative model. Specifically, we leverage measurements from the inverse problem to guide the initialization of the noise space, ensuring a more targeted exploration. For sampling, we employ Langevin dynamics directly within the noise space, taking advantage of the one-to-one mapping between noise and data spaces provided by models such as consistency models []. This deterministic mapping eliminates the need for approximating the measurement likelihood, and we establish a theoretical bound on the approximation error for posterior sampling.

Sampling in the noise space allows for a progressive accumulation of posterior samples, enabling efficient exploration and resulting in a diverse set of reconstructions, as demonstrated in Figure 1. Furthermore, Figure 2 illustrates the comparison of reconstruction times between our approach and DPS when generating different numbers of posterior samples per image. While the reconstruction time for DPS increases rapidly with the number of samples, our method incurs only a negligible increase, highlighting its computational efficiency. The key contributions of this work are summarized as follows:

- •
We present a posterior sampling method defined by Langevin dynamics in the noise space of a pre-trained generative model, enabling efficient accumulation of samples.
- •
We provide a theoretical guarantee on the posterior sampling approximation error, which is bounded by the approximation error of the prior by the pre-trained generative model.
- •
Our efficient accumulation of posterior samples facilitates exploration of the posterior, yielding high-fidelity and diverse samples. In experiments, we achieve comparable fidelity to diffusion model posterior sampling methods with superior sample diversity.

#### Notation.

We use $\propto$ to stand for the expression of a probability density up to a normalizing constant to enforce integral one,
e.g. $p(x)\propto F(x)$ means that $p(x)=F(x)/Z$ where $Z=\int F(x)dx$.
For a mapping $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and a distribution $P$, $T_{\#}P$ stands for the push-forwarded distribution,
that is $T_{\#}P(A)=P(T^{-1}A)$ for any measurable set $A$. When both $P$ and $T_{\#}P$ has density, $dP=pdx$, we also use $T_{\#}p$ to denote the density of $T_{\#}P$.

## 2 Background

#### Diffusion models.

Sampling from diffusion models (DMs) is accomplished via simulation of the reverse process corresponding to the forward-time, noising stochastic differential equation (SDE)
$dx_{t}=\mu(x_{t},t)dt+\beta(t)dW_{t}$ [],
where $W_{t}$ is the standard Brownian motion in $\mathbb{R}^{d}$ and $t\in[0,1]$. Initialized with data from a data-generating distribution $p_{\rm data}$, diffusion is typically parameterized such that the terminal distribution of the forward-time SDE is a tractable Gaussian distribution $\gamma$. This SDE shares marginal densities $p_{t}$ with the probability flow- (PF-)ODE:

$$ $dx_{t}=\left[\mu(x_{t},t)-\frac{1}{2}\beta(t)^{2}\nabla\log p_{t}(x_{t})\right ]dt.$ (1) $$

Score-based generative models are a class of DM which approximate $\nabla\log p_{t}(x_{t})$ with a neural network score model. Given such a model, (1) can be solved in reverse time using numerical ODE integration techniques [, ].

#### Deterministic diffusion solvers.

In contrast to stochastic DM samplers based on Markov chains [] and SDEs [], deterministic DM solvers primarily focus on simulating the PF-ODE (1). [] presented DDIM, an implicit modeling technique yielding a deterministic mapping between noise and data samples. Subsequent works considered alternate, higher-order solvers for the PF-ODE [], yielding high-quality samples in fewer function evaluations.

#### Flow models.

Continuous normalizing flows (CNFs) represent another class of ODE-based generative models, using neural networks to approximate the dynamics of a continuous mapping between noise and data []. Recent extensions focus on learning more direct trajectories [] and simulation-free training []. As with deterministic diffusion solvers corresponding to the PF-ODE, sampling via these methods requires numerical simulation of an ODE whose dynamics are defined by the neural network model.

#### Consistency models.

Efficient ODE simulation is of particular interest for efficient sampling from DMs [, ] and CNFs []. However, fast numerical ODE solvers still require tens of steps to produce high-fidelity samples [, ]. As a result, score model distillation techniques have arisen to yield fast, effective samplers from the PF-ODE. Consistency models (CMs) are a prominent class of distilled DMs that enable single- and few-step sampling []. CMs learn a mapping $f_{\theta}$ (parameterized by $\theta$) between a point $x_{t}$ along the PF-ODE trajectory to the initial state:

$$ $x_{0}=f_{\theta}(x_{t},t)\text{ for }t\in[0,1],$ (2) $$

where $x_{0}$ is a sample from $p_{\rm data}$. Therefore, single-step sampling can be achieved by sampling $x_{1}\sim\gamma$ and evaluating the CM at $x_{1}$. Multi-step sampling can be achieved by alternating denoising (via evaluation of the CM) and partial noising, trading off efficiency for fidelity.

## 3 Methodology

Assume that a pre-trained generative model is given, which provides a one-to-one mapping $\Phi$ from the noise space to the data space.
The data $x_{0}$ and noise $x_{1}$ both belong to $\mathbb{R}^{d}$, and
$x_{0}=\Phi(x_{1})$.
The observation is $y$, and the goal is to sample the data $x_{0}$ from the posterior distribution $p(x_{0}|y)$.
We derive the posterior sampling of the data vector $x_{0}$ via that of the noise vector $x_{1}$, making use of the mapping $\Phi$.

#### Likelihood and posterior.

We consider a general observation model where the conditional law $p(y|x_{0})$ is known and differentiable.
Define the negative log conditional likelihood as
$L_{y}(x_{0}):=-\log p(y|x_{0})$,
which is differentiable with respect to $x_{0}$ for fixed $y$. A typical case is the inverse problem setting: the forward model is

$$ $y=\mathcal{A}(x_{0})+n,$ (3) $$

where $\mathcal{A}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is the (possibly nonlinear) measurement operator, and $n$ is the additive noise. For fixed $y$, we aim to sample $x_{0}$ from
$p(x_{0}|y)={p(y|x_{0})p(x_{0})}/{p(y)}\propto p(y|x_{0})p(x_{0})$,
where $p(x_{0})$ is the true prior distribution of all data $x_{0}$, which we now denote as $p_{\rm data}$. We also call $p(x_{0}|y)$ the true posterior of $x_{0}$, donated as

$$ $p_{0,y}(x_{0}):=p(x_{0}|y)\propto p(y|x_{0})p_{\rm data}(x_{0}).$ (4) $$

#### Posterior approximated via generative model.

The true data prior $p_{\rm data}$ is nonlinear and complicated. Let $p_{\rm model}$ denote the prior distribution approximated by a pre-trained generative model $x_{0}=\Phi(x_{1})$, where $x_{1}\sim\gamma$.
A distribution from which samples are easily generated, such as the standard multi-variate Gaussian, is typically chosen for $\gamma$; we choose $\gamma=\mathcal{N}(0,I)$.
In other words,

$$ $p_{\rm data}\approx p_{\rm model}=\Phi_{\#}\gamma.$ (5) $$

Replacing $p_{\rm data}$ with $p_{\rm model}$ in (4) gives the model posterior of $x_{0}$, denoted $\tilde{p}_{0,y}$, which approximates the true posterior:

$$ $p_{0,y}(x_{0})\approx\tilde{p}_{0,y}(x_{0})\propto p(y|x_{0})\Phi_{\#}\gamma(x _{0}).$ (6) $$

Because $x_{0}=\Phi(x_{1})$, we have that
$\tilde{p}_{0,y}=\Phi_{\#}\tilde{p}_{1,y},$
where, by a change of variable from (6),

$$ $\tilde{p}_{1,y}(x_{1})\propto p(y|\Phi(x_{1}))\gamma(x_{1}).$ (7) $$

The distribution $\tilde{p}_{1,y}(x_{1})$ approximates the posterior distribution $p(x_{1}|y)$ in the noise space. When $p_{\rm data}=\Phi_{\#}\gamma$, we have $p_{0,y}=\tilde{p}_{0,y}$ and $p(\cdot|y)=\tilde{p}_{1,y}$.
When the generative model prior is inexact, the error in approximating the posterior can be bounded by that in approximating the data prior, see more in Section 4.

#### Posterior sampling by Langevin dynamics.

It is direct to sample the approximated posterior (7) in the noise space using Langevin dynamics.
Specifically, since we have
$\gamma(x_{1})\propto\exp(-\|x_{1}\|^{2}/2)$ and $\log p(y|\Phi(x_{1}))=-L_{y}(\Phi(x_{1}))$,
the following SDE of $x_{1}$ will have $\tilde{p}_{1,y}$ as its equilibrium distribution (proved in Lemma A.1):

$$ $dx_{1}=-(x_{1}+\nabla_{x_{1}}L_{y}(\Phi(x_{1})))dt+\sqrt{2}dW_{t}.$ (8) $$

The sampling in the noise space gives the sampling in the data space by the one-to-one mapping of the generative model, namely $x_{0}=\Phi(x_{1})$.

###### Example 3.1 (Inverse problem with Gaussian noise) .

For (3) with white noise, i.e., $n\sim\mathcal{N}(0,\sigma^{2}I)$, we have that, with a constant $c$ depending on $(\sigma,d)$,

$$ $L_{y}(x_{0})=-\log p(y|x_{0})=\frac{1}{2\sigma^{2}}\|y-\mathcal{A}(x_{0})\|_{2 }^{2}+c.$ $$

The noise-space SDE (8) can be written as

$$ $dx_{1}=-\left(x_{1}+\nabla_{x_{1}}\frac{\|\mathbf{y}-\mathcal{A}(x_{0})\|_{2}^ {2}}{2\sigma^{2}}\right)dt+\sqrt{2}dW_{t}.$ $$

Given $L_{y}(x_{0})$, standard techniques can be used to sample (overdamped) Langevin dynamics (8).
Evaluation of the gradient $\nabla_{x_{1}}L_{y}(x_{0})$ is the major computational cost, requiring differentiation through the model $\Phi$.
One technique to improve sampling efficiency is to employ a warm-start of the SDE integration by letting the minimization-only dynamics (using $\nabla_{x_{1}}L_{y}(x_{0})$) to converge to a minimum first, especially when the posterior concentrates around a particular point.
We postpone the algorithmic details to Section 5.

###### Example 3.1 (Inverse problem with Gaussian noise) .

## 4 Theory

In this section, we derive the theoretical guarantee of the model posterior $\tilde{p}_{0,y}$ in (6) to the true posterior $p_{0,y}$ in (4),
and also extend to the computed posterior $\tilde{p}_{0,y}^{S}$ by discrete-time SDE integration.
The analysis reveals a conditional number which indicates the intrinsic difficulty of the posterior sampling problem.
All proofs are in Appendix A.

### 4.1 Total Variation (TV) guarantee and condition number

Consider the approximation (5), that is, the pre-trained model generates a data prior distribution $\Phi_{\#}\gamma$ that approximates the true data prior $p_{\rm data}$.
We quantify the approximation in TV distance, namely

$$ $\operatorname{TV}(p_{\rm data},\Phi_{\#}\gamma)\leq\varepsilon.$ (9) $$

Generation guarantee in terms of TV bound has been derived in several flow-based generative model works, such as [, , ] on the PF-ODE of a trained score-based diffusion model [],
and [] on the JKO-type flow model [].
The following theorem proved in Appendix A shows that the TV distance between the model and true posteriors can be bounded proportional to that between the priors.

###### Theorem 4.1 (TV guarantee) .

Assuming (9), then
$\operatorname{TV}(p_{0,y},\tilde{p}_{0,y})\leq 2\kappa_{y}\varepsilon$,
where

$$ $\kappa_{y}:=\frac{\sup_{x_{0}}p(y|x_{0})}{\int p(y|x)p_{\rm data}(x)dx}.$ (10) $$

###### Remark 4.1 ( $\kappa_{y}$ as a condition number) .

The constant factor $\kappa_{y}$ is determined by the true data prior $p_{\rm data}$ and the conditional likelihood $p(y|x_{0})$ of the observation, and is independent of the flow model and the posterior sampling method.
Thus $\kappa_{y}$ quantifies an intrinsic “difficulty” of the posterior sampling, which can be viewed as a condition number of the problem.

###### Example 4.1 (Well-conditioned problem) .

Suppose $p(y|x_{0})\leq c_{1}$ for any $x_{0}$,
and on a domain $\Omega_{y}$ of the data space,

$$ $P_{\rm data}(\Omega_{y})\geq\alpha>0,\quad\text{ and }\quad p(y|x_{0})\geq c_{ 0}>0,\,\forall x_{0}\in\Omega_{y},$ $$

then we have $\int p(y|x)p_{\rm data}(x)dx\geq\int_{\Omega_{y}}p(y|x)p_{\rm data}(x)dx\geq
\alpha c_{0}$, and then

$$ $\kappa_{y}\leq\frac{1}{\alpha}\frac{c_{1}}{c_{0}}.$ $$

This shows that if the observation $y$ can be induced from some cohort of $x_{0}$ and this cohort is well-sampled by the data prior $p_{\rm data}$ (the concentration of $p_{\rm data}$ on this cohort is lower bounded by $\alpha$),
plus that the most likely $x_{0}$ is not too peaked compared to the likelihood of any other $x_{0}$ within this cohort (the ratio is upper bounded by $c_{1}/c_{0}$), then the posterior sampling is well-conditioned.

###### Example 4.2 (Ill-conditioned problem) .

Suppose $p(y|x_{0})$ is peaked at one data value $x_{0}^{\prime}$ and almost zero at other places, and this $x_{0}^{\prime}$ lies on the tail of the data prior density $p_{\rm data}$.
This means that the integral $\int p(y|x_{0})p_{\rm data}(x_{0})dx_{0}$ has all the contribution on a nearby neighborhood of $x_{0}^{\prime}$ on which $p_{\rm data}$ is small, resulting in a small value on the denominator of (10). Meanwhile, the value of $p(y|x_{0}^{\prime})$ is large.
In this case, $\kappa_{y}$ will take a large value, indicating an intrinsic difficulty of the problem.
Intuitively, the desired data value $x_{0}^{\prime}$ for this observation $y$ is barely represented within the (unconditional) data distribution $p_{\rm data}$, while the generative model can only learn from $p_{\rm data}$.
Since the pre-trained unconditional generative model does not have enough knowledge of such $x_{0}^{\prime}$, it is hard for the conditional generative model (based on the unconditional model) to find such a data value.

###### Theorem 4.1 (TV guarantee) .

###### Remark 4.1 ( κ y subscript 𝜅 𝑦 \kappa_{y} italic_κ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT as a condition number) .

###### Example 4.1 (Well-conditioned problem) .

###### Example 4.2 (Ill-conditioned problem) .

### 4.2 TV guarantee of the sampled posterior

Theorem 4.1 captures the approximation error of $\tilde{p}_{0,y}$ to the true posterior, where $\tilde{p}_{0,y}$ is the distribution of data $x_{0}$ when the noise $x_{1}$ in noise space achieves the equilibrium $\tilde{p}_{1,y}$ of the SDE (8).
In practice, we use a numerical solver to sample the SDE in discrete time.
The convergence of discrete-time SDE samplers to its equilibrium distribution has been established under various settings in the literature.
Here, we assume that the discrete-time algorithm to sample the Langevin dynamics of $x_{1}$ outputs $x_{1}\sim\tilde{p}_{1,y}^{S}$, which may differ from but is close to the equilibrium $\tilde{p}_{1,y}$.
Specifically, suppose $\operatorname{TV}(\tilde{p}_{1,y},\tilde{p}_{1,y}^{S})$ is bounded by some $\varepsilon_{S}$.

###### Lemma 4.2 (Sampling error) .

If $\operatorname{TV}(\tilde{p}_{1,y},\tilde{p}_{1,y}^{S})\leq\varepsilon_{S}$, then
$\operatorname{TV}(\tilde{p}_{0,y},\tilde{p}_{0,y}^{S})\leq\varepsilon_{S}$.

The lemma is by Data Processing Inequality, and together with Theorem 4.1 it directly leads to the following corollary on the TV guarantee of the sampled posterior.

###### Corollary 4.3 (TV of sampled posterior) .

Assuming (9)
and $\operatorname{TV}(\tilde{p}_{1,y},\tilde{p}_{1,y}^{S})\leq\varepsilon_{S}$,
then

$$ $\operatorname{TV}(p_{0,y},\tilde{p}_{0,y}^{S})\leq 2\kappa_{y}\varepsilon+ \varepsilon_{S}.$ $$

###### Lemma 4.2 (Sampling error) .

###### Corollary 4.3 (TV of sampled posterior) .

## 5 Algorithm

#### Numerical integration of the Langevin dynamics.

To numerically integrate the noise-space SDE (8), one can use standard SDE solvers. We adopt the Euler-Maruyama (EM) scheme. Let $\tau>0$ be the time step, and denote the discrete sequence of $x_{1}$ as $z^{i}$, $i=0,1,\cdots$. The EM scheme gives, with $\xi^{i}\sim\mathcal{N}(0,I)$,

$$ $z^{i+1}=(1-\tau)z^{i}-\tau{g}^{i}+\sqrt{2\tau}\xi^{i},\quad{g}^{i}:=\nabla_{x_ {1}}L_{y}(x_{0})|_{x_{1}=z^{i}}.$ (11) $$

See Algorithm 1 for an outline of our approach using EM. However, any general numerical scheme for solving SDEs can be applied; see Table A.4 in Appendix C for a comparison between our method using EM discretization and exponential integrator (EI) []. An initial value of $z^{0}$ in the noise space is also required. We adopt a warm-start procedure to initialize sampling; additional details are provided below.

Figure: Algorithm 1 Posterior Sampling in Noise Space

#### Computation of ∇ x 1 ⁢ L y ⁢ ( x 0 ) ∇ subscript 𝑥 1 subscript 𝐿 𝑦 subscript 𝑥 0 \nabla x_{1}L_{y}(x_{0}) ∇ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_L start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) .

The computation of the loss gradient depends on the type of generative model representing the mapping $\Phi$. For instance, if $\Phi$ is computed by solving an ODE driven by a normalizing flow, then its gradient can be computed using the adjoint sensitivity method []. If $\Phi$ is a DM or CM sampling scheme, one can backpropagate through the nested function calls to the generative model. Since we use one- or few-step CM sampling to represent $\Phi$ in the experiments, we take the latter approach to compute $\nabla x_{1}L_{y}(x_{0})$.

#### Choice of initial value and warm-start.

A natural choice for the initial noise space value $z^{0}$ can be a generic sample $z^{0}\sim\gamma$ (the noise space prior). However, while this will correspond to a high-likelihood sample according to the data prior (given a well-trained generative model), it may be far from the data posterior. As such, one may warm-start the sampler by optimizing $L_{y}(x_{0})$ with respect to $x_{1}$ using standard optimization techniques such as gradient descent or Adam. In all experiments, we warm start sampling using $K$ steps of Adam optimization, initializing EM sampling with $z^{0}$ being the optimization output. See Appendix B.1 for further detail.

#### Computational requirements.

The main computational burden is with respect to the computation of the loss gradient $\nabla_{x_{1}}L_{y}(x_{0})$, which requires differentiating through the mapping $\Phi$. This can be alleviated by choosing a $\Phi$ which consists of a small number of function evaluations (NFEs). Additional computational burden is due to burn-in/warm start to yield $z^{0}$, the initial value of EM simulation. Therefore, the total NFEs to simulate $N$ steps of EM (i.e., to yield $N$ samples) is $\eta\cdot(K+N)$, where $\eta$ is the NFEs required to evaluate $\Phi$. However, this burden is amortized over EM sampling, as progressive EM simulation yields increasingly fewer overall NFEs per sample, which asymptotically approaches $\eta$ (the NFEs required to compute $\Phi$). Therefore, we represent $\Phi$ using CM sampling, which can be accomplished for $\eta=1$ or 2. While multi-step ($\eta>1$) CM sampling is typically stochastic [], we fix the noise in each step to result in a deterministic mapping. See Appendix B.1 for details.

#### Role of step size τ 𝜏 \tau italic_τ .

The step size of EM, $\tau$, controls the time scales over which the Langevin dynamics are simulated with respect to the number of EM steps. Larger $\tau$ results in more rapid exploration of the posterior, potentially leading to more diverse samples over shorter timescales. However, $\tau$ must also be kept small enough to ensure the stability of EM sampling. Therefore, this hyper-parameter provides a degree of control over the diversity of samples provided by the proposed algorithm. Choosing large $\tau$ while maintaining stability can yield diverse samples, potentially revealing particularly uncertain semantic features within the posterior.

**Table 1: Quantitative comparison of linear image restoration tasks on LSUN-Bedroom (256 x 256) (top table) and ImageNet (64 x 64) (bottom table).**
| Method | 8x Super-resolution | Gaussian Deblur | 10$\%$ Inpainting |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | FID $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS$\downarrow$ | FID $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS$\downarrow$ | FID $\downarrow$ |  |
| DPS-DM | 20.4^∗ | 0.538^∗ | 0.470^∗ | 67.7^∗ | 22.1 | 0.589 | 0.407 | 65.3 | 22.4 | 0.634 | 0.417 | 67.7 |
| MPGD-DM | 19.2 | 0.338 | 0.689 | 288 | 23.6^∗ | 0.579 | 0.438 | 85.0 | 15.4 | 0.176 | 0.667 | 221 |
| LGD-DM | 20.1 | 0.529 | 0.483 | 69.3 | 22.2 | 0.590^∗ | 0.371^∗ | 60.1^∗ | 24.7^∗ | 0.742^∗ | 0.289^∗ | 47.3^∗ |
| DPS-CM | 10.7 | 0.077 | 0.758 | 307 | 11.2 | 0.092 | 0.735 | 279 | 19.9 | 0.454 | 0.517 | 128 |
| LGD-CM | 10.5 | 0.072 | 0.764 | 316 | 11.1 | 0.092 | 0.737 | 283 | 19.9 | 0.475 | 0.514 | 134 |
| CMEdit | N/A | N/A | 18.0 | 0.523 | 0.548 | 167 |  |  |  |  |  |  |
| Ours(1-step) | 20.4 | 0.535 | 0.418 | 71.1 | 22.4 | 0.598 | 0.368 | 70.6 | 23.8 | 0.682 | 0.358 | 72.9 |
| Ours(2-step) | 20.5 | 0.534 | 0.433 | 72.2 | 21.3 | 0.554 | 0.421 | 69.2 | 22.2 | 0.611 | 0.419 | 75.6 |

Figure: Figure 3: Image reconstructions for the linear and nonlinear tasks on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/lsun_fidelity.png

Figure: Figure 4: Image reconstructions for the linear tasks on ImageNet (64 x 64).
Refer to caption: extracted/5896580/imnet_fidelity.png

## 6 Experiments

#### Baselines.

We categorize the baselines into two groups. (1) DM-based methods: Diffusion Posterior Sampling (DPS) [], Loss-Guided Diffusion (LGD) [], and Manifold-Preserving Guided Diffusion (MPGD) []. These methods employ stronger priors compared to our approach, making them inherently stronger baselines and rendering the comparison across different backbones unfair.

To ensure a fairer comparison, we adopt a second set of baselines: (2) CM-based methods, where each DM-based method is adapted to use a consistency model (CM) backbone. Additionally, we include CMEdit, the modified CM sampler from [], for linear tasks. All DM baselines use the same EDM model from [], and all CM baselines use the corresponding LPIPS-distilled CM. Details and hyper-parameters for each baseline are outlined in Appendix B.2.

#### Datasets.

We include experiments on LSUN-Bedroom (256 x 256) [] and ImageNet (64 x 64) [], using 100 validation images for each dataset. All experiments are conducted using the pre-trained CMs from [], which were distilled using the LPIPS objective from pre-trained EDM models []. See Appendix B.1 for additional details regarding our method and hyper-parameters. We consider the following linear forward operators for inverse problem tasks: (i) for random mask inpainting, some percentage of the pixels are masked uniformly at random; (ii) for super-resolution, adaptive average pooling is applied; and (iii) for Gaussian deblurring, we use a kernel of $61\times 61$ pixels with standard deviation 3.0. We also consider nonlinear tasks: (i) nonlinear deblurring using a neural network forward model []; (ii) for phase retrieval, the magnitude of the Fourier coefficients is computed; and (iii) for high dynamic range (HDR) reconstruction, pixel values are multiplied by 2 and again truncated to [-1,1]. All experiments apply Gaussian noise with standard deviation $\sigma=0.1$ in the measurement space (except for phase retrieval experiments, which use $\sigma=0.05$). See Appendix B.3 for detailed descriptions of the forward operators. Additional experimental results can be found in Appendices C and D.

#### Metrics.

To assess reconstruction fidelity, we compare samples from each method using the Peak Signal-to-Noise Ratio (PSNR), Structural Similarity Index Metric (SSIM), Learned Perceptual Image Patch Similarity (LPIPS), and Fréchet Inception Distance (FID). To assess the diversity of samples, we consider the following metrics: (i) Diversity Score (DS), which is the ratio between the inter- and intra-cluster distances using 6 nearest neighbors clusters of ResNet-50 features, and (ii) Average CLIP Cosine Similarity (CS), which is the average cosine similarity between CLIP embeddings all sample pairs for a given image.

### 6.1 Image restoration results

#### Linear inverse problems.

We quantitatively compare the performance of the proposed approach to the baselines for point-estimate image restoration under linear forward models, where 10 samples are provided by each method for 100 images in the validation datasets. LSUN-Bedroom (256 x 256) results are reported in the top section of Table 1 and our approach is compared to the highest-fidelity baselines on ImageNet (64 x 64) in the bottom section of Table 1. Visual comparisons of point estimates are also visualized in the top three rows of Figure 3 (for LSUN) and in Figure 4 (for ImageNet). Compared with CM baselines, the proposed approach exhibits superior performance in producing high-fidelity candidate solutions to linear inverse problems. This corresponds to improved visual quality, as other CM approaches produce artifacts and poor reconstructions of the ground truth. The proposed method is also competitive against DM baselines, yielding samples of comparable quality both qualitatively and quantitatively.

**Table 2: Quantitative comparison of nonlinear image restoration tasks on LSUN-Bedroom (256 x 256).**
| Method | Nonlinear Deblur | Phase Retrieval | HDR Reconstruction |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | FID $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS$\downarrow$ | FID $\downarrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS$\downarrow$ | FID $\downarrow$ |  |
| DPS-DM | 21.6 | 0.586 | 0.413 | 75.7^∗ | 10.7 | 0.302 | 0.697^∗ | 90.1 | 21.7^∗ | 0.659^∗ | 0.396^∗ | 69.6^∗ |
| MPGD-DM | 17.0 | 0.194 | 0.683 | 259 | 9.96 | 0.271 | 0.728 | 118 | 20.5 | 0.586 | 0.408 | 73.2 |
| LGD-DM | 22.3^∗ | 0.632^∗ | 0.408^∗ | 106 | 10.8^∗ | 0.351^∗ | 0.709 | 82.0^∗ | 12.4 | 0.459 | 0.560 | 172 |
| DPS-CM | 17.7 | 0.303 | 0.574 | 137 | 10.1 | 0.197 | 0.726 | 195 | 13.5 | 0.405 | 0.597 | 173 |
| MPGD-CM | 13.1 | 0.100 | 0.762 | 306 | 9.39 | 0.111 | 0.786 | 312 | 11.7 | 0.296 | 0.638 | 223 |
| LGD-CM | 21.3 | 0.519 | 0.482 | 163 | 9.36 | 0.113 | 0.767 | 186 | 11.2 | 0.397 | 0.621 | 245 |
| Ours(1-step) | 20.3 | 0.566 | 0.440 | 76.7 | 10.3 | 0.315 | 0.709 | 82.9 | 19.6 | 0.599 | 0.436 | 88.0 |
| Ours(2-step) | 18.7 | 0.501 | 0.492 | 73.3 | 10.2 | 0.309 | 0.708 | 81.4 | 16.6 | 0.481 | 0.532 | 101 |
| Bold denotes the best CM method, underline denotes the second best CM method, and ^∗ denotes the best DM method. |  |  |  |  |  |  |  |  |  |  |  |  |

Figure: Figure 5: Posterior samples for the inpainitng (10%) (top three rows) and nonlinear deblur (bottom three rows) tasks on LSUN-Bedroom (256 x 256). Green boxes highlight low-uncertainty features and red boxes highlight highly uncertain features.
Refer to caption: extracted/5896580/lsun_diversity.png

#### Nonlinear inverse problems.

Quantitative comparisons for nonlinear tasks on 100 images from LSUN-Bedroom are displayed in Table 2, where metrics are again computed using 10 samples per image from each method. The proposed method is highly competitive against CM-backbone baselines in all tasks. Moreover, the performance is comparable to that of the DM-backbone baselines. Example reconstructions for each method are visualized in the bottom three rows of Figure 3. Other CM-based methods and MPGD-DM seemingly fail to remove the degradation and noise applied by the forward process, while the proposed method yields samples of visual quality comparable to that of DM baselines. Reconstructions generated using the proposed approach lack the artifacts of CM-backbone baselines while also capturing the fine details present in DM reconstructions. In particular, in the highly degraded and ill-posed phase retrieval task, our method yields samples that are markedly consistent with the ground truth, as PSNR and SSIM values are comparable to those of DM baselines.

### 6.2 Diversity of posterior samples

To assess the capacity of the proposed approach to generate diverse samples from the posterior, we conduct additional experiments comparing our method to the strongest baselines: DPS and LGD with a DM backbone. For each of the six (linear and nonlinear) tasks, we generate 25 samples for 100 images from the validation partition of LSUN-Bedroom (256 x 256) via each method. A quantitative comparison of the diversity of the samples from each method is shown in Table 3. Generally, the proposed approach provides competitive to superior performance in diversity metrics compared to DM baselines. Furthermore, visualizing a subset of the posterior samples in the inpainting (top three rows) and nonlinear deblurring (bottom three rows) tasks in Figure 5, one can observe that samples from our method have more clear visual diversity. High-level features of the scene, such as overall lighting or shading, are more variable across our posterior samples. Moreover, our method can identify certain and uncertain semantic features in the candidate reconstructions, as particular features such as windows and lamps have dramatic qualitative variation across the posterior samples from our approach.

**Table 3: Quantitative comparison of diversity metrics on linear and non-linear image restoration tasks on LSUN-Bedroom (256 x 256).**
| Method | SR(8x) | Gaussian Deblur | 10% Inpainting | Nonlinear Deblur | Phase Retrieval | HDR Reconstruction |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| DS $\uparrow$ | CS $\downarrow$ | DS $\uparrow$ | CS $\downarrow$ | DS $\uparrow$ | CS $\downarrow$ | DS $\uparrow$ | CS $\downarrow$ | DS $\uparrow$ | CS $\downarrow$ | DS $\uparrow$ | CS $\downarrow$ |  |
| DPS-DM | 2.14 | 0.843 | 2.10 | 0.938 | 2.33 | 0.876 | 2.22 | 0.924 | 2.42 | 0.809 | 2.25 | 0.873 |
| LGD-DM | 2.35 | 0.881 | 2.19 | 0.925 | 2.28 | 0.872 | 2.11 | 0.923 | 2.36 | 0.815 | 3.14 | 0.914 |
| Ours(1-step) | 3.01 | 0.879 | 3.26 | 0.997 | 3.15 | 0.869 | 2.80 | 0.912 | 3.08 | 0.914 | 3.09 | 0.927 |
| Ours(2-step) | 2.67 | 0.919 | 2.62 | 0.866 | 2.48 | 0.864 | 2.69 | 0.885 | 2.89 | 0.862 | 3.23 | 0.904 |
| Bold denotes the best method, underline denotes the second best method. |  |  |  |  |  |  |  |  |  |  |  |  |

## 7 Related works

#### Posterior sampling with generative models.

Diffusion-based inverse problem solvers consist of task-specific frameworks [, , ], optimized approaches [, , ], and training-free techniques leveraging pre-trained diffusion priors [, , , , , , , , , ]. Early training-free methods for solving inverse problems utilize measurement-space projection [, ], while others addressed noisy problems via consistency in the spectral domain [, , ] or using manifold constraints [, ]. Recent works consider general noisy and nonlinear inverse problems using an approximation of the measurement likelihood in each generation step [, , ]. An emerging area of interest focuses on developing diffusion posterior sampling techniques with provable guarantees [, ]. For instance, [] develop an alternating measurement projection/guided diffusion approach for which they provide asymptotic convergence guarantees, while [] utilize tilted transport in linear inverse problems which provably samples the posterior under certain conditions. Diffusion-base posterior sampling works can also be adapted to flow-based models, e.g., [] adapt $\Pi$GDM [] to CNFs. These existing works modify the sampling trajectory of generative priors, requiring repeated simulation of the entire sampling process to produce multiple posterior samples, hindering scalability to many samples. The proposed sampling in the noise space of one- or few-step mappings enables the efficient generation of many posterior samples.

#### Guided generation via noise space iteration.

For generative models that provide deterministic mappings between a latent noise space and data, such as GANs [], flows [], and CMs [], optimization of noise can guide generation towards conditional information [, , , , , ]. In the GAN literature, this is primarily addressed using text-to-image guided synthesis [, ] or task-specific objectives []. This type of approach has also been used to solve inverse problems using flow-based models [, ]; for instance, D-Flow [] optimizes with respect to the noise input to CNFs. Our method also iterates in the noise space, simulating Langevin dynamics for posterior sampling instead of optimizing to yield a point estimate. Computing gradients through CNFs is expensive [], requiring at least tens of function evaluation per ODE solution [, ]. The use of CMs in our approach facilitates computation of the gradient in as few as one call to the neural network, enabling the progressive accumulation of posterior samples during Langevin dynamics simulation.

## 8 Discussion

We have outlined a novel approach for posterior sampling via Langevin dynamics in the noise space of a generative model. Using a CM mapping from noise to data, our posterior sampling provides solutions to general noisy image inverse problems, demonstrating superior reconstruction fidelity to other CM methods and competitiveness with diffusion baselines. A primary limitation of our approach is the low visual quality in some posterior samples. Fidelity drawbacks can be attributed to a relatively poor approximation of the prior by CMs. Future work will focus on improving fidelity of diverse samples, perhaps by using more accurate prior models and adaptive simulation of the SDE. Regardless, our method produces highly diverse samples, representing meaningful semantic uncertainty of data features within the posterior.

## Acknowledgments

MR, YX, and XC were partially supported by National Science Foundation (NSF) DMS-2134037.
YX was also partially supported by
NSF DMS-2220495.
XC was also partially supported by
NSF DMS-2237842
and Simons Foundation (grant ID: 814643).
JL was partially supported by NSF DMS-2309378.

## Appendix A Proofs

###### Lemma A.1 .

The equilibrium distribution of SDE (8) is $\tilde{p}_{1,y}$.

###### Proof of Lemma A.1 .

Under generic condition, the Langevin dynamics

$$ $dX_{t}=-\nabla U(X_{t})dt+\sqrt{2}dW_{t}$ $$

have the equilibrium $\rho_{\infty}\propto e^{-U}$.
For $\tilde{p}_{1,y}$ in (7) to be the equilibrium, it suffices to verify that

$$ $\nabla\log\tilde{p}_{1,y}=-(x_{1}+\nabla_{x_{1}}L_{y}(\Phi(x_{1}))).$ $$

This follows by that $\log\gamma(x_{1})=-\|x_{1}\|^{2}/2+c$ and
$\log p(y|\Phi(x_{1}))=-L_{y}(\Phi(x_{1}))$.
∎

###### Proof of Theorem 4.1 .

By (4) and (6), we have

$$ $p_{0,y}(x_{0})=\frac{1}{Z_{y}}p(y|x_{0})p_{\rm data}(x_{0}),\quad\tilde{p}_{0, y}(x_{0})=\frac{1}{\tilde{Z}_{y}}p(y|x_{0})\Phi_{\#}\gamma(x_{0}),$ $$

where

$$ $Z_{y}:=\int p(y|x_{0})p_{\rm data}(x_{0})dx_{0},\quad\tilde{Z}_{y}:=\int p(y|x _{0})\Phi_{\#}\gamma(x_{0})dx_{0}.$ $$

Then, we have

$$ $\displaystyle 2\operatorname{TV}(p_{0,y},\tilde{p}_{0,y})$ $\displaystyle=\int|p_{0,y}(x_{0})-\tilde{p}_{0,y}(x_{0})|dx_{0}$ $\displaystyle\leq\int\frac{1}{Z_{y}}p(y|x_{0})\left|p_{\rm data}(x_{0})-\Phi_{ \#}\gamma(x_{0})\right|dx_{0}+\left|\frac{\tilde{Z}_{y}-Z_{y}}{Z_{y}}\right|.$ (A.1) $$

By definition of $\kappa_{y}$ in (10), we have
$\frac{1}{Z_{y}}p(y|x_{0})\leq\kappa_{y}$, $\forall x_{0}$,
and thus

$$ $\int\frac{1}{Z_{y}}p(y|x_{0})\left|p_{\rm data}(x_{0})-\Phi_{\#}\gamma(x_{0}) \right|dx_{0}\leq\kappa_{y}\int\left|p_{\rm data}(x_{0})-\Phi_{\#}\gamma(x_{0} )\right|dx_{0}.$ $$

Meanwhile,
$\tilde{Z}_{y}-Z_{y}=\int p(y|x_{0})(\Phi_{\#}\gamma(x_{0})-p_{\rm data}(x_{0})
)dx_{0},$
and then

$$ $\displaystyle\frac{|\tilde{Z}_{y}-Z_{y}|}{Z_{y}}$ $\displaystyle\leq\int\frac{1}{Z_{y}}p(y|x_{0})|\Phi_{\#}\gamma(x_{0})-p_{\rm data }(x_{0})|dx_{0}$ $\displaystyle\leq\int\kappa_{y}|\Phi_{\#}\gamma(x_{0})-p_{\rm data}(x_{0})|dx_ {0}.$ $$

Putting back to (A.1), we have

$$ $2\operatorname{TV}(p_{0,y},\tilde{p}_{0,y})\leq 2\kappa_{y}\int\left|p_{\rm data }(x_{0})-\Phi_{\#}\gamma(x_{0})\right|dx_{0}=4\kappa_{y}\operatorname{TV}(p_{ \rm data},\Phi_{\#}\gamma),$ $$

which proves the theorem under (9).
∎

###### Proof of Lemma 4.2 .

By that $\tilde{p}_{0,y}=\Phi_{\#}\tilde{p}_{1,y}$,
$\tilde{p}_{0,y}^{S}=\Phi_{\#}\tilde{p}_{1,y}^{S}$, and Data Processing Inequality.
∎

###### Proof of Corollary 4.3 .

By Theorem 4.1, Lemma 4.2, and triangle inequality since TV is half of the $L^{1}$ norm between two densities.
∎

###### Lemma A.1 .

###### Proof of Lemma A.1 .

###### Proof of Theorem 4.1 .

###### Proof of Lemma 4.2 .

###### Proof of Corollary 4.3 .

## Appendix B Experimental details

### B.1 Details of the proposed approach

#### Consistency model generative process.

To represent the map $\Phi$ from noise space to data space, we utilize the pre-trained CMs of [] with a 1- or 2-step sampler. For the 2-step sampler, we use standard multistep consistency sampling (Algorithm 1, []), i.e.,

$$ $x_{0}=f_{\theta}\left(f_{\theta}(x_{T},T)+\sqrt{t^{2}-\epsilon^{2}}z,t\right),$ $$

where $f_{\theta}$ is the pre-trained CM, $x_{T}\leftarrow x_{1}$, $T=80$, $\epsilon=2\times 10^{-3}$ is a small noise offset , and $t$ is an intermediate "time step" along the PF-ODE trajectory (the "halfway" point). In [], $z$ is sampled from the standard Gaussian for each call to $\Phi$. In this work, we sample $z$ once and fix it for all future calls to $\Phi$, which we observe to empirically improve performance.

#### Warm-start initialization and sampling.

The posterior sampling process begins with a warm-start initialization consisting of $K$ steps of Adam optimization with learning rate, $\beta_{1}$, and $\beta_{2}$ for each experiment outlined in Tables A.1, A.2, and A.3. This is followed by $N$ steps of Langevin dynamics simulation (via EM discretization in the main-text experiments) using step size $\tau$. The NFEs per sample can be computed as $\eta(K+N)/N$, where $\eta$ is the number of steps used for CM generation. All experiments are implemented in PyTorch and are run on a system with NVIDIA A100 GPUs.

See below for a pseudo-code implementation of one iteration of our sampling procedure:

**Table A.1: Hyper-parameters for linear and nonlinear image restoration tasks on LSUN-Bedroom (256 x 256).**
| Method | 8x Super-resolution | Gaussian Deblur | 10$\%$ Inpainting | Nonlinear Deblur | Phase Retrieval | HDR Reconstruction |
| --- | --- | --- | --- | --- | --- | --- |
| DPS-DM | $\zeta=25,N=100$ | $\zeta=7,N=100$ | $\zeta=25,N=100$ | $\zeta=15,N=100$ | $\zeta=10,N=100$ | $\zeta=5,N=100$ |
| MPGD-DM | $\zeta=25,N=100$ | $\zeta=15,N=100$ | $\zeta=25,N=100$ | $\zeta=7,N=100$ | $\zeta=1,N=100$ | $\zeta=5,N=100$ |
| LGD-DM | $\zeta=25,M=1,N=100$ | $\zeta=25,M=10,N=100$ | $\zeta=7,M=25,N=100$ | $\zeta=9,M=10,N=100$ | $\zeta=1,M=10,N=100$ | $\zeta=30,M=10,N=100$ |
| DPS-CM | $\zeta=25,N=100$ | $\zeta=7,N=100$ | $\zeta=25,N=100$ | $\zeta=8,N=100$ | $\zeta=9,N=100$ | $\zeta=4,N=100$ |
| MPGD-CM | N/A | N/A | N/A | $\zeta=15,N=100$ | $\zeta=3,N=100$ | $\zeta=30,N=100$ |
| LGD-CM | $\zeta=25,M=1,N=100$ | $\zeta=7,M=1,N=100$ | $\zeta=5,M=1,N=100$ | $\zeta=15,M=10,N=100$ | $\zeta=0.5,M=10,N=100$ | $\zeta=15,M=10,N=100$ |
| Ours(1-step) | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=200,{\rm lr}=1\times 10^{-3}$ | $K=800,{\rm lr}=5\times 10^{-3}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=1\times 10^{-5}$ | EM: $N=10,\tau=1\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-5}$ | EM: $N=10,\tau=5\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-6}$ |  |
| Ours(2-step) | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=500,{\rm lr}=5\times 10^{-3}$ | Adam: $K=500,{\rm lr}=1\times 10^{-3}$ | Adam: $K=500,{\rm lr}=5\times 10^{-3}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=1\times 10^{-5}$ | EM: $N=10,\tau=1\times 10^{-7}$ | EM: $N=10,\tau=1\times 10^{-5}$ | EM: $N=10,\tau=5\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-6}$ |  |

**Table A.2: Hyper-parameters for linear image restoration tasks on ImageNet (64 x 64).**
| Method | 4x Super-resolution | Gaussian Deblur | 20$\%$ Inpainting |
| --- | --- | --- | --- |
| DPS-DM | $\zeta=20,N=100$ | $\zeta=15,N=100$ | $\zeta=30,N=100$ |
| LGD-DM | $\zeta=3,M=10,N=100$ | $\zeta=1,M=10,N=100$ | $\zeta=5,M=10,N=100$ |
| DPS-CM | $\zeta=30,N=100$ | $\zeta=30,N=100$ | $\zeta=25,N=100$ |
| LGD-CM | $\zeta=3,M=10,N=100$ | $\zeta=7,M=10,N=100$ | $\zeta=6,M=10,N=100$ |
| Ours(1-step) | Adam: $K=800,{\rm lr}=1\times 10^{-2}$ | Adam: $K=800,{\rm lr}=1\times 10^{-2}$ | $K=800,{\rm lr}=1\times 10^{-2}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=5\times 10^{-4}$ | EM: $N=10,\tau=3\times 10^{-5}$ | EM: $N=10,\tau=1\times 10^{-4}$ |  |
| Ours(2-step) | Adam: $K=500,{\rm lr}=5\times 10^{-2}$ | Adam: $K=500,{\rm lr}=5\times 10^{-2}$ | Adam: $K=500,{\rm lr}=5\times 10^{-2}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=1\times 10^{-4}$ | EM: $N=10,\tau=3\times 10^{-5}$ | EM: $N=10,\tau=1\times 10^{-4}$ |  |

**Table A.3: Hyper-parameters for linear and nonlinear diversity experiments on LSUN-Bedroom (256 x 256).**
| Method | 8x Super-resolution | Gaussian Deblur | 10$\%$ Inpainting | Nonlinear Deblur | Phase Retrieval | HDR Reconstruction |
| --- | --- | --- | --- | --- | --- | --- |
| DPS-DM | $\zeta=7,N=100$ | $\zeta=7,N=100$ | $\zeta=7,N=100$ | $\zeta=5,N=100$ | $\zeta=5,N=100$ | $\zeta=1,N=100$ |
| LGD-DM | $\zeta=15,M=1,N=100$ | $\zeta=5,M=1,N=100$ | $\zeta=15,M=1,N=100$ | $\zeta=4,M=10,N=100$ | $\zeta=0.5,M=10,N=100$ | $\zeta=10,M=10,N=100$ |
| Ours(1-step) | Adam: $K=400,{\rm lr}=5\times 10^{-3}$ | Adam: $K=600,{\rm lr}=5\times 10^{-3}$ | $K=600,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=200,{\rm lr}=1\times 10^{-3}$ | $K=800,{\rm lr}=5\times 10^{-3}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=4\times 10^{-4}$ | EM: $N=10,\tau=1\times 10^{-6}$ | EM: $N=10,\tau=1\times 10^{-4}$ | EM: $N=25,\tau=7.5\times 10^{-6}$ | EM: $N=25,\tau=3\times 10^{-6}$ | EM: $N=25,\tau=3\times 10^{-6}$ |  |
| Ours(2-step) | Adam: $K=600,{\rm lr}=5\times 10^{-3}$ | Adam: $K=600,{\rm lr}=5\times 10^{-3}$ | Adam: $K=800,{\rm lr}=5\times 10^{-3}$ | Adam: $K=500,{\rm lr}=5\times 10^{-3}$ | Adam: $K=500,{\rm lr}=1\times 10^{-3}$ | Adam: $K=500,{\rm lr}=5\times 10^{-3}$ |
| $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ | $\beta_{1}=0.9,\beta_{2}=0.999$ |  |
| EM: $N=10,\tau=4\times 10^{-4}$ | EM: $N=10,\tau=1\times 10^{-5}$ | EM: $N=10,\tau=1\times 10^{-4}$ | EM: $N=25,\tau=7.5\times 10^{-6}$ | EM: $N=25,\tau=3\times 10^{-6}$ | EM: $N=25,\tau=3\times 10^{-6}$ |  |

### B.2 Details of the baselines

The baseline methods conduct $t=1,\dots,N$ Euler steps for sampling. All methods require a denoiser to provide $x_{0}\approx\hat{x}_{0}(x_{t})$ at each sampling step $t$, which is achieved using either a pre-trained EDM [] or CM [], both obtained from [] for each dataset.

#### Diffusion Posterior Sampling (DPS).

DPS [] utilizes the denoiser corresponding to a pre-trained DM to approximate the measurement likelihood gradient at each step of DM sampling. At each state $x_{t}$ along the diffusion sampling trajectory, a score-base diffusion model can provide a predicted $\hat{x}_{0}(x_{t})$, which can be used to compute $\nabla_{x_{t}}p(y|\hat{x}_{0})$ via differentiation through the score-based model. In DPS, each step of diffusion sampling is adjusted by this gradient with weight $\zeta$, i.e., $x_{t-1}\leftarrow x_{t-1}-\zeta\nabla_{x_{t}}p(y|\hat{x}_{0})$.

#### Manifold Preserving Guided Diffusion (MPGD).

MPGD [] computes the gradient of the measurement likelihood in the denoised space rather than with respect to $x_{t}$ at each step, taking a gradient step in $\hat{x}_{0}$ before updating the diffusion iterate. That is, MPGD conducts the update $\hat{x}_{0}\leftarrow\hat{x}_{0}(x_{t})-\zeta\nabla_{\hat{x}_{0}}p(y|\hat{x}_{
0}(x_{t}))$, which can then be use to yield $x_{t-1}$ at each step. MPGD also provides an optional manifold projection step which utilizes pre-trained autoencoders to ensure $\hat{x}_{0}$ remains on the data manifold. For a fair comparison, we only consider MPGD without manifold projection in this work.

#### Loss Guided Diffusion (LGD).

LGD [] aims to improve the approximation of $p(y|x_{0})$ at each step along the sampling trajectory via a Monte Carlo approach. Viewing $p(y|\hat{x}_{0})$ in DPS as a delta distribution approximation of $p(y|x_{0})$ about $\hat{x}_{0}$, LGD instead computes the log-mean-exponential of $p(y|\hat{x}^{(m)}_{0})$ for $m=1,\dots,M$ perturbed copies of $\hat{x}_{0}$. That is, $p(\hat{x}_{0}|x_{t})\sim\mathcal{N}(\hat{x}_{0}(x_{t}),r_{t}^{2}I)$, where $r_{t}=\beta_{t}/\sqrt{1+\beta_{t}^{2}}$. The weighted (by $\zeta$) Monte Carlo gradient $\nabla_{x_{t}}\log\left(\frac{1}{M}\sum_{m=1}^{M}\exp\left(p\left(y|\hat{x}^{(
m)}_{0}\right)\right)\right)$ is then used to adjust $x_{t-1}$, as in DPS.

### B.3 Degradations and forward operators

In all experiments, pixel values are scaled from [-1, 1] (as in []) before application of forward operators. The details of the measurement likelihoods corresponding to each forward operator are outlined below. All methods use $\sigma=0.1$, except for phase retrieval, which uses $\sigma=0.05$.

#### Super-resolution.

The super-resolution task is defined by the following measurement likelihood:

$$ $y\sim\mathcal{N}(y|{\rm AvgPool}_{f}(x),\sigma^{2}I),$ $$

where ${\rm AvgPool}$ represents 2D average pooling by a factor $f$.

#### Gaussian deblur.

Gaussian blur is defined by a block Hankel matrix $C^{\psi}$ representing convolution of $x$ with kernel $\psi$:

$$ $y\sim\mathcal{N}(y|C^{\psi}x,\sigma^{2}I).$ $$

We consider a 61 x 61 Gaussian kernel with standard deviation of 3.0, as in [].

#### Inpainting.

The measurement likelihood corresponding to $p\%$ inpainting is a function of a mask $P$ with (1-$p$)% uniformly random 0 values:

$$ $y\sim\mathcal{N}(y|Px,\sigma^{2}I).$ $$

#### Nonlinear deblur.

Following [], the forward nonlinear blur operator is a pre-trained neural network $\mathcal{F}_{\phi}$ to approximate the integration of non-blurry images over a short time frame given a single sharp image []. Therefore, the measurement likelihood is as follows:

$$ $y\sim\mathcal{N}(y|\mathcal{F}_{\phi}(x),\sigma^{2}I).$ $$

#### Phase retrieval.

The forward operator of the phase retrieval task takes the absolute value of the 2D Discrete Fourier Transform $F$ applied to $x$: $|Fx|$. However, since this task is known to be highly ill-posed [, ], an oversampling matrix $P$ is also applied (with oversampling ratio 1 in this work):

$$ $y\sim\mathcal{N}(y||FPx|,\sigma^{2}I).$ $$

#### High dynamic range reconstruction.

In the HDR forward model, pixel values are scaled by a factor of 2 before truncation back to the range [-1, 1]. Therefore, the measurement likelihood is as follows:

$$ $y\sim\mathcal{N}(y|{\rm clip}(2x,-1,1),\sigma^{2}I),$ $$

where ${\rm clip}(\cdot,-1,1)$ truncates all input values to the range [-1, 1].

## Appendix C Additional experiments

#### Numerical SDE solver comparison.

Alternative numerical methods to EM (11) can be applied to discretize the Langevin dynamics SDE, such as the exponential integrator (EI) []. The EI scheme discretizes the nonlinear drift term $g^{i}=\left.\nabla_{x_{1}}L_{y}(x_{0})\right|_{x_{1}=z^{i}}$ and integrates the continuous-time dynamics arising from the linear term:

$$ $z^{i+1}=e^{-\tau}z^{i}-(1-e^{-\tau})g^{i}+\sqrt{1-e^{-2\tau}}\xi^{i},$ $$

where $\xi^{i}\sim\mathcal{N}(0,I)$. In Table A.4, quantitative comparison between our method using EM versus EI is shown on generating 10 samples for 100 images from the LSUN-Bedroom validation dataset, where the forward operator is nonlinear blurring. The same hyper-parameters are used for both methods, which are outlined in Table A.1. In this case, there is a marginal improvement in most metrics when using the EI scheme.

**Table A.4: Comparison between our method with EM and EI integration on the nonlinear deblur task on LSUN-Bedroom (256 x 256).**
| Method | PSNR $\uparrow$ | SSIM $\uparrow$ | LPIPS $\downarrow$ | FID $\downarrow$ |
| --- | --- | --- | --- | --- |
| Ours-EM(1-step) | 20.3 | 0.566 | 0.440 | 76.7 |
| Ours-EM(2-step) | 18.7 | 0.501 | 0.492 | 73.3 |
| Ours-EI(1-step) | 20.5 | 0.569 | 0.437 | 76.3 |
| Ours-EI(2-step) | 18.7 | 0.504 | 0.491 | 74.2 |

## Appendix D Additional qualitative results

Visualizations of additional reconstructions from our method corresponding to the linear and nonlinear experiments from Section 6.1 can be found in Figures A.1, A.2, A.3, A.4, and A.5. Additionally, diverse sets of samples from our one-step / two-step CM method corresponding to the experiments of Section 6.2 are visualized in Figures A.6, A.7, A.8, A.9, A.10, and A.11. Finally, diverse samples via the linear tasks on ImageNet (64 x 64) are shown in Figures A.12, A.13, and A.14. In these experiments, we use the one-step CM sampler with the same hyper-parameters as in Table A.2, but with $\tau=4\times 10^{-4}$ for inpainting, $\tau=9\times 10^{-4}$ for super-resolution, and $\tau=5\times 10^{-5}$ for Gaussian deblur.

Figure: Figure A.1: Additional image reconstructions for inpainting (left) and 8x super-resolution (right) on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/supp_plots-inpaint_supp_and_SR8x_fidelity.jpg

Figure: Figure A.2: Additional image reconstructions for Gaussian Deblurring on LSUN-Bedroom (256 x 256) (left) and ImageNet (64 x 64) (right).
Refer to caption: extracted/5896580/supp_plots-GD_supp_fidelity.jpg

Figure: Figure A.3: Additional image reconstructions for inpainting (left) and 4x super-resolution (right) on ImageNet (64 x 64).
Refer to caption: extracted/5896580/imnet_inpaint_sr_supp.jpg

Figure: Figure A.4: Additional image reconstructions for nonlinear deblur (left) and HDR reconstruction (right) on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/nlblur_hdr_supp.jpg

Figure: Figure A.5: Additional image reconstructions for phase retrieval on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/phase_supp.jpg

Figure: Figure A.6: Additional sets of samples for Inpainting (10%) on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/supp_plots-Inpaint_diversity.jpg

Figure: Figure A.7: Additional sets of samples for SR (8x) on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/supp_plots-SR8x_diversity.jpg

Figure: Figure A.8: Additional sets of samples for SR (8x) on LSUN-Bedroom (256 x 256) for 2-step method.
Refer to caption: extracted/5896580/supp_plots-GD_diversity.jpg

Figure: Figure A.9: Additional sets of samples for nonlinear deblur on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/nlblur_div_supp.jpg

Figure: Figure A.10: Additional sets of samples for HDR reconstruction on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/hdr_div_supp.jpg

Figure: Figure A.11: Additional sets of samples for phase retrieval on LSUN-Bedroom (256 x 256).
Refer to caption: extracted/5896580/phase_div_supp.jpg

Figure: Figure A.12: Sets of samples for 20% inpainting on ImageNet (64 x 64).
Refer to caption: extracted/5896580/imnet_inpaint_div_supp.jpg

Figure: Figure A.13: Sets of samples for 4x super-resolution on ImageNet (64 x 64).
Refer to caption: extracted/5896580/imnet_sr_div_supp.jpg

Figure: Figure A.14: Sets of samples for Gaussian deblurring on ImageNet (64 x 64).
Refer to caption: extracted/5896580/imnet_gblur_div_supp.jpg

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: * Can you show the inpainting mask of samples in Figure 1 and Figure 5 in the updated draft?
soundness: 3
strengths: * The paper is well written and easy to follow. 
* Projecting the posterior sampling to a noise space is a good idea, as the posterior distribution in the noise space is close to a single modal Gaussian distribution, that is more friendly to MCMC.
* Theoretical analysis has been provided for the proposed method.
summary: The paper proposes a method to do posterior sampling of samples, given certain conditions or partial information of the samples. The method assume there is a deterministic mapping between Gaussian noise and data, parametrized by e.g. a consistency model or a flow-based model. The posterior sampling is projected to and happens in the noise space with Langevin dynamics, and then the noise is projected back to the data space. Empirical results show that the proposed method leads to diverse and high quality samples in solving linear and nonlinear inverse problems.
weaknesses: * My main concern is that the paper fails to put itself to the current position in the literature. This type of projecting MCMC sampling to a more MCMC friendly noise space has been well established in [1], and has been later on adapted to generative modeling regime in e.g. [2, 3], with the deterministic mapping being a VAE or a flow-based model. The contribution of this paper, positioned in these literation, is that it adapted the sampling to a posterior distribution, and leverages a CM with fixed noise as the deterministic mapping. In that case, I think the novelty is limited. 

* The paper claims that the accumulation of samples leads to diverse samples. This needs to be further justified by analyzing the convergence behavior of the sampling chains. How do you make sure the samples from adjacent sampling steps are not correlated to each other (this could be guaranteed by other baseline methods that always start the sampling from independently sampled noise).

* The method assume that the sampling from the CM model is with fixed noise, which is lack of justification, and might partially explain why the sampling quality is suboptimal. 

* Why are 1-step CM results better than 2-step CM empirically in general? 

* Empirical results are not convincing enough. E.g., some samples in the top right row of figure 1 look oversaturated, which might indicate the sampling chain is not stable or mixing. 


[1] NeuTra-lizing Bad Geometry in Hamiltonian Monte Carlo Using Neural Transport

[2] VAEBM: A Symbiosis between Variational Autoencoders and Energy-based Models

[3] MCMC Should Mix: Learning Energy-Based Model with Neural Transport Latent Space MCMC

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: 1. What is the difference between the proposed method with an Energy-Based Model?
2. Except for image reconstruction, could the proposed model be applied in other task?
soundness: 2
strengths: This paper has proposed an efficient posterior sampling method that could avoid running the full sampling chain.
summary: This paper has proposed an efficient posterior sampling by simulating Langevin dynamics with a pre-trained generative model.
weaknesses: 1. The proposed method of posterior sampling using Langevin dynamics has been extensively studied and applied within the context of energy-based models. As a result, this approach lacks sufficient novelty for this work, given the established body of research already dedicated to similar methodologies. 
2. The proposed method of sampling by posterior sampling by Langevin dynamics has been early explored in many EBM works such as [1] - [4] for multiple kinds of tasks, such as image generation, translation and saliency prediction. While there is no discussion about these work either in background or experiments.

3. Efficient sampling based on molecular dynamics have been well studies in many literatures such as [5] - [7], while these study lack of discussion in the paper.
4. The experiment is limited to the image reconstruction tasks, which raises concerns about the method’s generalizability to other domains or applications. Without testing on a broader range of scenarios, it is difficult to assess the model’s robustness and adaptability to different types of data or tasks.

[1] Xie, Jianwen, et al. "A theory of generative convnet." International conference on machine learning. PMLR, 2016.

[2] Xie, Jianwen, et al. "Cooperative learning of energy-based model and latent variable model via MCMC teaching." Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 32. No. 1. 2018.

[3] Zhang, Jing, et al. "Learning generative vision transformer with energy-based latent space for saliency prediction." Advances in Neural Information Processing Systems 34 (2021): 15448-15463.

[4] Zhao, Yang, and Changyou Chen. "Unpaired image-to-image translation via latent energy transport." Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2021.

[5] Gao, Ruiqi, et al. "Learning energy-based models by diffusion recovery likelihood." arXiv preprint arXiv:2012.08125 (2020).

[6] Gao, Ruiqi, et al. "Flow contrastive estimation of energy-based models." Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2020.

[7] Zhu, Yaxuan, et al. "Learning energy-based models by cooperative diffusion recovery likelihood." arXiv preprint arXiv:2309.05153 (2023).

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: As stated in weaknesses.
soundness: 1
strengths: - Well-motivated problem, the authors did a good literature review that lists relevant works.
- The derviation of the framework is based on establishing theories of SMC/denoising diffusion models. A theoretical analysis is always welcomed.
- The method is relatively straight-forward and easy to implement, and work on both linear and non-linear inverse problem settings.
summary: This paper introduced a new diffusion solver for inverse problem by using posterior sampling in the noise space. The method is training-free, and, therefore, is more efficient compared to some other related training-based frameworks. The authors also provide theoretical guarantees for the approximation error in total variation distance; along with empirical benchmarks to demonstrate the improvements in realistic image restoration/inpainting/super-resolution task.
weaknesses: - Huge doubt about practical performance (inference time): while the authors report low NFE, the total runtime of the sampling framework is not reported. The backpropagation through the whole pretrained consistency model $\phi$ for each steps are costly in both memory and computational time. I think it will be fairer to compare total wall-time with other baselines instead of just listing the NFEs as stated in the paper 
- Unclear about the advantages of the proposed methods vs. baselines used in the benchmark: I also disagree on the argument of fair comparison to replace diffusion backbone by consistency models (CM) backbone to DPS and LGD  make them as fair as current CM-based framework in this paper (written around line 370-374). If DPS-DM and LGD-DM make faster inference *and* better quality reconstructed images, is it necessary calling them unfair comparison?
- Is the method completely training free? One of the key trick to make this method works, IMO, is the warm-start of the initial noise $z^0$, detailed around line 294-300, and on Appendix B.1. This is also related to the above point: the authors should also take into account runtime of the warm-start step and report it as requested on the point above.
- Missing strong baseline on linear inverse problem: the authors should include FPS (Dou & Song 2024) to the linear inverse posterior sampling baselines.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: I don't have specific questions and overall I recognize the contributions of this paper while I feel the novelty is not significant. So I prefer a boarderline accept.
soundness: 3
strengths: - The paper is well-written and easy to follow. The idea is simple but effective.
- The leveraging of consistency models significantly reduces the sampling costs, which is better than the previous data-prediction diffusion models.
- The experiments are solid.
summary: This paper proposes to use Langevin dynamics defined in latent space to perform posterior sampling of diffusion models, where the likelihood is defined in the clean data space. To evaluate the energy function, this paper proposes to use consistency models for significantly reducing the sampling cost. This paper also proposes a theoretical guarantee with mild assumptions. The empirical results show the effectiveness of the proposed method.
weaknesses: This paper doesn't have major weakness, though there are two minor aspects:

- The method is slightly lack of novelty since the combination of consistency models and Lagevin dynamics is a direct generation of previous method such as DPS.
- The proposed method can only tackle with simple likelihood functions such as inpainting / deblurring / super resolution. It is unclear whether the proposed method can do complicated posterior inference in traditional Bayesian inference settings.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: **Q1** : in Figure 2, how do the authors check if posterior samples are distinct? For instance, for DPS, do distinct posterior samples mean samples initialized from distinct prior noise?

**Q2** : what are the unconditional FIDs / Inception Scores of diffusion models and consistency models used in Section 6?

**Q3** : have the authors tried using other types of gradient-based MCMC such as Hamiltonian Monte Carlo (HMC)? HMC can mix faster than Langevin Monte Carlo, so perhaps it could be combined with adjoint method for diffusion prior based posterior sampling.
soundness: 2
strengths: **Originality** : to the best of my knowledge, this paper is the first to perform posterior sampling in the latent space of consistency models.

**Quality** : the proposed method is theoretically sound, and claims for diverse posterior sampling is supported with experiments on ImageNet 64x64 and LSUN bedroom 256x256.

**Clarity** : the proposed method and main results are clearly presented. I had no problem following the exposition.

**Significance** : while current diffusion-based image restoration methods often require large number of function evaluations, this paper introduces some ideas for fast image restoration.
summary: This paper proposes a posterior sampling scheme for solving inverse problems in the latent space of consistency models.
weaknesses: **Limited Originality** : running MCMC in the latent space of one-step pushforward generative models is not new. Similar ideas are already explored in works such as [1,2] (none of which are cited in this paper). I also feel the paper is limited in its technical novelty, as it does not provide any insight into efficient posterior sampling in the latent space of more challenging multi-step priors such as diffusion models or hierarchical variational autoencoders.

**Overstatements** : I feel that the paper overstates the extent of its contributions. First of all, the proposed method *is not applicable to all types of generative priors*, but is compatible only with implicit generative models / pushforward generative models. Second, the proposed algorithm is, practically speaking, *not applicable to diffusion priors*, as it would require backpropagation through up to thousands to neural net compositions. The authors do mention adjoint sensitivity methods for backprop through diffusion in Section 5, but experiment results are not presented. Perhaps it is because it is too computationally expensive to run multiple steps of Langevin MCMC until mixing occurs.

**Missing Ablations** : the paper is missing ablations w.r.t. design choices in Section 5. For instance, how does the performance vary w.r.t. Euler Maruyama step size $\tau$? Is optimal $\tau$ consistent across dataset and data dimension?

[1] Neutra-lizing Bad Geometry in Hamiltonian Monte Carlo using Neural Transport

[2] MCMC Should Mix: Learning Energy-based Model with Neural Transport Latent Space MCMC