SNIPS: Solving Noisy Inverse Problems Stochastically
Year: 2021

Abstract:
In this work we introduce a novel stochastic algorithm dubbed SNIPS, which draws samples from the posterior distribution of any linear inverse problem, where the observation is assumed to be contaminated by additive white Gaussian noise. Our solution incorporates ideas from Langevin dynamics and Newton’s method, and exploits a pre-trained minimum mean squared error (MMSE) Gaussian denoiser. The proposed approach relies on an intricate derivation of the posterior score function that includes a singular value decomposition (SVD) of the degradation operator, in order to obtain a tractable iterative algorithm for the desired sampling. Due to its stochasticity, the algorithm can produce multiple high perceptual quality samples for the same noisy observation. We demonstrate the abilities of the proposed paradigm for image deblurring, super-resolution, and compressive sensing. We show that the samples produced are sharp, detailed and consistent with the given measurements, and their diversity exposes the inherent uncertainty in the inverse problem being solved.

heading: 1 Introduction
text: Many problems in the field of image processing can be cast as noisy linear inverse problems. This family of tasks includes denoising, inpainting, deblurring, super resolution, compressive sensing, and many other image recovery problems. A general linear inverse problem is posed as
y = Hx + z, (1)
where we aim to recover a signal x from its measurement y, given through a linear degradation operator H and a contaminating noise, being additive, white and Gaussian, z ∼ N ( 0, σ20I ) . In this work we assume that both H and σ0 are known.
Over the years, many strategies, algorithms and underlying statistical models were developed for handling image restoration problems. A key ingredient in many of the classic attempts is the prior that aims to regularize the inversion process and lead to visually pleasing results. Among the various options explored, we mention sparsity-inspired techniques [13, 55, 11], local Gaussian-mixture modeling [57, 63], and methods relying on non-local self-similarity [6, 9, 36, 51]. More recently, and with the emergence of deep learning techniques, a direct design of the recovery path from y to an estimate of x took the lead, yielding state-of-the-art results in various linear inverse problems, such as denoising [25, 59, 61, 52], deblurring [22, 48], super resolution [10, 17, 54] and other tasks [29, 28, 19, 16, 37, 58].
Despite the evident success of the above techniques, many image restoration algorithms still have a critical shortcoming: In cases of severe degradation, most recovery algorithms tend to produce washed out reconstructions that lack details. Indeed, most image restoration techniques seek a reconstruction that minimizes the mean squared error between the restored image, x̂, and the unknown original one, x. When the degradation is acute and information is irreversibly lost, image reconstruction becomes a highly ill-posed problem, implying that many possible clean images could explain the given
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
measurements. The MMSE solution averages all these candidate solutions, being the conditional mean of the posterior of x given y, leading to an image with loss of fine details in the majority of practical cases. A recent work reported in [5] has shown that reconstruction algorithms necessarily suffer from a perception-distortion tradeoff, i.e., targeting a minimization of the error between x̂ and x (in any metric) is necessarily accompanied by a compromised perceptual quality. As a consequence, as long as we stick to the tendency to design recovery algorithms that aim for minimum MSE (or other distances), only a limited perceptual improvement can be expected.
When perceptual quality becomes our prime objective, the strategy for solving inverse problems must necessarily change. More specifically, the solution should concentrate on producing a sample (or many of them) from the posterior distribution p (x|y) instead of its conditional mean. Recently, two such approaches have been suggested – GAN-based and Langevin sampling. Generative Adversarial Networks (GANs) have shown impressive results in generating realistically looking images (e.g., [14, 35]). GANs can be utilized for solving inverse problems while producing high-quality images (see e.g. [2, 31, 34]). These solvers aim to produce a diverse set of output images that are consistent with the measurements, while also being aligned with the distribution of clean examples. A major disadvantage of GAN-based algorithms for inverse problems is their tendency (as practiced in [2, 31, 34]) to assume noiseless measurements, a condition seldom met in practice. An exception to this is the work reported in [33], which adapts a conditional GAN to become a stochastic denoiser.
The second approach for sampling from the posterior, and the one we shall be focusing on in this paper, is based on Langevin dynamics. This core iterative technique enables sampling from a given distribution by leveraging the availability of the score function – the gradient of the log of the probability density function [38, 3]. The work reported in [44, 20, 46] utilizes the annealed Langevin dynamics method, both for image synthesis and for solving noiseless inverse problems.1 Their synthesis algorithm relies on an MMSE Gaussian denoiser (given as a neural network) for approximating a gradually blurred score function. In their treatment of inverse problems, the conditional score remains tractable and manageable due to the noiseless measurements assumption.
The question addressed in this paper is the following: How can the above line of Langevin-based work be generalized for handling linear inverse problems, as in Equation 1, in which the measurements are noisy? A partial and limited answer to this question has already been given in [21] for the tasks of image denoising and inpainting. The present work generalizes these ([44, 20, 46, 21]) results, and introduces a systematic way for sampling from the posterior distribution of any given noisy linear inverse problem. As we carefully show, this extension is far from being trivial, due to two prime reasons: (i) The involvement of the degradation operator H, which poses a difficulty for establishing a relationship between the reconstructed image and the noisy observation; and (ii) The intricate connection between the measurements’ and the synthetic annealed Langevin noise. Our proposed remedy is a decorrelation of the measurements equation via a singular value decomposition (SVD) of the operator H, which decouples the dependencies between the measurements, enabling each to be addressed by an adapted iterative process. In addition, we define the annealing noise to be built as portions of the measurement noise itself, in a manner that facilitates a constructive derivation of the conditional score function.
Following earlier work [44, 20, 46, 21], our algorithm is initialized with a random noise image, gradually converging to the reconstructed result, while following the direction of the log-posterior gradient, estimated using an MMSE denoiser. Via a careful construction of the gradual annealing noise sequence, from very high values to low ones, the entries in the derived score switch mode. Those referring to non-zero singular values start by being purely dependent on the measurements, and then transition to incorporate prior information based on the denoiser. As for entries referring to zero singular values, their corresponding entries undergo a pure synthesis process based on the prior-only score function. Note that the denoiser blends values in the evolving sample, thus intermixing the influence of the gradient entries. Our derivations include an analytical expression for a positiondependent step size vector, drawing inspiration from Newton’s method in optimization. This stabilizes the algorithm and is shown to be essential for its success.
We refer hereafter to our algorithm as SNIPS (Solution of Noisy Inverse Problems Stochastically). Observe that as we target to sample from the posterior distribution p (x|y), different runs of SNIPS on the same input necessarily yield different results, all of which valid solutions to the given inverse
1The work reported in [18, 42, 26] and [15, 23] is also relevant to this discussion, but somewhat different. We shall specifically address these papers’ content and its relation to our work in section 2.
problem. This should not come as a surprise, as ill-posedness implies that there are multiple viable solutions for the same data, as has already been suggested in the context of super resolution [31, 2, 34]. We demonstrate SNIPS on image deblurring, single image super resolution, and compressive sensing, all of which contain non-negligible noise, and emphasize the high perceptual quality of the results, their diversity, and their relation to the MMSE estimate.
To summarize, this paper’s contributions are threefold:
• We present an intricate derivation of the blurred posterior score function for general noisy inverse problems, where both the measurement and the target image contain delicately inter-connected additive white Gaussian noise.
• We introduce a novel stochastic algorithm – SNIPS – that can sample from the posterior distribution of these problems. The algorithm relies on the availability of an MMSE denoiser.
• We demonstrate impressive results of SNIPS on image deblurring, single image super resolution, and compressive sensing, all of which are highly noisy and ill-posed.
Before diving into the details of this work, we should mention that using Gaussian denoisers iteratively for handling general linear inverse problems has been already proposed in the context of the Plugand-Play-Prior (PnP) method [53] and RED [39], and their many followup papers (e.g., [60, 30, 1, 49, 7, 50, 40, 4]). However, both PnP and RED are quite different from our work, as they do not target sampling from the posterior, but rather focus on MAP or MMSE estimation.

