We thank the reviewer for their comments and thoughtful engagement with our work. The concerns expressed in the review are addressed below.

• Algorithm 1. Thanks for catching this, there were indeed some typos in Algorithm 1 (note that the analysis was done on the correct algorithm). The corrected algorithm uses a schedule for noise levels $\sigma_k$ and weights $\alpha_k$. We will also fix the overloading of variable $k$ for the iteration index in the main text which is confusing: we will use $i$ for the PGD algorithm loop and $k$ for the MMSE Averaging loop (it is fixed in the algorithm below). The correct algorithm looks like this:

  $\\hat{x}^{(0)} \\gets y$

  for $i = 0, 1, \\dots$ do

  $\\quad \\hat{x}\_0^{(i+1)} = \\hat{x}^{(i)} - \\tau \\lambda \\nabla f(\\hat{x}^{(i)})$

  $\\quad \\text{for}~k = 0, 1, \\dots, n_{i+1} - 1$ do

  $\\qquad\\sigma\_k \\gets \\sqrt{\\frac{\\tau}{k + 1}}$

  $\\qquad\\alpha_k \\gets \\frac{k+1}{k+2}$

  $\\qquad\\hat{x}\_{k+1}^{(i+1)} \\gets \\alpha_k \\mathrm{MMSE}\_{\\sigma\_k}(\\hat{x}\_k^{(i+1)} ) + (1 - \\alpha\_k) \\hat{x}\_{0}^{(i+1)}$

  $\\quad\\hat{x}^{(i+1)} \\gets \\hat{x}^{(i+1)}\_{n\_{i+1}}$

• Continuous-time version. It is indeed possible to derive a continuous-time analogue of our gradient descent recursion on the smoothed proximal objectives. Specifically, one can write
  \\dot{x}_t = - \\gamma\_t \nabla F\_{\\sigma\_t}(x\_t) = - \gamma\_t (x\_t - y) + \\gamma\_t \\tau \\nabla \\ln p\_{\\sigma\_t}(x\_t) \

  where $\\sigma_t$ and $\\gamma_t$ are time-dependent parameters. Now assuming that $\\gamma_t = (1 + \\tau / \\sigma_t^2)^{-1}$, this can be equivalently rewritten as
  \\dot{x}\_t = - x\_t + \\gamma\_t y + (1 - \\gamma\_t) \\mathrm{MMSE}\_{\\sigma\_t}(x\_t) \

  which would correspond to the continuous-time equivalent of our MMSE averaging sequence. Assuming that $p(x)$ is a Gaussian (hence the MMSE has a closed form), and that $\sigma_t$ decreases as $\tau / t$, we can solve the ODE in closed form and show that it converges to the proximal solution. However, while the ODEs have a similar shape to the one derived in Appendix B of Delbracio & Milanfar (2023), the two algorithms have different objectives and are not entirely comparable: the proximal solution for us, samples from $p(x)$ for InDI. Of course the type of convergence proved in InDi is weaker than what we prove in that it has a) restrictive assumptions requiring Gaussian $p(x)$ and b) is based on a continuous limit which can't be achieved in practice, while ours is based on a discrete recursion.

• Significance. We agree that the assumptions (log-concavity, differentiability) do limit direct applicability. However, as with many theoretical papers, we believe this work serves as a first step toward understanding the behaviour of practical denoising-based schemes from a variational optimisation perspective, especially for a problem where no convergence rates have been proven so far. Besides, even under strong assumptions, establishing convergence to the true proximal operator and MAP estimator is highly non-trivial. We hope this theoretical foundation can guide future extensions to more general priors.

• Lack of experiments. Our primary goal in this paper is to provide theoretical foundations for a class of heuristics already used in practice—such as in Cold Diffusion—which have demonstrated strong empirical performance. Nevertheless, we agree that a thorough empirical study is very valuable, and we believe it would hold as a separate work by itself.

We thank the reviewer for their time and thoughtful feedback. Below we respond to the points raised, highlighting the changes we intend to introduce in the manuscript in response to the comments above.

