Thank you for your positive feedback and insightful comments. We will make sure to take your invaluable suggestions into account in revising the paper and in our future work. Below we address the two questions you raise.

Connections with projection-type approaches.

• Thank you for your thoughtful question. We are willing to clarify our connection and difference with projection-type approaches as follows.

While our algorithm draws inspiration from projection-type approaches, the essential difference between our algorithm and these approaches persists even for linear problems. Roughly speaking, our proximal consistency step can be reduced (cf. Line 274-277 in our paper) to a slightly more complicated, yet theoretically rigorous version of the (proximal) projection step in these approaches for linear problems, but the same cannot be said about our denoising diffusion step.

A manifestation of such essential difference is that our denoising diffusion step requires multiple steps of reverse diffusion in a single iteration of outer loop, and has an accurate, simple expression for the (stationary) distribution it aims to sample from. In contrast, for projection-type approaches, it is typical that only one step of reverse diffusion is run in a single iteration of outer loop, and the distribution of its output, even in the ideal continuous setting, cannot be expressed in a simple way, which significantly hinders theoretical analysis. From a big picture, by running more steps of reverse diffusion every time after we do proximal consistency sampling, we stabilize the process since we better compensate for the drifts from the data manifold induced by (proximal) projection onto the null space of consistency equations. This is the heuristic reason for the improved robustness of our algorithm, which is backed by rigorous mathematical proofs.

Order of proximal consistency step and denoising diffusion step.

• Thank you for your sharp observation. The first reason for such an order is that the algorithm would have visually better outputs if ending with the denoising diffusion step instead of the proximal consistency step, as denoising diffusion by design outputs visually satisfactory images while proximal consistency sampler only cares about consistency. Such a difference was not present in previous projection-type approaches, as the denoising stage there runs only one step of reverse diffusion, thus it is insignificant whether it comes as the penultimate step or the last step. This again constitutes an example of the difference of our algorithm. The second reason is pedagogical, since this order better parallels that of projected gradient descent, which may facilitate the understanding of the logic of our algorithm.

Thank you for your positive feedback on our theoretical contribution and for your insightful comments. It seems that the reference items [1-2] in the review are unfortunately missing, so this rebuttal will be based on our understanding of general related literatures. We would really appreciate it and would be happy to discuss in more details if you could provide these references during discussion period.

Distinction from previous work and new insights.

• All of the previous works that used both score-based model and plug-and play, e.g. [RS18, LBA+22, BB23, CDC23] are more or less heuristic and do not have provable consistency or robustness. While some of them realized the potential of split Gibbs sampling as a plug-and-play posterior sampler, none of them found a way to solve the denoising step rigorously. One of our new insights is to show that the denoising step can be solved in full rigor with a reverse diffusion process using only unconditional score function, and it is even possible to make this diffusion process a DDIM one which enables application of many acceleration techniques. Moreover, the gain of our denoising diffusion step is not only mathematical rigor, but also improved robustness, since by running a full, multiple-step reverse diffusion process instead of a one-step heuristic denoiser as in previous works, we are able to better compensate for the drifts from the data manifold induced by the consistency step. It is with this insight that we are able to establish the first provably consistent and robust diffusion posterior sampling algorithm.

Demonstration of convergence.

• To further demonstrate the convergence of our algorithm as shown by our theory, we run a proof-of-concept experiment with data generated from a Gaussian mixture model, so that the convergence of distribution can be examined directly. The detailed settings can be found in General Response to All Reviewers, and the results can be found in the rebuttal pdf. It can be clearly verified from these results that our algorithm does converge to the true posterior distribution.

Additional noise in DDS-DDPM.

• Thank you for your sharp observations. As far as we can see, only the second row in Table 5 has noticeable noise for DPnP-DDPM. However, we also note that in this row, DDS-DDPM also recovers visibly finer details and textures of the board and the text in the original image. This tradeoff between the capability of reconstructing finer details and the risk of introducing additional noise is indeed a general phenomenon that has been observed in previous works, e.g., in [SKZ+23, page 21].