heading: 2 Background
text: The Langevin dynamics algorithm [3, 38] suggests sampling from a probability distribution p (x) using the iterative transition rule
xt+1 = xt + α∇xt log p (xt) + √ 2αzt , (2)
where zt ∼ N (0, I) and α is an appropriately chosen small constant. The added zt allows for stochastic sampling, avoiding a collapse to a maximum of the distribution. Initialized randomly, after a sufficiently large number of iterations, and under some mild conditions, this process converges to a sample from the desired distribution p (x) [38].
The work reported in [44] extends the aforementioned algorithm into annealed Langevin dynamics. The annealing proposed replaces the score function in Equation 2 with a blurred version of it, ∇x̃t log p (x̃t), where x̃t = xt + n and n ∼ N ( 0, σ2I ) is a synthetically injected noise. The core idea is to start with a very high noise level σ and gradually drop it to near-zero, all while using a step size α dependent on the noise level. These changes allow the algorithm to converge much faster and perform better, because it widens the basin of attraction of the sampling process. The work in [20] further develops this line of work by leveraging a brilliant relation attributed to Miyasawa [32] (also known as Stein’s integration by parts trick [47] or Tweedie’s identity [12]). It is given as
∇x̃t log p (x̃t) = D (x̃t, σ)− x̃t
σ2 , (3)
where D (x̃t, σ) = E [x|x̃t] is the minimizer of the MSE measure E [ ‖x−D (x̃t, σ) ‖22 ] , which can be approximated using a denoising neural network. This facilitates the use of denoisers in Langevin dynamics as a replacement for the evasive score function.
When turning to solve inverse problems, previous work suggests sampling from the posterior distribution p (x|y) using annealed Langevin dynamics [20, 46, 21] or similar methods [15, 18, 42, 26], by replacing the score function used in the generation algorithm with a conditional one. As it turns out, if limiting assumptions can be posed on the measurements formation, the conditional score is tractable, and thus generalization of the annealed Langevin process to these problems is within reach. Indeed, in [44, 20, 46, 42, 26] the core assumption is y = Hx for specific and simplified choices of H and with no noise in the measurements. The works in [15, 23] avoid these difficulties altogether by returning to the original (non-annealed) Langevin method, with the unavoidable cost of becoming extremely slow. In addition, their algorithms are demonstrated on inverse problems in which the additive noise is restricted to be very weak. The work in [21] is broader, allowing for an arbitrary additive white Gaussian noise, but limits H to the problems of denoising or inpainting. While all these works demonstrate high quality results, there is currently no clear way for deriving the blurred score function of a general linear inverse problem as posed in Equation 1. In the following, we present such a derivation.