• Originality of MMSE Averaging. We fully agree that the algorithm is not novel, and we try to clearly state this in the paper, also explicitly linking it to Cold Diffusion (Bansal et al., 2023) and Indi (Delbracio & Milanfar, 2023). We would like to thank the reviewer for pointing out the paper of Feng et al., 2024 which also provides interesting insights connecting iterative denoising to diffusion models. However, our main contribution is to establish rigorous convergence guarantees under specific choices of noise levels and averaging weights of the iterative denoising framework. It can be seen as a theoretical justification for many iterative algorithms which have been proposed with empirical and heuristic bases, but with no rigorous justification.

  While the algorithm itself has been widely used before, the fact that the averaging recursion converges towards the proximal operator is not obvious, and requires noting its equivalence with gradient descent on a sequence of smoothed objectives (Prop. 1 in the paper). Likewise, Algorithm 1 is indeed a PnP method and it highlights one of the potential applications which is enabled by having access to a converging approximation to the proximal operator of $\\log p$ (MMSE Averaging). Contrary to many successful PnP methods which approximate the proximal operator without guarantees on its accuracy (e.g. directly replacing it with a denoiser), we provide a theoretical proof that — under our assumptions — the PnP algorithm we propose converges to the MAP estimate.

• Significance under assumptions. Concerns about the strength of the assumptions are indeed valid. Nevertheless, our goal was to provide a clear and rigorous theoretical foundation in a setting where the analysis is quite delicate, even under strong assumptions. Establishing convergence to the true proximal operator and MAP estimator is highly non-trivial and not common in the literature. We believe this is an important first step towards understanding and extending such results to more realistic, potentially non smooth, settings (see also response to Rev. LqTm on the potential extensions to non-convex settings).

  Note that the constraints on the prior distribution do not directly influence the constraints on the network used for MMSE denoising: while ReLU networks are not smooth, they can be used to approximate smooth functions effectively and hence used in our framework; we will clarify this distinction in the paper as it can indeed be confusing.

  To give an example of a more realistic prior which would satisfy our assumptions, we can start with the TV prior, commonly used for images as it penalizes discontinuities in the data. The TV prior for a discrete 1D signal would be $p(x) \\propto \\exp(-\\sum\_{i=1}^{n-1} |x\_{i+1} - x\_i|)$. Since we need the prior to be differentiable, we can use a smooth approximation to the absolute value $|x| \\approx \\sqrt{x^2+\\delta^2}-\\delta$ for small $\\delta$. Then we can calculate the second and third derivatives to verify that the distribution is indeed log concave, and that the bound on the third derivative depends on $\\delta$. This prior is certainly more realistic than the Gaussian prior for image data and is (or was, before deep learning priors emerged) widely used in inverse problems.

• Theorem 2 convergence to the average. Although last-iterate convergence would be great to have, it typically requires a noticeably more complex analysis. For a good discussion on the difficulties of last-iterate analysis see "Last iterate convergence of SGD for Least-Squares in the Interpolation regime", A. Varre, L. Pillaud-Vivien, N. Flammarion, 2021. In general average-iterate analysis is quite standard in the context of analyzing optimization algorithms, for example see "Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization", M. Schmidt, N. Le Roux, F. Bach, 2011 who prove a similar theorem in a slightly different setting using the average iterate. Note that best-iterate convergence can be obtained easily from average iterate convergence. We will highlight this point in the revised manuscript.

• Typo in $\\alpha\_k$. Thanks for catching this, as you said the correct value is $\\alpha\_k = \\frac{k+1}{k+2}$. Algorithm 1 contained a few inaccuracies, we have restated it in the rebuttal to reviewer 3djj.

• Statement of assumptions. We will modify the theorem statements to be explicit about the required assumptions in the revised manuscript.

We wish to thank the reviewer for their comments and thoughtful engagement with our work. The questions raised in the review in are addressed below.

• Compare to existing inexact proximal algorithms. We of course do not claim novelty in using an inner loop for approximate proximal computations, but we will make sure this is clear in the revised manuscript. As noted by the reviewer, this approach is classical, going back to Schmidt et al. (2011) and earlier cited in Section 4. As discussed in lines 299–302, our analysis differs in that the approximation guarantees for the underlying proximal solver are on the iterates (from Th. 1), not the objective values. This distinction slightly alters the analysis, particularly in controlling error propagation. Nevertheless, we do not wish to claim this as a major contribution but more as a natural consequence of how our approximation of the proximal operator can be integrated into standard PGD frameworks.

