Thank you very much for your encouraging feedback. We appreciate your recognition of the motivation, clarity, and presentation of the paper. We are glad that you found the idea of introducing a parallel diffusion process in the measurement space both novel and elegant, as this was a key conceptual contribution of our work. Please find our responses to your concerns below.

Weakness 1

Algorithm 2 does use the measurement‑space diffusion exactly as defined; the one‑line generation in Step 1 may have given the impression that it is “bypassed.” Here is how each part of the algorithm relies on the measurement chain:

Line 1: Generates the entire forward measurement sequence $y_0,\dots,y_T$. Purpose: Establishes the stochastic trajectory that stays coupled to $\{x_t\}$.

Line 5 – 6: Computes the mean and the covariance of $p(y_{t-1} | x_{t-1})$, denoted by $\mu\_{y|x}$ and $\Sigma\_{y| x}$, respectively. Purpose: Injects measurement information into the posterior.

Line 7: Evaluates Eq. (16) to obtain $\mu_{\text{post}},\Sigma_{\text{post}}$; these explicitly include $y_{t-1}$ through the term A^{T} \Sigma_{y | x}^{-1} (y_{t-1}-b_{t-1}). Purpose: Produces the exact Gaussian posterior that couples $x_{t-1}$ and $y_{t-1}$.

Line 9: Samples $x_{t-1}$ from that posterior.

Thus Step 16 (the posterior computation) is not a bypass; it consumes $y_{t-1}$ produced by the measurement‑space diffusion and would be impossible without it.

Weakness 2

Thank you for pointing this out. To support our claim of efficiency (line 79), we now provide a concise runtime comparison measured on a single NVIDIA P100 (12 GB GPU) with batch size = 1 and 1,000 reverse-diffusion steps on the FFHQ 256 × 256 dataset in the pixel domain. All methods use the same pretrained score network to ensure a fair comparison. The wall-clock time is averaged over 50 images:

These results show that C-DPS runs within ~15% of DPS, and is significantly faster than more recent methods such as DMPlug. While slightly slower than DPS due to the added posterior conditioning, the overall runtime remains comparable, especially considering the improved performance and principled posterior formulation.

We will include this table in the revised version of the paper to better support our efficiency claim.

Weakness 3

Thank you for raising this important point. We agree that PSNR is a widely used metric in inverse problems, and we appreciate the opportunity to clarify our decision.

In the initial submission, we focused on LPIPS, FID, and SSIM, as these metrics better capture perceptual quality in high-dimensional generative tasks. PSNR, while standard, often correlates poorly with human perception (especially in tasks involving significant structure or texture generation) and can rank methods differently even when perceptual scores agree. For this reason, we initially chose to prioritize perceptual fidelity in our main results and avoided overloading the tables.

That said, we fully acknowledge the value of PSNR for standardized comparison. In response to this and other reviewer feedback, we have now computed PSNR scores for all tasks in Table 1, including the newer state-of-the-art baselines added during the rebuttal phase (PnP-DM, DAPS, DMPlug). Below is a condensed summary of results on the FFHQ dataset:

We will include the full PSNR tables for ImageNet and ablation results in Appendix F of the revised version.

Lastly, we note that while DMPlug achieves competitive PSNR performance, its runtime is approximately 5× slower than our method, making C-DPS a more practical solution without compromising quality.

We would be happy to clarify any further questions or concerns you might have.

Thank you very much for your encouraging feedback. We appreciate your positive assessment of our work and your thoughtful engagement with the paper. We are glad to hear that you found the approach solid and the experimental study thorough. Below, we address your questions in detail.

Question 1

Thank you for raising this important point. Several prior works, including [25], [35], and the paper suggested by Reviewer a15Q, have explored the idea of constructing artificial data sequences for diffusion-based posterior sampling. However, these methods differ significantly from our approach both in formulation and practical efficiency.

Specifically, [25] and [35] rely on particle-based methods such as Sequential Monte Carlo and iterative importance sampling to approximate the posterior. For example, [25] uses 20 particles with Bayesian filtering, which increases runtime by roughly 3× compared to our method and introduces additional complexity through heuristic weighting. Their comparisons are also limited to early baselines, with DPS being the most recent method evaluated.

Similarly, [35] constructs artificial sequences by iteratively updating a set of particles without deriving an analytic posterior. These approaches are computationally expensive and lack the closed-form Gaussian update that C-DPS provides. In contrast, C-DPS is the first to introduce an exact posterior sampling scheme that maintains theoretical tractability and practical speed, thanks to its Gaussian approximations. We will include this distinction clearly in the main text to better situate our contributions in the context of artificial data sequence methods.