heading: 3 The Proposed Approach: Deriving the Conditional Score Function

heading: 3.1 Problem Setting
text: We consider the problem of recovering a signal x ∈ RN (where x ∼ p (x) and p (x) is unknown) from the observation y = Hx + z, where y ∈ RM ,H ∈ RM×N ,M ≤ N, z ∼ N ( 0, σ20I ) , and H and σ0 are known.2 Our ultimate goal is to sample from the posterior p (x|y). However, since access to the score function∇x log p(x|y) is not available, we retarget our goal, as explained above, to sampling from blurred posterior distributions, p (x̃|y), where x̃ = x + n and n ∼ N ( 0, σ2I ) , with noise levels σ starting very high, and decreasing towards near-zero.
As explained in the supplemental material, the sampling should be performed in the SVD domain in order to get a tractable derivation of the blurred score function. Thus, we consider the singular value decomposition (SVD) of H, given as H = UΣVT , where U ∈ RM×M and V ∈ RN×N are orthogonal matrices, and Σ ∈ RM×N is a rectangular diagonal matrix containing the singular values of H, denoted as {sj}Mj=1 in descending order (s1 > s2 > · · · > sM−1 > sM ≥ 0). For convenience of notations, we also define sj = 0 for j =M + 1, . . . , N . To that end, we notice that
p (x̃|y) = p ( x̃|UTy ) = p ( VT x̃|UTy ) . (4)
The first equality holds because the multiplication of y by the orthogonal matrix UT does not add or remove information, and the second equality holds because the multiplication of x̃ by VT does not change its probability [24]. Therefore, sampling from p ( VT x̃|UTy ) and then multiplying the result by V will produce the desired sample from p (x̃|y). As we are using Langevin dynamics, we need to calculate the conditional score function∇VT x̃ log p ( VT x̃|UTy ) . For simplicity, we denote hereafter yT = UTy, zT = UT z, xT = VTx, nT = ΣVTn, and x̃T = VT x̃. Observe that with these notations, the measurements equation becomes
y = Hx + z = UΣVTx + z,
and thus
UTy = ΣVTx + UT z = ΣVT (x̃− n) + UT z = ΣVT x̃−ΣVTn + UT z, where we have relied on the relation x̃ = x + n. This leads to
yT = Σx̃T − nT + zT . (5)
In this formulation, which will aid in deriving the conditional score, our aim is to make design choices on nT such that zT − nT has uncorrelated entries and is independent of x̃T . This brings us to the formation of the synthetic annealed noise, which is an intricate ingredient in our derivations.
2We assume M ≤ N for ease of notations, and because this is the common case. However, the proposed approach and all our derivations work just as well for M > N .
We base this formation on the definition of a sequence of noise levels {σi}L+1i=1 such that σ1 > σ2 > · · · > σL > σL+1 = 0, where σ1 is high (possibly σ1 > ‖x‖∞) and σL is close to zero. We require that for every j such that sj 6= 0, there exists ij such that σijsj < σ0 and σij−1sj > σ0. This implies ∀i : σisj 6= σ0, which helps ease notations. SNIPS works just as well for σisj = σ0.
Using {σi}L+1i=1 , we would like to define {x̃i} L+1 i=1 , a sequence of noisy versions of x, where the noise level in x̃i is σi. One might be tempted to define these noise additions as independent of the measurement noise z. However, this option leads to a conditional score term that cannot be calculated analytically, as explained in the supplemental material. Therefore, we define these noise additions differently, as carved from z in a gradual fashion. To that end, we define x̃L+1 = x, and for every i = L,L− 1, . . . , 1: x̃i = x̃i+1 + ηi, where ηi ∼ N ( 0, ( σ2i − σ2i+1 ) I ) . This results in
x̃i = x + ni, where ni = ∑L k=i ηk ∼ N ( 0, σ2i I ) .
And now we turn to define the statistical dependencies between the measurements’ noise z and the artificial noise vectors ηi. Since ηi and z are each Gaussian with uncorrelated entries, so are the components of the vectors ΣVTηi, ΣV
Tni, and zT . In order to proceed while easing notations, let us focus on a single entry j in these three vectors, for which sj > 0, and omit this index. We denote these entries as ηT,i, nT,i and zT , respectively. We construct ηT,i such that
E [ηT,i · zT ] =  E [ η2T,i ] for i ≥ ij E [( zT − nT,ij )2] for i = ij − 1
0 otherwise.
This implies that the layers of noise ηT,L+1, . . . ηT,ij are all portions of zT itself, with an additional portion being contained in ηT,ij−1. Afterwards, ηT,i become independent of zT . In the case of sj = 0, the above relations simplify to be E[ηT,i · zT ] = 0 for all i, implying no statistical dependency between the given and the synthetic noises. Consequently, it can be shown that the overall noise in Equation 5 satisfies
( ΣVTni − zT ) j = nT,i − zT ∼
{ N ( 0, s2jσ 2 i − σ20 ) if σisj > σ0
N ( 0, σ20 − s2jσ2i ) otherwise.
(6)
The top option refers to high values of the annealed Langevin noise, in which, despite the possible decay caused by the singular value sj , this noise is stronger than zT . In this case, nT,i contains all zT and an additional independent portion of noise. The bottom part assumes that the annealed noise (with the influence of sj) is weaker than the measurements’ noise, and then it is fully immersed within zT , with the difference being Gaussian and independent.

heading: 3.2 Derivation of the Conditional Score Function
text: The above derivations show that the noise in Equation 5 is zero-mean, Gaussian with uncorrelated entries and of known variance, and this noise is independent of x̃i. Thus Equation 5 can be used conveniently for deriving the measurements part of the conditional score function. We denote x̃T = V
T x̃i, x̃ = x̃i, n = ni for simplicity, and turn to calculate ∇x̃T log p (x̃T |yT ). We split x̃T into three parts: (i) x̃T,0 refers to the entries j for which sj = 0; (ii) x̃T,< corresponds to the entries j for which 0 < σisj < σ0; and (iii) x̃T,> includes the entries j for which σisj > σ0. Observe that this partition of the entries of x̃T is non-overlapping and fully covering. Similarly, we partition every vector v ∈ RN into v0,v<,v>, which are the entries of v corresponding to x̃T,0, x̃T,<, x̃T,>, respectively. Furthermore, we define v 60,v 6<,v6> as all the entries of v except v0,v<,v>, respectively. With these definitions in place, the complete derivation of the score function is detailed in the supplemental material, and here we bring the final outcome. For x̃T,0, the score is independent of the measurements and given by
∇x̃T,0 log p (x̃T |yT ) = ( VT∇x̃ log p (x̃) ) 0 . (7)
For the case of x̃T,>, the expression obtained is only measurements-dependent, ∇x̃T,> log p (x̃T |yT ) = ( ΣT ( σ2iΣΣ T − σ20I )† (yT −Σx̃T ) ) > . (8)
Lastly, for the case of x̃T,<, the conditional score includes two terms – one referring to the plain (blurred) score, and the other depending on the measurements,
∇x̃T,< log p (x̃T |yT ) = ( ΣT ( σ20I− σ2iΣΣ T )† (yT −Σx̃T ) ) < + ( VT∇x̃ log p (x̃) ) < . (9)
As already mentioned, the full derivations of equations 7, 8, and 9 are detailed in the supplemental material. Aggregating all these results together, we obtain the following conditional score function:
∇x̃T log p (x̃T |yT ) = ΣT ∣∣∣σ20I− σ2iΣΣT ∣∣∣† (yT −Σx̃T ) + (VT∇x̃ log p (x̃))∣∣ 6> , (10)
where (v)| 6> is the vector v, but with zeros in its entries that correspond to v>. Observe that the first term in Equation 10 contains zeros in the entries corresponding to x̃T,0, matching the above calculations. The vector ∇x̃ log p (x̃) can be estimated using a neural network as in [44], or using a pre-trained MMSE denoiser as in [20, 21]. All the other elements of this vector are given or can be easily obtained from H by calculating its SVD decomposition once at the beginning.

heading: 4 The Proposed Algorithm
text: Armed with the conditional score function in Equation 10, the Langevin dynamics algorithm can be run with a constant step size or an annealed step size as in [44], and this should converge to a sample from p (x̃T |yT ). However, for this to perform well, one should use a very small step size, implying a devastatingly slow convergence behavior. This is mainly due to the fact that different entries of x̃T advance at different speeds, in accord with their corresponding singular values. As the added noise in each step has the same variance in every entry, this leads to an unbalanced signal-to-noise ratio, which considerably slows down the algorithm.
In order to mitigate this problem, we suggest using a step size vector αi ∈ RN . We denote Ai = diag (αi), and obtain the following update formula for a Langevin dynamics algorithm:
VT x̃i = V T x̃i−1 + c ·Ai · ∇VT x̃i log p ( VT x̃i|yT ) + √ 2 · cA 1 2 i · zi, (11)
where the conditional score function is estimated as described in subsection 3.2, and c is some constant. For the choice of the step sizes in the diagonal of Ai, we draw inspiration from Newton’s method in optimization, which is designed to speed up convergence to local maximum points. The update formula in Newton’s method is the same as Equation 11, but without the additional noise zi, and with Ai being the negative inverse Hessian of log p ( VT x̃i|yT ) . We calculate a diagonal approximation of the Hessian, and set Ai to be its negative inverse. We also estimate the conditional score function using Equation 10 and a neural network. Note that this mixture of Langevin dynamics and Newton’s method has been suggested in a slightly different context in [43], where the Hessian was approximated using a Quasi-Newton method. In our case, we analytically calculate a diagonal approximation of the negative inverse Hessian and obtain the following:
(αi)j =  σ2i , sj = 0 σ2i − σ20 s2j , σisj > σ0 σ2i · ( 1− s2j σ2i σ20 ) , 0 < σisj < σ0.
(12)
The full derivations for each of the three cases are detailed in the supplemental material. Using these step sizes, the update formula in Equation 11, the conditional score function in Equation 10, and a neural network s (x̃, σ) that estimates the score function ∇x̃ log p (x̃),3 we obtain a tractable iterative algorithm for sampling from p (x̃L | y), where the noise in x̃L is sufficiently negligible to be considered as a sampling from the ideal image manifold.
Algorithm 1: SNIPS
Input: {σi}Li=1, c, τ , y, H, σ0 1 U,Σ,V← svd (H) 2 Initialize x0 with random noise U [0, 1] 3 for i← 1 to L do 4 (Ai)0 ← σ2i I 5 (Ai)< ← σ2i · ( I− σ 2 i
σ20 Σ<Σ<
T )
6 (Ai)> ← σ2i I− σ20Σ † >Σ †T > 7 for t← 1 to τ do 8 Draw zt ∼ N (0, I)
9 dt ← ΣT · ∣∣∣σ20I− σ2iΣΣT ∣∣∣† · (UTy −ΣVTxt−1)+ (VT · s (xt−1, σi))∣∣ 6>
10 xt ← V · ( VTxt−1 + cAidt + √ 2cA 1 2 i zt ) 11 end 12 x0 ← xτ 13 end
Output: x0
Note that when we set H = 0 and σ0 = 0, implying no measurements, the above algorithm degenerates to an image synthesis, exactly as in [44]. Two other special cases of this algorithm are obtained for H = I or H = I with some rows removed, the first referring to denoising and the second to noisy inpainting, both cases shown in [21]. Lastly, for the choices of H as in [20] or [44, 46] and with σ0 = 0, the above algorithm collapses to a close variant of their proposed iterative methods.

heading: 5 Experimental Results
text: In our experiments we use the NCSNv2 [45] network in order to estimate the score function of the prior distribution. Three different NCSNv2 models are used, each trained separately on the training sets of: (i) images of size 64× 64 pixels from the CelebA dataset [27]; (ii) images of size 128× 128 pixels from LSUN [56] bedrooms dataset; and (iii) LSUN 128× 128 images of towers. We demonstrate SNIPS’ capabilities on the respective test sets for image deblurring, super resolution, and compressive sensing. In each of the experiments, we run our algorithm 8 times, producing 8 samples for each input. We examine both the samples themselves and their mean, which serves as an approximation of the MMSE solution, E [x|y].
3Recall that s (x̃, σ) = (D (x̃, σ)− x̃) /σ2, being a denoising residual.
For image deblurring, we use a uniform 5 × 5 blur kernel, and an additive white Gaussian noise with σ0 = 0.1 (referring to pixel values in the range [0, 1]). Figure 1 demonstrates the obtained results for several images taken from the CelebA dataset. As can be seen, SNIPS produces visually pleasing, diverse samples.
For super resolution, the images are downscaled using a block averaging filter, i.e., each nonoverlapping block of pixels in the original image is averaged into one pixel in the low-resolution image. We use blocks of size 2× 2 or 4× 4 pixels, and assume the low-resolution image to include an additive white Gaussian noise. We showcase results on LSUN and CelebA in Figures 2, 4, and 5.
For compressive sensing, we use three random projection matrices with singular values of 1, that compress the image by 25%, 12.5%, and 6.25%. As can be seen in Figure 3 and as expected, the more aggressive the compression, the more significant are the variations in reconstruction.
We calculate the average PSNR (peak signal-to-noise ratio) of each of the 8 samples in our experiments, as well as the PSNR of their mean, as shown in Table 1. In all the experiments, the empirical conditional mean presents an improvement of around 2.4 dB in PSNR, even though it is less visually appealing compared to the samples. This is consistent with the theory in [5], which states that the difference in PSNR between posterior samples and the conditional mean (the MMSE estimator) should be 3 dB, with the MMSE estimator having poorer perceptual quality but better PSNR.
A comparison of our deblurring results to those obtained by RED [39] is detailed in the supplemental material. We show that SNIPS exhibits superior performance over RED, achieving more than 11% improvement in PSNR and more than 58% improvement in LPIPS [62], a perceptual quality metric.

heading: 5.1 Assessing Faithfulness to the Measurements
text: A valid solution to an inverse problem should satisfy two conditions: (i) It should be visually pleasing, consistent with the underlying prior distribution of images, and (ii) It should be faithful to the given measurement, maintaining the relationship as given in the problem setting. Since the prior distribution
is unknown, we assess the first condition by visually observing the obtained solutions and their tendency to look realistic. As for the second condition, we perform the following computation: We degrade the obtained reconstruction x̂ by H, and calculate its difference from the given measurement y, obtaining y −Hx̂. According to the problem setting, this difference should be an additive white Gaussian noise vector with a standard deviation of σ0. We examine this difference by calculating its empirical standard deviation, and performing the Pearson-D’Agostino [8] test of normality on it, accepting it as a Gaussian vector if the obtained p-value is greater than 0.05. We also calculate the Pearson correlation coefficient (denoted as ρ) among neighboring entries, accepting them as uncorrelated for coefficients smaller than 0.1 in absolute value. In all of our tests, the standard deviation matches σ0 almost exactly, the Pearson correlation coefficient satisfies |ρ| < 0.1, and we obtain p-values greater than 0.05 in around 95% of the samples (across all experiments). These results empirically show that our algorithm produces valid solutions to the given inverse problems.

heading: 6 Conclusion and Future Work
text: SNIPS, presented in this paper, is a novel stochastic algorithm for solving general noisy linear inverse problems. This method is based on annealed Langevin dynamics and Newton’s method, and relies on the availability of a pre-trained Gaussian MMSE denoiser. SNIPS produces a random variety of high quality samples from the posterior distribution of the unknown given the measurements, while guaranteeing their validity with respect to the given data. This algorithm’s derivation includes an intricate choice of the injected annealed noise in the Langevin update equations, and an SVD decomposition of the degradation operator for decoupling the measurements’ dependencies. We demonstrate SNIPS’ success on image deblurring, super resolution, and compressive sensing.
Extensions of this work should focus on SNIPS’ limitations: (i) The need to deploy SVD decomposition of the degradation matrix requires a considerable amount of memory and computations, and hinders the algorithm’s scalability; (ii) The current version of SNIPS does not handle general content images, a fact that is related to the properties of the denoiser being used [41]; and (iii) SNIPS, as any other Langevin based method, requires (too) many iterations (e.g., in our super-resolution tests on CelebA, 2 minutes are required for producing 8 sample images), and means for its acceleration should be explored.

heading: 7 Funding Transparency Statement
text: This research was partially supported by the Israel Science Foundation (ISF) under Grant 335/18 and the Technion Hiroshi Fujiwara Cyber Security Research Center and the Israel Cyber Bureau. Bahjat Kawar’s scholarship was partially provided by Li Ka Shing Fellowships and the Planning and Budgeting Committee of the Israel Council for Higher Education.

Review summaries:
This manuscript describes a novel method for sampling from the posterior in a linear inverse problem subjected to noise for an arbitrary measurement matrix. This is achieved by combining a black-box prior (such as that derived from a minimum MSE Gaussian denoiser) whose score function is easily computable with an annealed Langevin dynamics framework. The novelty is the construction of the Langevin noise term, which guarantees that the conditional score function is tractable. In addition, the authors propose a heuristic for step size calculation inspired by Newton's method in optimization. The method is evaluated on a variety of linear inverse problems: image deblurring, super resolution, and compressive sensing. For each task, realistic results are achieved and the residual is shown to conform to the imposed statistical model. 

This paper proposes an algorithm called SNIPS that solves noisy linear inverse problems via an approximate MMSE estimator which makes use of pre-trained Gaussian denoisers. The estimator is based on MCMC sampling which is driven by a synthetically annealed version of Langevin dynamics. The motivations and background leading up to the algorithms derivation are well discussed. Experiments show interesting results with regard to the diversity and faithfulness of the images generated.   ---Edit after author responses---  Having read the other reviews and the responses posted by the authors, I am reducing my original scores from 6 --> 4.  My main issues are that the submission lacks (i) a comparison with an experimental baseline, and (ii) clarity in explaining the technical sections of the paper.

The authors address the problem of conditional sampling using Langevin dynamics. It is known that annealed Langevin dynamics speeds up convergence, and prior work that uses Langevin dynamics for conditional sampling does not correctly model the effects of annealing. Specifically, since there is some annealing noise added to the estimate, the likelihood distribution of the measurements given an estimate will change, and the authors derive approximations to the correct functional form. Experiments show that their new estimator can be implemented in practice. 

This paper proposes a technique based on annealed Langevin dynamics to solve noisy linear inverse problems. The technique is built on top of [18, 19] and proposes to perform the Langevin dynamics in the spectral space of the linear degradation operator $H$. To improve the performance of the proposed method, the authors use coordinate wise step sizes obtained by approximating the hessian with its diagonal. Qualitative results from the proposed method are presented for compress sensing, inpainting, super-résolution and deblurring.

The paper describes an apparently new approach to solving noisy inverse problems.  It shows how the approach can be applied to three different noisy inverse problems: which it calls image deblurring, super-resolution, and compressive sensing.   It validates the approach on synthetic data which is synthesized (eg blurred and noise added) starting from publicly available imagery.  It is carefully written and I find it readable.