• Denoiser approximation error. Indeed in practice we only have access to approximations of the MMSE denoiser. Our analysis can be extended to take into account this additional source of error: if each time the denoiser is applied, the score estimate is off by $\\varepsilon_k$ (i.e. we get $\\nabla \\log p\_{\\sigma\_k}(x\_k) + \\varepsilon\_k$), and if $|\\varepsilon\_k| \\leq \\varepsilon$ is uniformly bounded, then we can prove convergence to a neighborhood of radius $O(\\varepsilon)$ around the true proximal point at a $1/k$ rate. We will include this extension as a comment in the main text, referring to the appendix for a formal statement.

• Dealing with non-convexity. Our theoretical results rely on the convexity of the negative log prior. Handling the general non-convex setting is outside the scope of our work as it would certainly lead to several additional technical hurdles. One possible way forward is to note that convolution with a Gaussian intuitively "convexifies" a probability distribution to some extent. This could significantly help during optimization in practice, and could potentially be useful for theoretical analysis: there exists a minimum $k$ (dependent on the prior $p$) for which $p\_{\\sigma\_k}$ is log-concave. For example, Stochastic Localization via Iterative Posterior Sampling, L. Grenioux et al., 2024 show in a sampling setting that this convexification effect can be used to help sample from non-convex distributions belonging to a restricted class. Another possible direction will be to exploit existing tools allowing to leverage convergence results on convex auxiliary functions, as done by Paquette et al. in Catalyst Acceleration for Gradient-Based Non-Convex Optimization. These directions are perfectly relevant, but beyond the scope of our paper. As noted in the discussion with Rev. wDZG, there are of course limitations to our work, but, as far as we know, this is the first convergence result with known rate obtained for such optimization schemes with non-trivial (non-Gaussian) data distributions.

• Experiments. We agree that the toy example in Section 5 is minimal and serves only illustrative purposes, to convey intuition about the conditioning improvements and the trajectory of the iterates. Our primary goal was indeed to provide convergence results and the connections to imaging arised mainly from noting the similarity between MMSE Averaging and PnP algorithms in imaging settings (as well as notably Cold Diffusion), where they have shown good empirical performance. We will expand the introduction by adding more context about optimization aspects of the problem. In particular, note that the gradient of (Proximal Objective) is not Lipschitz, hence gradient descent converges as $O(\\frac{1}{\\sqrt{k}})$, and acceleration will not improve this rate. Splitting methods on the other hand will require the prox of $\\log p$, which results in a chicken and egg problem. Finally, higher order methods are also intractable since they will require for example the Hessian of $\\log p$ which is not generally known. Hence the rate of $O(\\frac{1}{k})$ obtained in Theorem 1 is a strict improvement for the specific problem we're dealing with.

• Convergence guarantees in PnP. In the introduction and in the related works section, we noted that there exist many PnP methods with convergence guarantees. However, we also noted that they do not converge towards the MAP solution, but rather to the fixed point of some other operator (different from the proximal operator of $-\\log p$). Our goal is to bridge this gap by showing that one can compute the true MAP minimiser via MMSE denoisers. The analysis in Provable Convergence of Plug-and-Play Priors with MMSE Denoisers, Xu et al., 2020, instead follows a different path: the assumptions are much weaker (no convexity required), but convergence is only shown to a stationary point of a function which is not the proximal operator (note the regularizer explicitly stated in Eq. 8 is not equal to $-\\log p(x)$). We will add comparison with this paper to the related work section, since it is indeed relevant.

• Is it necessary to work with $\log p$ (i.e., a log-prior), rather than a general concave function? Yes. Our work requires access to the score function of $-\\log p\_\\sigma$. If instead one starts from a general function $h$, the analogous quantity would be $\\log (\\exp(-h) * \\mathcal{N})$, which is generally intractable as it doesn't have a natural Tweedie-based reformulation. Indeed, as is often the case when dealing with score functions, being able to apply the Tweedie formula is fundamental.