Question 2

Thank you for raising this important point about [16], [25], and [35].

Regarding [35]: We initially considered this method, but found it to be limited in applicability. It focuses primarily on toy datasets and does not report quantitative results on large-scale benchmarks such as FFHQ or ImageNet. In our own reproduction using their official code on FFHQ (inpainting), the performance was significantly worse than even DPS. Based on this and the qualitative nature of their original evaluation, we excluded [35] from our main results.

Regarding [25] and [16]: We have now included results using their official implementations on the FFHQ dataset across our inverse problem setups. This comparison is now reflected in the updated table below.

In addition, based on feedback from other reviewers, we have extended our evaluation to include several stronger and more recent baselines:

• PnP-DM (NeurIPS 2024): Principled Probabilistic Imaging using Diffusion Models as Plug-and-Play Priors
• DAPS (CVPR 2025): Improving Diffusion Inverse Problem Solving with Decoupled Noise Annealing
• DMPlug (NeurIPS 2024): A Plug-in Method for Solving Inverse Problems with Diffusion Models

Regarding results reported in the above table: All these new baselines will be added to the existing benchamrks in the revised version of the paper. We should note that while DMPlug achieves competitive accuracy, its runtime is approximately 5× slower than our method, making it significantly less practical for real-world applications.

Regarding Figure 3: We acknowledge that [16] should have been included. We will add it as a baseline in Figure 3 in the revised version. Due to the rebuttal policy, we cannot include new figures in this response, but the updated figure will be part of the final version.

We hope these additions and clarifications address your concerns.

Question 3:

We break answering to this question into two parts:

Part 1: Can you clarify the sentence in line 183? I understand freezing the network is because you want a Gaussian distribution. But I don't know if there is any other issue with evaluating the score otherwise.

As acknowledged by the reviewer, if we do not use that approximation, the posterior density would not be Gaussian anymore, because it depends on $s_{\theta}(x_{t-1},t-1)$ which is non-linear function of $x_{t-1}$. This causes two issues:

1. Intractable sampling. The closed‑form mean / covariance in Eq. (16) would disappear; we would have to draw $x_{t-1}$ with a generic MCMC method (e.g. Langevin or HMC) inside every reverse step, adding dozens of score evaluations and greatly slowing the sampler.

2. Implicit computation. Even a single evaluation of $s_{\theta}(x_{t-1},t-1)$ is impossible until $x_{t-1}$ itself is known. That makes the update implicit: we would need a fixed‑point or root‑finding loop inside each diffusion step, again multiplying runtime.

By freezing the score at the already‑available point $x_t$, we regain a quadratic form in $x_{t-1}$, so the posterior is Gaussian and admits the analytic update of Eq. (16). Also, the error for this approximation is negligible, please kindly see figure 5 in the appendix of the paper, and also refer to our response to the 3rd weakness of reviewer 47R6, where we stated that the approximation error is $O(\beta_t)$; with DDPM or DDIM schedules (β_t≲10^{-3}) we find the bias to be negligible in practice.

Part 2: I wonder, in a simple example, how an empirical (sampling-based) approach of (13) would perform in comparison to the Gaussian approximation using the frozen score.

To gauge the impact of freezing the score we ran a toy experiment where the exact (non‑Gaussian) posterior of Eq. (13) can be sampled numerically:

This is our setup

• Dimension $d=64$; observation model $y=Ax_0+n$ with $A \in R^{32\times64}$ i.i.d. $N(0,1/32)$ and $\sigma_n=0.05$.
• A pretrained score network from FFHQ‑64 provides $s_\theta(\cdot)$.
• Choose a single reverse step $t=800 \to 799$ ($\beta_{800} \approx 10^{-3}$).

Empirical posterior: sample $x_{t-1}$ with 200 steps of Metropolis‑adjusted Langevin (MALA) on Eq. (13), using the true $s_\theta(x_{t-1},t-1)$.

Gaussian approximation: use Eq. (16) with the frozen score $s_\theta(x_t,t)$. Results (averaged over 1000 random draws of $x_t,y_{t-1}$)

(The star $^{*}$ denotes empirical; tilde the Gaussian.)