Typos.

• Thank you for your careful reading. We will proofread our paper more carefully and correct this typo among a few others in the final version.

Thank you for your careful review and insightful suggestions. Below we address your concerns in a point-to-point manner.

Comparison with ReSample [SKZ+23].

• Thank you for your suggestion. There are a few recent algorithms that included a stochastic correction step for better robustness. The LGD-MC algorithm [SZY+23], which we compared against in our paper, is also one of them. We chose it due to its simplicity for implementation, but we are also more than happy to include comparison with the non-latent version of ReSample in the final version; comparison of latent versions, on the other hand, requires developing a latent version of our algorithm, which we feel more appropriate to leave to future work as discussed in General Response to All Reviewers. There are some preliminary comparisons in General Response to All Reviewers and in the rebuttal pdf, and we plan to include a more detailed comparison in the final version.

Computational cost.

• Thank you for your suggestion. In our paper, "~3000 NFEs" is the computational cost of the DDPM version of DPnP; we also proposed a DDIM version in the paper, which costs ~1500 NFEs, only 1.5x the cost of DPS (1000 NFEs). Moreover, the DDIM version allows to incorporate acceleration techniques for unconditional DDIM sampling, e.g. DPMSolver++, which can reduce the computation significantly. We would be more than happy to add more discussions in the final version.

Performance scaling with NFEs.

• Thank you for raising this good point. We have included a few samples with different NFEs in the rebuttal, Figure 3, and are willing to investigate it more systematically in the final version if possible.

The name DDS for Denoising Diffusion Sampler.

• Thank you for your suggestion. We will change it to an unambiguous name in the final version.

Notation $G$ in line 136.

• Thank you for your very careful reading and sharp observation. This is a typo for $M$, the forward diffusion process. We will proofread the paper carefully and make corrections in the final version.

Thank you for your detailed review and valuable comments. Below we address your comments in a point-to-point manner.

Examples for Assumption 1 on the forward model and the likelihood (point 1,2).

• A typical scenario is where the measurement noise $\xi$ in the measurement model $y=\mathcal{A}(x^\star)+\xi$ is Gaussian, e.g., $\xi\sim\mathcal N(0,\sigma^2I_m)$. In such case, we have$\mathcal L(x;y)=\log p(y|x^\star=x)=-\frac1{2\sigma^2}\Vert y-\mathcal A(x)\Vert^2 - \frac{m}2\log(2\pi\sigma^2).$ Therefore Assumption 1 holds for all applications where the forward model $\mathcal A$ is almost everywhere differentiable, in particular, if $\mathcal A$ is linear. This includes inpainting, deblurring, super-resolution, phase retrieval, etc. It is noteworthy that most of the existing works (e.g. [CKM+23, SKZ+23, SZY+23, DS24, MSKV24], and [1-4] you listed), have the same differentiability assumption or requires the stronger condition that $\mathcal A$ is linear.

The noise level $\eta$ (point 3).

• The noise level $\eta$ in the DPnP framework is a hyperparameter used for annealing. It is not related to, and should not be confused with, the measurement noise level (the variance of $\xi$ in the measurement model $y=\mathcal A(x^\star)+\xi$). The annealing noise level at the $k$-th iteration, denoted by $\eta_k$ in Algorithm 1, is set manually, hence is always known.

How to know the likelihood term (point 4).

• As the example above (please refer to our response to point 1 and 2) demonstrates, the likelihood can be computed if the forward model $\mathcal{A}$ and the noise distribution is known. This setting is adopted by most of the existing works (please refer to our response to point 1,2), including [1-4] you listed (in fact, most of them assumed $\mathcal{L}(x;y)=-\frac1{2\sigma^2}\Vert y-\mathcal A(x)\Vert^2 + \rm const$, which is a special case of ours with Gaussian measurement noise). It might also be helpful to note that our log-likelihood function $\mathcal L$ is fully, simply determined by the measurement model, and should not be confused with $\log p(y|x_t)$, which was also called likelihood by some of the previous works.