• Non-Gaussian priors. While strong, the assumptions needed for our results to hold do not rule out several interesting non-Gaussian priors. Two examples are

  1. the smoothed TV prior (where the smoothing can be a simple approximation of the absolute value $|x| \\approx \\sqrt{x^2 + \\delta^2} - \\delta$, or the Pseudo-Huber loss) which can be used as a prior for images penalizing discontinuities in their gradient,
  2. a uniform distribution over an interval, convolved with a Gaussian to smooth out the discontinuities at the edges.

  Both satisfy Assumptions 1 and 2, and the third derivative is bounded by a quantity which depends on $\\delta$ for 1. and on the standard deviation of the Gaussian in 2. Unlike in the Gaussian case the 3rd derivative is not uniformly lower than the 2nd on the whole of $\\mathbb{R}$. However, note that this does not make the naive GD on the unsmoothed proximal objective a good algorithm, because we cannot compute $\\nabla \\log p(x)$, but only $\\nabla \\log p\_\\sigma(x)$ for $\\sigma > 0$.

• Comparison to simulated annealing. Both simulated annealing (SA) and the MMSE Averaging algorithm involve decreasing a "temperature" or noise level which decreases over time to ensure convergence. In our case this serves to deal with non-smoothness of the objective, while in the case of SA it helps deal with non-convexity (or even non-continuity issues). While the principle of reducing the noise is the same, the context of the two algorithms remains quite different: MMSE Averaging is based on gradient descent while SA performs Gibbs sampling at each step to decide the next iterate. Furthermore, while $\\sigma\_k$ in MMSE Averaging does not introduce additional noise, and thus the algorithm is deterministic, SA is inherently a stochastic algorithm. We will add a short comment to clarify the connections with SA.

Concerning the minor remarks:

• We agree that "ill-conditioned log-prior" is not standard and we will clarify this. It is taken from the standard use of "ill-conditioned objective function" from the optimization literature.

• We will revise the sentence about “a new wave of approaches” and cite earlier works such as Louchet & Moisan (2013), clarifying that this perspective has an older history which has recently gained renewed attention. Thanks for pointing this out.

• We will restate all theoretical results with full hypotheses for clarity, and similarly recall the assumptions behind the Beck and Teboulle (2009) result.

• Regarding the comparison between $O(\\log(1/\\varepsilon))$ and $O(1/\\varepsilon)$: our point concerns the impact of conditioning. We will clarify this in the discussion.

• We will fix missing notations like $\\gamma_k$ and J^\\star, and rewrite unclear sentences accordingly.

Thank you for your additional comments and for your positive feedback. Below are some clarifications in response to the points you raised:

On the approximate denoiser setting: As is often the case when writing a paper, some natural extensions—while conceptually straightforward—are not initially considered or prioritised. We fully agree that including a convergence result in the presence of denoiser approximation error is relevant. We are happy to provide this result in the revised version and appreciate your suggestion.

On non-Gaussian priors: After re-reading the corresponding paragraph in our rebuttal, we recognise that it was quite unclear—we apologise for the confusion, and would like to clarify our point here.

The two examples we gave (a smoothed TV prior and the convolution of a uniform distribution with a Gaussian) were meant to illustrate non-Gaussian priors that nonetheless satisfy our assumptions. In particular, the convolution of a uniform distribution with a Gaussian is indeed log-concave due to the preservation of log-concavity under convolution.

However, the final sentence of that paragraph— “note that this does not make the naive GD on the unsmoothed proximal objective a good algorithm, because we cannot compute $\nabla \ln p$, but only $\nabla \ln p_\sigma$ for $\sigma > 0$”—was poorly phrased and misleading. What we intended to convey is the following:

• While for some toy examples (like the uniform or smoothed TV prior) $\nabla \log p$ may be tractable, in practical settings involving learned denoisers, we do not have access to $\nabla \ln p$, but only to the smoothed scores $\nabla \ln p_\sigma$ for $\sigma > 0$. Therefore, naive gradient descent on the original objective (involving $\log p$) cannot be directly implemented, whereas our proposed method remains applicable.

We also wish to emphasise that the primary motivation of our work is to address priors that do not have a closed-form expression and are instead learned from data, as is typically the case in modern applications.

We hope this better explains our intent. Please don’t hesitate to reach out if this remains unclear or raises further questions.

We wish to thank the reviewer for their comments and thoughtful engagement with our work. The questions raised in the review in are addressed below.