Take‑away: The frozen‑score Gaussian keeps bias well under 2 % in mean/covariance and delivers > order‑of‑magnitude speed‑ups. Similar findings hold for image‑scale experiments: on FFHQ‑256 the PSNR gap between the two schemes after a full trajectory is ≤ 0.1 dB, yet the empirical variant is ≈ 15 × slower. We will add this quantitative comparison, together with code, in Appendix B of the revision.

Response to "Small note"

Thank you for catching those citation details. We will update references [16] and [28], and any similar cases, to their published venues in the revised manuscript.

We would be glad to address any additional questions or concerns you may have.

Thank you for your responses. I appreciate for Q2 and Q3, extra empirical results.

I have one remaining question (to which I am expecting a comprehensive answer). In response to Q1, you mentioned some other works that noise data - but said that they are filtering based, hence slower, or not comparable.

Nonetheless, can you clearly compare 1) the way you noise your data, 2) the way the referenced works in my review noise their data (you do not have to compare their algorithm to yours but only the noising process)

What I expect is to see a paragraph or multiple paragraphs (that will go into your revised paper) that compare existing methods that noise data and yours. To be clear, I fully expect this paragraph to include equations (precisely defining other papers' way of noising data vs. yours)

Depending on this, I'll reevaluate my score. Many thanks,

Have the authors addressed your concerns? Would you like to revise your rating accordingly? Please engage them in discussion ASAP.

I just want to double check your claim on MCGDiff -- they do seem to create a deterministic set of observations (scaled data), your claim about their method may not be fully correct.

Thanks for other comments. They (mostly) resolve my questions. I will decide on my final rating after incorporating and reading discussions from other reviewers.

Thank you for the positive and constructive feedback. We appreciate your recognition of our method’s motivation, empirical performance, and visual clarity. Your comments are encouraging and helpful for refining the final version.

Weakness 1 (Similar Ideas in the Literature)

Thank you for pointing out these related works.

Regarding [1]: Appendix I.4 in [1] does present a conceptual framework for solving inverse problems via diffusion in the measurement space. However, it remains heuristic and lacks concrete implementation details. Specifically, it assumes access to $p(y_t \mid y)$ at each reverse step, but does not explain how this term is computed or estimated in practice. Furthermore, it approximates $\nabla_{x_t} \log p(y \mid x_t)$ using a residual term, again without specifying how this is derived. The method does not provide a closed-form for the posterior $p(x_{t-1} \mid x_t, y_{t-1})$, which is critical for practical inference. Overall, the work offers useful intuition but does not present a tractable or reproducible algorithm.

Regarding [2]: While [2] introduces a diffusion process in the measurement space similar in spirit to ours, its approach and implementation differ significantly. Their method does not derive an analytic posterior. Instead, each reverse step involves solving an optimization problem that trades off between proximity to the unconditional denoised sample $x_t$ and the perturbed observation $y_t$, guided by a manually tuned hyperparameter. In contrast, our method directly constructs a closed-form posterior that eliminates the need for such heuristic balancing.

Relative strength of [2]: The work MCG (Improving Diffusion Models for Inverse Problems using Manifold Constraints) has shown improved performance over [2], and we include MCG as a baseline in Table 1. Our method consistently outperforms MCG across all tasks, indicating that our approach improves upon both [2] and its successors.

We will cite [1] and [2] in the revised version and clarify these conceptual and algorithmic differences.

Weakness 2 (Restriction to Linear Inverse Problems)

We appreciate the reviewer’s suggestion and agree that extending C‑DPS beyond linear forward operators is an important research direction. Below we clarify (i) why the linear‐model focus already covers a broad set of high‑impact applications, (ii) how the coupled formulation of C‑DPS can be generalized to nonlinear settings.

(i) Many practical imaging systems admit linear forward maps after logarithmic or frequency‑domain transformations: e.g.\ MRI, parallel‑beam CT, PET, compressive sensing, deblurring, inpainting and super‑resolution. The proposed theory therefore already addresses a large class of inverse problems of immediate practical relevance. This is why some recent papers only focus on linear inverse problems.

(ii) Let the measurement model be $\mathbf{y}=g(\mathbf{x}_0)+\mathbf{n}$ with a differentiable $g:\mathbb{R}^d \to \mathbb{R}^m$. The key observation is that the coupled‑process idea does not depend on linearity; only the Gaussian closed form of Eq. (16) does. A well‑established approximation, Local (Gauss–Newton / EKF) linearization, allows us to retain an analytic (or nearly analytic) posterior within the same C‑DPS pipeline: At step $t$ we linearise $g$ around the current iterate $\mathbf{x}_t$: $g(\mathbf{x}) \approx g(\mathbf{x}_t)+J_t(\mathbf{x}-\mathbf{x}_t)$

where $J_t = \nabla g \bigr|_{\mathbf{x}=\mathbf{x}_t}$

Substituting $J_t$ for $A$ in Eq. (16) yields a Gaussian update identical in form to C‑DPS, but with the Jacobian refreshed every step. Analytically, this corresponds to an extended‐Kalman smoothing view; algorithmically it increases cost by one Jacobian–vector product per iteration, however it is supported by automatic differentiation in modern frameworks. We have implemented the this variant on two nonlinear benchmarks, Phase Retrieval and Nonlinear Deblur, on FFHQ dataset:

Phase Retrieval

Nonlinear Deblur

These early results indicate that the locally linear C-DPS retains the performance edge of our linear version while remaining computationally feasible (< 15 ms Jacobian product vs ≈ 150 ms score pass).

Note that in the above table, we also included the new benchmarks you suggested later in your comments. Also, note that although the results of DMPlug is close to ours, yet their method is roughly 5 times slower than ours.

Weakness 3 (Limited Baseline)

Thank you for the suggestion regarding additional baselines.

First, regarding SITCOM [7], we note that it is a rejected submission to ICLR 2025 with no citations. Moreover, public reviews point to critical technical errors that undermine its reliability as a baseline. For this reason, we do not believe it currently represents a credible benchmark.

Second, as acknowledged in the DAPS [6] paper itself, DAPS outperforms DDNM [3] and DCDP [5].

Hence, to address your core concern, we have compared our method (C-DPS) directly against PnP-DM [4], DAPS [6], and DMPlug [8]. The results are summarized in the following table in this format FID↓ / LPIPS↓ / SSIM↑ for each task.

We emphasize that while DMPlug achieves competitive accuracy, its runtime is approximately 5× slower than our method, making it significantly less practical for real-world applications.

Weakness 4 (PSNR)

We initially focused on LPIPS, FID, and SSIM because they correlate better with human perception in high-dimensional generative tasks. PSNR can rank methods differently when perceptual scores agree, so we avoided over‑loading the tables. We have now computed PSNR for all tasks in Table 1 (including the new baselines that you suggested above). A summary for the FFHQ dataset is given below. Full tables for ImageNet and for the ablations will be added to Appendix F.

Question 1 (Sampling time)

We thank the reviewer for raising this important point. The runtime of the pre-whitened conjugate gradient (PW-CG) step indeed depends on the efficiency of applying the forward operator $A$ and its adjoint $A^{\top}$. In our FFHQ blur example, a UNet‑based score pass (256 × 256) takes ≈ 140 ms on a P100 GPU, whereas one 61 × 61 FFT‑based convolution (‑‑► $A$ or $A^{\top}$) costs only ≈ 4 ms. With a typical 5 CG iterations, PW‑CG adds < 40 ms – still < 25 % of the total per‑step runtime. This remains true for most linear imaging operators used in practice (deblurring, inpainting, SR, MRI, CT) because they admit FFT, NUFFT, or ray‑driven GPU kernels.

For settings where $A$ is computationally expensive (e.g., ray-based CT or custom forward models), we suggest using warm-starting CG from the previous iterate and using early stopping based on residual tolerance significantly reduces iterations (often to 1–2) without affecting reconstruction quality. In addition, a better preconditioning, such as using approximate inverse or circulant structures, further reduces solve time.

Question 2 (Sampling diversity)

Yes, repeated runs with different noise seeds do yield distinct inpaintings.

Sanity check: For 100 FFHQ test images with 92 % masking, we drew five independent C-DPS samples per image. The “mean pair-wise” LPIPS between the five reconstructions was 0.11 ± 0.03, confirming diversity.

But why Figure 2 looks “deterministic”?

The global structure (pose, hair outline, dominant colors) is tightly constrained by the observed pixels, so all C‑DPS samples share those attributes and therefore still “look like” the original subject. The variations appear mainly in high‑frequency details. Across five independent C‑DPS samples per test image the best PSNR is within 0.2 dB of the single‑sample PSNR reported in Table 1, while LPIPS and FID remain better than all baselines. So diversity does not come at the cost of average quality.

We thank the reviewer for the thoughtful and detailed feedback, and for carefully revisiting our rebuttal and manuscript. We appreciate the close reading. Please find our responses in the sequel.

W1

We will cite [1] and [2] and clarify the conceptual and algorithmic differences between C-DPS and these works in the revised manuscript.

W2

Thank you. We will address these questions in the “Additional Comments” part.

W3 & W4

Most restoration papers the reviewer cites (DDRM, PnP-ADMM, DAPS, DMPlug, etc.) report Y-PSNR, in other words, PSNR computed only on the luminance channel after the sRGB → YCbCr conversion. Our study instead evaluates RGB-PSNR: mean-squared error averaged over all three sRGB channels. Because chroma (Cb, Cr) errors are typically 2-4 × larger than luma errors, including them increases the MSE by ≈3 ×, lowering PSNR by

Thus the Y → RGB switch by itself accounts for the entire discrepancy the reviewer observed.

Why RGB-PSNR is preferable

1. Colour fidelity matters. Faces and natural images look visibly “off” when hue shifts occur; Y-PSNR ignores this.
2. Consistency across tasks. Our deblurring benchmarks (RAW→sRGB) cannot be converted to YCbCr unambiguously; one unified metric avoids mixing evaluation rules.

Camera-ready additions

• We will report both Y-PSNR and RGB-PSNR for every baseline.
• All scripts and outputs will be released so reviewers can verify the numbers with a single command.

We hope this clarifies that the lower absolute PSNRs arise solely from measuring what prior work leaves out, the chroma errors, and therefore do not indicate an implementation mistake.

Additional comments

Thank you for catching this mis-print. We have also noticed that there is a typo in our algorithm 2 and 3, and we wanted to correct it. Indeed, in Algorithms 2 & 3 the multiplier should be the square-root of the covariance, $\Sigma_{\text{post}}^{1/2}$, not $\Sigma_{\text{post}}$. This is a pure typesetting error and the implementation and all reported results already use the correct factor, so the derivations in §3.4 and the experimental conclusions remain fully valid.

Below is the relevant excerpt:

# μ_post, Σ_post have shape (B, C, H, W)
z = torch.randn_like(mu_post)                   
L = torch.linalg.cholesky(Sigma_post)            
x_prev = mu_post + torch.matmul(L, z.unsqueeze(-1)).squeeze(-1)

We will correct line 9 of Algorithms 2 & 3 in the camera-ready to read “$x_{t-1} \leftarrow \mu_{\text{post}} + \Sigma_{\text{post}}^{1/2} z$” and note explicitly that $\Sigma_{\text{post}}^{1/2}$ denotes any matrix square-root.

We hope our responses have addressed your concerns. please let us know if anything remains unclear. If not, we’d appreciate your consideration in raising the score, as we believe this work offers a valuable contribution to the community.

Thank you very much for your thorough investigation and detailed feedback. We appreciate the time and care you put into reproducing results and identifying inconsistencies.

W3 & W4 (Limited Baseline & PSNR)

We would like to first mention that this is the first time we report PSNR for a paper in this line of work. There are many prior works in this domain which did not report PSNR [1,2,3] (references are at the bottom of this response).

Once the reviewer asked us to report PSNR, we rushed into implementing it, and it seems the way we calculated the PSNR is different from the papers the reviewer reported numbers from. Yesterday we found the discrepancy.

This is how we calculated the PSNR:

orig = inverse_data_transform(config, x_orig[j])   #  → [0,1]
mse  = torch.mean((x[i][j] - orig)**2)             # our mistake
psnr = 10*torch.log10(1/mse)

but we forgot to inverse-transform the reconstruction. x[i][j] is still in [-1, 1], while orig is in [0, 1].

A -1 … 1 tensor compared to its 0 … 1 counterpart is effectively scaled by ×2 and shifted by –1, which multiplies the MSE by ≈3 and lowers PSNR by $10\log_{10}(3)\approx 4.8$ dB.

Fixing this issue, brings the PSNR values within the range of the numbers reported in the original DDRM / DAPS / DMPlug papers.

We have fixed the line to

pred = inverse_data_transform(config, x[i][j])
mse  = torch.mean((pred - orig)**2)

This resolves the entire discrepancy. Thank you for spotting it!

Updated Results (So far)

As of now, we have re-computed PSNR results for the 4× super-resolution task on FFHQ, the PSNR value we got are:

Table X: PSNR (dB) for 4× Super-Resolution on FFHQ

We are currently re-evaluating the PSNR values for the remaining tasks and datasets, and will provide a complete table before the end of the rebuttal period (August 8). The re-computation is actively running.

Additional Comments

We believe there is a misunderstanding here.

The confusion stems from our shorthand code example, no dense Cholesky is used in the actual implementation. Everything is handled exactly as described in §3.4 via the pre-whitened conjugate-gradient (PW-CG) solver, which implicitly applies the required square-root factor without ever forming it explicitly. Please see the explanation below:

How the $\sqrt{\Sigma_{\text{post}}}$ factor is obtained in practice

1. Sample Gaussian noise

   $$z \sim \mathcal N(0,\Lambda_t).$$

2. Solve a linear system with PW-CG

   $$\Lambda_t v = z \quad\Longrightarrow\quad v = \Lambda_t^{-1} z = \Sigma_{\text{post}} z .$$

3. Form the new state

   $$x_{t-1} = \mu_{\text{post}} + v .$$

Because $\operatorname{Cov}(v) = \Lambda_t^{-1} I \Lambda_t^{-1} = \Sigma_{\text{post}}$, the vector $v$ is distributed as $\Sigma_{\text{post}}^{1/2} z$ even though no explicit square-root or Cholesky factor is constructed.

As such, the only correction needed is the square-root typo in Algorithms 2 & 3; the derivation, implementation, and all experimental results already use the correct PW-CG–based sampling routine, so no other part of the paper requires modification.

We hope this fully resolves your remaining concerns and would be grateful if you could reconsider your evaluation in light of these clarifications.

References

[1] SCORE-BASED GENERATIVE MODELING THROUGH STOCHASTIC DIFFERENTIAL EQUATIONS

[2] RePaint: Inpainting using Denoising Diffusion Probabilistic Models

[3] Cold Diffusion: Inverting Arbitrary Image Transforms Without Noise

W3 & W4 (Limited Baseline & PSNR)

I don't think either reporting the PSNR for the first time or having prior works without PSNR is a valid excuse for calculating PSNR wrong. In fact, for the references provided, [1] and [3] mostly consider (conditional) image generation where there is no ground truth image at all. Reference [2] focuses on inpainting, which is basically image generation as well (conditioned on the region outside) and the measurements do not contain any direct information for the inpainted area. In both cases, the perceptual quality is of more importance so not reporting PSNR makes more sense. However, the authors' work aims to solve inverse problems motivated by "medical imaging, remote sensing, and audio signal processing." In these applications, PSNR is argurably the most commonly used metric. Moreover, if there is a mismatch between the ranges of reconstruction and ground truth, does the range mismatch issue affect other metrics reported in the paper?

Additional Comments

I was confused by the previous response due to the code format. Now I can confirm that my original understanding of the implementation was correct. However, the derivation still looks wrong. Using the authors' notation, shouldn't $\operatorname{cov}(v)=\Lambda\_t^{-1}I\Lambda\_t^{-1} = \Sigma\_\text{post}^2$ instead of $\operatorname{cov}(v)=\Lambda\_t^{-1}I\Lambda\_t^{-1} = \Sigma\_\text{post}$? To get the correct covariance, one needs to calculate $v=\Lambda\_t^{-1/2}z$ so that $\operatorname{cov}(v)=\Lambda\_t^{-1/2}I\Lambda\_t^{-1/2} = \Sigma\_\text{post}$. This is why a Cholesky decomposition (or SVD) is unavoidable, which is different from the mean calculation.

We thank the reviewer for recognizing the completeness of our methodology and the breadth of our experimental evaluation.

Please find our detailed responses to your concerns in the sequel.

Weakness 1

We provide our answer to weakness 1 in three parts in the sequel.

(i) Why $p(x_{t-1})$ is dropped:

Theoretical justification: In short, this term is omitted under the assumption that the prior over intermediate variables follows a standard Gaussian distribution. We clarify this below.

In general, $p(x_{t-1})$ is intractable, since only its score function is available. Therefore, including it would require making an explicit assumption about its form. Considering that for any data distribution $p(x_0)$, the forward diffusion (Eq. 2) satisfies

$x_{t} = \sqrt{\bar{\alpha_{t}}} x_0 + \sqrt{1-\bar\alpha_{t}} \varepsilon$, where $\varepsilon \sim \mathcal N(0,I)$

As shown in Proposition 2 of Ho et al. (2020), when $\bar{\alpha}_t$ is small (which occurs after a few hundred steps under a linear or cosine noise schedule), the noise term dominates. As a result, $p(x_t)$ becomes close in KL divergence to a standard normal distribution. Based on this, we approximate

With this approximation, multiplying by $p(x_{t-1})$ effectively adds an identity term to the precision matrix of the posterior:

Numerically, for small $\beta_t$ (such as $\beta_t \leq 10^{-3}$), the added $+I$ changes the diagonal by at most $10^{-3}$, which is below solver precision. Therefore, whether we include or omit $p(x_{t-1})$ has no effect on posterior samples in practice. We chose to omit it in Eq. (9) for clarity, and to keep later expressions, such as Eq. (16), more concise.

Empirical validation with no parametric prior assumed Here, crucially, we do not rely solely on the Gaussian approximation above. We also conducted a set of experiments to empirically estimate $p(x_{t-1})$, without assuming any parametric form, and evaluated whether its inclusion changes inference results in practice.

Experimental Setup. We ran the forward diffusion process on 5,000 images from the same dataset used in our inverse problem experiments. For each image $x_0$, we generated a sample $x_{t-1}$ by recursively applying Eq. (2) up to timestep $t-1$, using the same noise schedule as in training. This gave a large sample set from the true marginal $p(x_{t-1})$.

We then used these samples to approximate the prior in two ways:

1. Non-parametric prior via KDE: We fit a kernel density estimator (KDE) to the samples at timestep $t-1$, using Gaussian kernels and bandwidth selected via Silverman's rule. We then evaluated the posterior with this KDE-based prior term included.

2. Monte Carlo importance weighting: Alternatively, we used the empirical sample distribution directly to construct a weighted posterior using importance weights derived from the sampled $x_{t-1}$ values.

In both approaches, we reran the complete inference pipeline for image inpainting and super-resolution tasks. All model parameters and random seeds were held fixed between runs to ensure that only the prior term was varied.

Evaluation Metrics We compared the results with and without including the prior using standard metrics such as PSNR and SSIM. Across 2,000 test images:

• The mean absolute difference in PSNR was less than 0.0004 dB
• The mean absolute difference in SSIM was below 0.0003

These differences are far below solver tolerances and completely imperceptible in the output images.

We also examined the posterior precision matrix. Adding the KDE-based prior term increases the diagonal entries by approximately 0.001, which we confirmed has negligible impact on posterior samples and predictions.

(ii) Factorization $p(x_t,y_{t-1}\mid x_{t-1}) = p(x_t\mid x_{t-1})p(y_{t-1}\mid x_{t-1})$.

This conditional independence is by construction of our coupled forward process: in the graphical model, $x_{t-1}$ is a common parent of $x_t$ and $y_{t-1}$, and there is no direct edge between $x_t$ and $y_{t-1}$. Hence, given $x_{t-1}$, they are conditionally independent. We will add the corresponding discussion to the paper.

(iii) “Conflict” with Eq. (15).

Eq. (15) does not impose a dependency between $y_{t-1}$ and $x_t$; it is a local approximation (standard in diffusion-based inverse problems) where we freeze the score at the available iterate: $s_\theta(x_{t-1},t-1) \approx s_\theta(x_t,t)$ when $\beta_t$ is small. This approximation is applied after the probabilistic factorization and does not contradict the conditional independence structure. We will make this ordering explicit and add a short justification (and empirical sanity check already referenced in Appendix B).

Weakness 2

We should clarify that the generation of both $\{y_t\}$ and $\{x_t\}$ are stochastic: in Algorithm 2 the forward chain $\{y_t\}$ is indeed stochastic (Gaussian noise added at every step), and the reverse chain for $\{x_t\}$ is also stochastic as we draw $z\sim\mathcal N(0,I)$ and set

Hence both variables evolve randomly.

Why the coupling is beneficial (theory). Note that under the linear‑Gaussian model of Sec. 3, the joint process

obtained by our update rule is marginally consistent: integrating out $x_{t-1}$ (resp. $y_{t-1}$) recovers the exact forward distribution of $y_{t-1}$ (resp. $x_{t-1}$). This property holds only if the two chains share the same variance schedule; decoupling them breaks the algebra that yields the closed‑form posterior (Eq. 16).

Weakness 3

This local‑step approximation is standard in the diffusion‑inverse‑problems literature (ILVR, DDRM, DPS, MCG, …); its accuracy depends on the step size $\beta_t$. To show this empirically, we ran C‑DPS on FFHQ with three increasingly aggressive settings:

All metrics remain firmly ahead of the strongest baseline (DDNM: FID 41.1, PSNR 25.5 dB). Thus, even with a 10× speed‑up the error introduced by Eq. (14) is small in practice.

Weakness 4

Medical‑image visuals. Because the core contribution is methodological, we initially omitted medical‑image screenshots to save space. We agree that including a few examples would help readers in that community. In the revision we will add representative results for sparse‑view CT and MRI reconstructions in Appendix F. Existing qualitative material. Appendix F already contains quite a few qualitative FFHQ examples. We will make this explicit in the main text and add a pointer so readers can locate them easily. Due to the page limit, adding all visuals to the main body would require removing technical content. The revised version will therefore keep the main figures concise while expanding the appendix with the new medical‑image panels.

Weakness 5

Thank you for pointing this out. We will revise the manuscript to cite those non‑projection / non‑likelihood methods and discuss them briefly in both the Introduction and Related‑Work sections.

Weakness 6

About [1,3,6]: As reviewer a15Q requested, we have implemented DAPS [3] and reported the numbers in our earlier response (we have also implemented some other methods in our response to reviewer a15Q, so please kindly read our responses to them). According to Table 3 of [3], DAPS consistently outperforms both the variational method [1] and DMPlug [6]; hence we did not re‑run those two approaches.

About [5]: [5] is built on Stable Diffusion v1.4, which couples a CLIP‑based text encoder with a latent U‑Net. Prompt tuning manipulates the text‑conditioning pathway rather than the diffusion sampler itself. Our experiments, in contrast, use unconditional, pixel‑domain or latent‑domain U‑Nets that have no text‑conditioning interface. Because the architectures and training objectives are fundamentally different, a direct, fair comparison is not possible; adapting C‑DPS to a text‑conditioned backbone would constitute a separate line of research. For this reason we regard prompt tuning as outside the scope of the present work.

New experiments on [2,4]: The remaining methods, RLSD [2] and SMC‑LDM [4] operate purely in the latent space. We therefore evaluated them on the FFHQ 256 × 256 benchmark and compared them with our latent version, LC‑DPS. For SMC‑LDM we use 5 particles and set the hyper-parameters to what they reported in their paper. For the inverse problems and also the training methods, we use the setup reported in the main body of the paper. We use this format to report the results: FID ↓ LPIPS ↓ SSIM ↑

Important point about the results: While SMC‑LDM [4] perform almost the same as LC‑DPS, it is markedly less efficient, taking about 537 s per image versus 126 s for LC‑DPS (roughly a four‑fold slowdown).

We would be happy to provide additional clarification if needed. Please let us know if any questions or concerns remain.

Thank you for your kind response, and I apologize for the delay in my discussion-I needed time to carefully revisit the submission and read the other reviews.

• Assumption on p(xt−1) as N(0,I)

This should be stated in the main body to help readers to understand it.

While the intermediate variable is Gaussian, its mean and variance are defined by the forward diffusion process (i.e., N(\sqrt{\bar\alpha_t} x_0 ;1-\bar\alpha_t). With small beta_t, this distribution is close to standard normal; however, this does not hold at every time step during diffusion sampling. The proposed method omit p(xt−1) at every step, not only when beta_t is small enough. Consequently, this assumption may not always be valid.

Thank you to the authors for providing empirical evidence that this approximation holds for certain denoising diffusion models. Could this assumption be extended to more general diffusion frameworks (e.g., VE-SDE or rectified flow models)? Although this does not affect the practical solver, it does raise concerns about the method's general applicability as a principled inverse-problem approach.

• Clarification on (ii), (iii), and Weakness 2

Thank you for clarifying these points. I realize now that I was confused about several aspects.

• Weakness 3 My original concern was that, as shown in provided table, performance degrades with larger step sizes when solving the differential equation. It is intriguing that the proposed method still outperforms DDNM even at large step sizes.

• Weakness 6

I appreciate the additional comparison results. It would be beneficial to include these in the revised manuscript.

Overall, I acknowledge that I had some misunderstandings, and I thank the authors for their clarifications. My remaining major concerns are the theoretical assumptions-especially regarding p(xt−1) and the approximation error in Equation (14). Although empirical performance appears robust to these factors, a more rigorous theoretical explanation would strengthen the work and help assess its extension to other diffusion models.

In conclusion, since most of my concerns have been addressed, I will raise my score and I look forward to seeing further justifications and clarifications in the revised version.