Computation cost of inverting $\Sigma$ (point 5).

• While most of the existing works (including [1-4]) simply assumed (implicitly or explicitly) $\Sigma=\frac1{\sigma^2}I$ is scalar (please refer to our response to point 4), it is possible to handle general $\Sigma$ with standard numerical tricks in our work. For example, one can pre-compute, once and for all, the vector $A^\top\Sigma^{-1}y$ (which amounts to solving a linear equation) and a prefactorizization of the PSD matrix $A^\top\Sigma^{-1}A$. Using such precomputed information, each iteration only needs to compute several matrix-vector products, which has a negligible computation cost compared to Neural Function Evaluations. The cost of computing the prefactorization of $A^\top\Sigma^{-1} A$ has also been observed to be acceptable [SVMK22, KEES22]. Another option is to use MALA for the PCS step as we did for general non-linear inverse problems, so that there is no need to invert any matrix.

Discretization error and TV bound for subsamplers (point 6,7).

• As Theorem 2 indicates, discretization error affects our algorithm only through total variation error of subsamplers DDS and PCS, so we respond to point 6 and 7 together. The reason we did not include an analysis of TV error of subsamplers is, as explained in the paper, that both subsamplers can be analyzed with existing techniques in the references we provided. For example, applying the results in [LWCC23], one may show that DDS-DDPM has total variation error $\varepsilon_{\sf DDS-DDPM}\le \tilde C\frac{d^2}{\sqrt{T'}}+\tilde C\sqrt d\varepsilon_{\sf score},$ where $\tilde C$ hides logarithmic dependence on parameters, $T'$ is the number of steps (cf. Algorithm 2), $\varepsilon_{\sf score}$ is score estimation error defined in [LWCC23]. Similar bound for DDS-DDIM can also be proved using the results therein. On the other hand, utilizing results on the convergence of MALA [CLA+21], one may prove that PCS has total variation error $\varepsilon_{\sf PCS}\le C\exp(-C\sqrt{N}\eta^{-1/6}d^{-1/4}).$

Comparison with works other than DPS (point 8).

• Thank you for the references. As we mentioned in General Response to All Reviewers, we compare our algorithm not only with DPS, but also with more recent, highly competitive algorithm like LGD-MC [SZY+23], and plan to add comparison with ReSample (pre version available in the rebuttal pdf). Note that comparison with such stronger baseline as we have done is actually not too common in the literature, given the highly active status of the field and the well-established position of DPS, e.g. in all the above references and [1-4] you listed, the strongest non-latent baseline is DPS (latent models are discussed in the next paragraph). In addition, none of these previous works come with a provable robustness guarantee, thus our algorithm is also of independent theoretical interest.
• Concerning [1-4], all of them except [2] are confined to linear inverse problems, while our focus is on general, non-linear inverse problems (e.g. phase retrieval and quantized sensing). [1-3] rely on the use of latent diffusion models ([2] further requires prompting the model) which is again a narrower setting than ours as discussed in General Response to All Reviewers. Unavailability of code also hinders comparison with [2-4].

Whether the dataset is the same with DPS (point 9).

• We have strictly followed the original DPS paper on the usage of dataset. We used the same dataset without picking smaller subset than DPS. Our experiments are done for different inverse problems, though, as our focus is on more non-linear problems. In doing so, we made our best effort to ensure fair comparison by fine-tuning all competing algorithms and listing explicitly our problem parameters in the paper.

Dear Reviewer 43oG,

We've taken your initial feedback into careful consideration in our response. Could you kindly confirm whether our responses have appropriately addressed your concerns?

If you find that we have properly addressed your concerns, could you please kindly consider increasing your initial score accordingly? Please let us know if you have further comments.

Thank you for your time and effort in reviewing our work!

Many thanks, Authors