• Intractable $\\log p$. We agree and would like to clarify that, in practice, we do not compute $\\log p\_\\sigma$ directly. Instead, modern denoising models (e.g., trained MMSE denoisers, score-based diffusion models, flow-based models) are explicitly trained to approximate the score $\\nabla \\log p\_\\sigma$ for a range of noise levels. These models can be trained on randomly sampled noise levels, which enables them to interpolate effectively for intermediate values, including those required by our method. Hence, our analysis applies to practical implementations where this score is approximated via learned denoisers.

• Comparison with using $\\nabla \\log p(x, \\sigma = \\mathrm{small})$ directly. Using a fixed small noise level to approximate the score of the original prior is a valid approach, but as discussed in the paper, it requires performing very small step sizes that can lead to arbitrarily slow convergence to the MAP of a smoothed density. In contrast, the scheme we analyze leverages large $\\sigma$ in early iterations and gradually anneals to smaller $\\sigma$ to reach the true proximal point with $O(1/k)$ rate (cf. Theorem 1 and paragraph right after, lines 202-218). Importantly, this scheme does not rely on any smoothness parameter about $\\log p$ (e.g., Lipschitz constants), which are typically unknown in practice. There are of course limitations due to the assumptions made about $\\log p$, but, as far as we know, this is the first convergence result with known rate obtained for such optimization schemes with non-trivial (non-Gaussian) data distributions.

• Sensitivity to the noise schedule $\\sigma\_k$ in Theorems 1/2. Our analysis can indeed easily be adapted to a more general class of decreasing noise schedules. In particular, for sequences $\\sigma\_k$ such that $\\sum\_k \\sigma\_k = \\infty$ and corresponding step sizes $\\gamma_k = 1 / (1 + \\tau / \\sigma_k^2)$, we obtain a general convergence rate of the form: $O \\Big( \\min(\\sigma\_k, \\exp( - \\sum\_{j = 1}^k \\sigma_j / \\tau ) \\Big)$ for the smoothed gradient descent recursion of Proposition 1. Thus, any $\\sigma\_k$ decreasing at most as $\\tau / k$ ensures a $O(1/k)$ rate. The specific choice $\\sigma_k^2 = \\tau / (k+1)$ is made to emphasize the connection with the MMSE Averaging scheme $x\_{k+1} = \\alpha\_k \mathrm{MMSE}\_{\\sigma\_k}(x\_k) + (1 - \\alpha\_k) y$, while other choices of noise schedules and averaging weights will lose this equivalent formulation. We will clarify this in the paper.

• Comparison to preconditioned GD when $\\log p$ is tractable. If $\\log p$ and its gradient were tractable, it is very likely that better optimization schemes could be used, but it would fall beyond the typical use case where $p$ represents an unknown data distribution. That being said, theoretical convergence with a preconditioner $P$ would depend on the condition number $\\kappa_P = \\sup\_x \\lambda\_{\\mathrm{max}}(P \\nabla^2 \\log p(x)) / \\inf\_x \ \\lambda\_{\\mathrm{min}}(P \\nabla^2 \\log p(x))$ with a rate of $O(\\kappa\_P \\log(1/\\varepsilon))$. In the best cases this can lead to a considerable speed-up with $\\kappa_P$ close to $1$ but in practice for very badly conditioned functions $\\log p$ the matrix inversion (the optimal $P$ being $(\\nabla^2 \\log p(x))^{-1}$) could lead to severe numerical errors. In contrast, the MMSE averaging scheme does not require any knowledge of $\\log p$ and does not explicitly require any sort of matrix inversion, as long as the MMSE estimator can be computed/well approximated.

• Heuristics for the number of iterations in Algorithm 1. A full empirical study is beyond the scope of our paper, but it is indeed an important future direction to investigate. Existing works such as Delbracio & Milanfar (2023) suggest that 5–20 inner denoising steps are often sufficient for high-quality results. We expect our method to behave similarly: the inner loop (MMSE Averaging) benefits from rapid convergence due to improved conditioning in early iterations. We will add a comment to the paper with guidance and references.

• $\\alpha\_k$ Correction. Thank you for catching this. Indeed, there is a typo: the correct choice used throughout the theoretical analysis is $\alpha_k = \\frac{k+1}{k+2}$, not $\frac{1}{k+2}$. We will correct this to be consistent in the paper.