SNIPS: Solving Noisy Inverse Problems Stochastically
Authors: Bahjat Kawar
Gregory Vaksman
Michael Elad
Venue: NeurIPS 2021 Poster (2021)
Keywords: Inverse Problems
Image Restoration
Denoising
Langevin Dynamics
Diffusion
Super Resolution
Deblurring
Compressive Sensing

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.

SNIPS: Solving Noisy Inverse Problems Stochastically
Bahjat Kawar, Gregory Vaksman, Michael Elad Computer Science Department, Technion, Haifa, Israel {bahjat. kawar, grishav, elad}-Ccs technion ac.
Abstract
In this work we introduce 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 ex- ploits pre-trained minimum mean squared error (MMSE) Gaussian denoiser. The proposed approach relies on an intricate derivation of the posterior score function that includes 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.
1
Introduction
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
Hx + 2,
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, ~ N (0, ofT). In this work we assume that both H and ao 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-inspires techniques [13] 55 [T]. local Gaussian-mixtun modeling [57] 63) and methods relying non-local self-similarity 6 9 36 [T]. More recently, deep learning techniques, a direct design of the recovery path from y to estimate of x took the lead, yielding ofate-of-the-rr results in various linear inverse problems,
and with the emergence
such as denoising
25
59 61 55), deblurring 22 48
super resolution
17 54
and other
tasks 29 28 99
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, &, and the unknown original one, x. When the degradation acute and information is irreversibly lost, image reconstruction becomes highly ill-posed problem, implying that many possible clean images could explain the given
35th Conference on Neural Information Processing Systems (NcurIPS 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-distortionn tradeoff, i.e., targeting a minimization of the error between & and x (in any metric) is necessarily accompanied by a compromised perceptual quality. As consequence, as long as we stick to the tendency to design recovery algorithms that aim for minimum MSE (or other distances), only 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 (xly) 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 D. GANs can be utilized for solving inverse problems while producing high-quality images 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 22 31 34) ) to assume noiseless measurements, 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.
(see e.g. 2
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 [8 3. The work reported in 44
20
46) utilizes the annealed
Langevin dynamics method, both for image synthesis and for solving noiseless inverse problems [
Their synthesis algorithm relies on an MMSE Gaussian denoiser (given as a neural network) for approximating 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- work be generalized for handling linear inverse problems, as in Ln[Eqattion I] in which the measurements are noisy? A partial and limited answer to this question has already been given in [1 for the tasks of results, and
image denoising and inpainting. The present work generalizes these (44 20 46 21
introduces 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 manner that facilitates constructive derivation of the conditional score function.
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 the pnior-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 position- dependent step size vector, drawing inspiration from Newton's 's method in optimization. This stabilizes the algorithm and is shown to be essential for its success.
Following earlier work 44
20
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 (xly), different runs of SNIPS on the same input necessarily yield different results, all of which valid solutions to the given inverse
42 26 and 15 23 is also relevant to this discussion, but somewhat different. We shall specifically address these papers' content and its relation work in secttotinn 2
'The work reported in 18

---

Original
Blurred
Samples from our algorithm
Mean
std
Figure I: Deblurring results on CelebA [27] images (uniform x 5 blur and an additive noise with
op = 0.1). Here and in all other shown figures, the standard deviation image is scaled by 4 for better visual inspection.
problem. This should not come surprise, ill-posedness implies that there multiple viable solutions for the same data, as has already been suggested the context of super resolution [31] 2 34). We demonstrate SNIPS 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 '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 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, should mention that using Gaussian denoisers iteratively for handling general linear inverse problems has been already proposed in the context of the Plug- and-Play-Priod (PnP) method (53) and RED (39. and their many followup papers (e-g., [60] 30 49
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 MAP or MMSE estimation.
2
Background
The Langevin dynamics algorithm using the iterative transition rule
B 38 suggests sampling from a probability distribution p (x)
X x + V lx.) + V2ozz
where N (0, I) and a is an appropriately chosen small constant. The added allows for stochastic sampling, avoiding a collapse to a maximum the distribution. Initialized randomly, after sufficiently large number of iterations, and under some mild conditions, this process converges to 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 [Equation 2 with a blurred version of it, N (0, o is synthetically injected noise. The core
vx, log p (X+), where x = + n and idea to start with a very high noise level a and gradually drop it to near-zero. all while using a step size dependent
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 brilliant relation attributed to Miyasawa 33 (also
known as Stein's 's integration by parts trick [47 or Tweedie's 's identity [12). It is given as
D (3,, 0) X2
Vx, log (x.,
(3)
3

---

where D 0) [x] is the minimizer of the MSE measure E be approximated using a denoising neural network. This facilitates the use of denoisers in Langevin dynamics as a replacement for the evasive score function.
=
||x D (R,, [[]] which can
When turning solve inverse problems, previous work suggests sampling from the posterior distri-
bution (xly) 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 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 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 [ 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 [2]] 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 I
In the following, we
present such a derivation.
3 The Proposed Approach: Deriving the Conditional Score Function
3.1 Problem Setting
We consider the problem of recovering a signal x E RN (where x p (x) and p (x) is unknown) from the observation y Hx + z, where y E RRM, H E R . < N, ~ N' (i) and H and are known Our ultimate goal is to sample from the posterior p (xly). However, since access to the score function Vx log p(xly) is not available, we retarget our goal, as explained above, to sampling from blurred posterior distributions, p (xly), where x = n and n N (0,o²1). with noise levels a 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 = UEVT, where U E R M and V R Nx. are orthogonal matrices, and Σ e M XN is a rectangular diagonal matrix containing the singular values of H. denoted as } in descending order (s1 > 82
> SM-1 > SM 2 0). For convenience of notations, we also define 8) = 0 for j = M + 1, ,1...,,,,, N. To that end, we notice that
p (x = (x = (V"xU"y).
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 V does not change its probability 224 Therefore, sampling from p (VT X 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 V log (V x/*x|U" U For simplicity, we denote
hereafter yr U 27 UTa XT VTx, ny EV- n, and Xy- VT'. Observe that with these notations, the measurements equation becomes
y Hx + 2 = UEET x2
and thus
d'y= EV- (x-n) n) - EV- UTZ
where we have relied on the relation X = x + This leads to yr ΣX- nr Zr-
(5)
In this formulation, which will aid in deriving the conditional score, our aim is to make design choices on n7 such that 27 ny has uncorrelated entries and is independent of ž7. This brings us to the formation of the synthetic annealed noise, which is an intricate ingredient in our derivations.
We 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 sequence of noise levels {oil= 0i > d2 >... or. 01+1 = O, where or is high (possibly 0i > ||x and σT, is close zero. We require that for every j such that 8j * O, there exists i, such that o4 < on and - > do- This implies Vi : 0,8; + og. which helps ease notations. SNIPS works just as well for 0:8) = do-
L+1 such that
1++ Using {o we would like to define {x, a sequence of noisy versions of x, where the noise level in X; is o., One might be tempted to define these noise additions as independent of the measurement noise 2. 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 2 in a gradual fashion. To that end, we define X = x, and for every i L, I...., 1: x; = x + 77, where 7, N (o, (o of+1) I). This results in x, = n, where = Σi= 77 ~ N (..
And now we turn to define the statistical dependencies between the measurements' noise and the artificial noise vectors 77;- Since n) and z are each Gaussian with uncorrelated entries, so are the components of the vectors EV n7, EV'n; and 2T. In order proceed while easing notations, let us focus on a single entry j in these three vectors, for which > 0, and omit this index. We denote these entries i, nT, and 27, respectively. We construct 7/,,, such that
E
Tr 2T nT.
for i 2 1
E Mmns 2r|
E
for in -1 otherwise.
This implies that the layers of noise +-1,,+1,.... Mr,, are all portions of 24 itself, with an additional portion being contained in 11 Afterwards, RI, become independent of 2r. In the case of 8j = 0, the above relations simplify to be Fl 2r] = 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
[N soo - og)) if f. > do N (0. - 3707) otherwise.
(IV'n: n. - - NT.
- -T: -
(6)
The top option refers to high values of the annealed Langevin noise, in which, despite the possible decay caused by the singular value sy, this noise is stronger than 2T. In this case. Rr,, contains all 2T and an additional independent portion of noise. The bottom part assumes that the annealed noise (with the influence of 8j) is weaker than the measurements' noise, and then it is fully immersed within 27, with the difference being Gaussian and independent.
3.2 Derivation of the Conditional Score Function
The above derivations show that the noise in Equation is zero-mean, Gaussian with uncorrelated entries and of known variance, and this noise is independent of x, Thus Equation 5) can be used conveniently for deriving the measurements part of the conditional score function. We denote X7 V x,, X x, n n, for simplicity, and turn to calculate vx, log (Krlyr). We split XT into three parts: (i) NT,0 refers to the entries j for which sj = 0; (ii) ëT,< corresponds to the entries j for which 0 < 0:8; doi and (ii) Åτ.> includes the entries j for which 0;8j > do- Observe that this partition of the entries of Xy is non-overlapping and fully covering. Similarly, partition every vector v R/ into Vo, Ve, Vs, which are the entries of v corresponding to respectively. Furthermore, we define Vy, Vz, Vy as all the entries of v except Vo. Ve. 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 &r,o. the score is independent of the measurements and given by
XT,0,
V log (Xrlyr) = (V" Vx logr (x)),
(7)
For the case of XT,>, the expression obtained is only measurements-deeeedeen
V k logp (xr
Στ OF -g(1 (yr EXT)
(8)
5

---

Original
Low-res
Samples from our algorithm
Mean
std
Super resolution results on LSUN bedroom [36 images (downscaling 4 by plain averaging and adding noise with op 0.04).
Figure 2:
Lastly, for the case of & the conditional score includes two terms one referring to the plain (blurred) score, and the other depending on the measurements,
&T. log (Srlyr) = ΣT ( - -.1-, (yr - -ExE + VVVVxlogp(ee
As already mentioned, the full derivations of equations 7 and arre detailed in the supplemental material. Aggregating all these results together, we obtain the following conditional score function:
Vxr log (RTlyt) > -0
(yT ΣXT) + (v" V log (x)
x
(10)
where (v)ly 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 Xy.u, matching the above calculations. The vector V & log, (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.
4 The Proposed Algorithm
Armed with the conditional score function in Equation 10 the Langevin dynamics algorithm can be with constant step size an annealed step size in [44. and this should converge a sample from p (&rlyr). However, for this to perform well, one should use a very small step size, implying a devastatingly slow convergence behavior. This is mainly due the fact that different entries of XT 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, suggest using step size vector a R.N. We denote A; = diag (a,), and obtain the following update formula for a Langevin dynamics algorithm:
VT x = V c.-1 VVTR, log (V + V2. cA 2,,
(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 A,, 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 II] but without the additional noise 2,, and with A, being the negative inverse Hessian of logp [ x:(yT). We calculate diagonal approximation of the Hessian, and set A, to be its negative inverse. We also estimate the conditional score function using Equation 10 and neural network. Note that this mixture of Langevin dynamics and Newton's method has been suggested in slightly different context in (4). where the Hessian was approximated using Quasi-Newte method. In our case, we analytically calculate diagonal approximation of the negative inverse Hessian and obtain the following:
0
8j =
(ax)) =
0,
0,8; > 60
(12)
o7
+
< 0,8, < 00-
$j
0

---

25%
12.5%
0.25%
Original
Degraded
Sample
Mean
Figure 3: Compressive sensing results CelebA
27
image with additive noise of do 0.1.
The full derivations for each of the three cases are detailed in the supplemental material. Using these step sizes, the update formula in Equation I1] the conditional score function in Equation 10 and neural network (x, o) that estimates the score function V log (x) we obtain tractable iterative algorithm for sampling from (X y), where the noise in X is sufficiently negligible to be considered sampling from the ideal image manifold.
Algorithm 1: SNIPS
Input:
{o:) C, y, H, og U, E. V < sddd (H)
2 Initialize Xo with random noise U [0, 1]
for i -- to L do
(Ai), - I².
(Ai) + of
of E<
(At) + - xt
7
for 1 to do
Draw 24 ~N((0.)
d E" ofll of ال
(U"y -V="-VV'x--> (VL s ox)=1
x, V. (V'x,- x c,d, (2CA32))
end
12
Xg < X+
13 end
Output: X0
Note that when we set H O and = 0. implying measurements, the above algorithm degenerates to an image synthesis, exactly 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 noisy inpainting, both cases shown in [2T] Lastly, for the choices H in 20 44 46) and with 0, the above algorithm collapses to close variant of their proposed iterative methods.
Experimental Results
network in order to estimate the score function of In our experiments we use the NCSNv2 45
the prior distribution. Three different NCSNv2 models are used, each trained separately on the training sets of: (i) images of size 64 x 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. each of the experiments, we run our algorithm 8 times, producing samples for each input. We examine both the samples themselves and their mean, which serves as an approximation of the MMSE solution, [xly]-
Recall that (x, σ) (8, f being denoising residual.
7

---

Original
Low-res
Samples from our algorithm
Mean
std
Figure 4: Super resolution results on CelebA 27 images (downscaling 1 by plain averaging and adding noise with do 0.1).
2
Original
Low res
Samples from our algutithm
Mean
sid
Figure 5: Super resolution results on CelebA 27 images (downscaling 1 by plain averaging and adding noise with do 0.1).
For image deblurring, we use uniform x 5 blur kernel, and an additive white Gaussian noise with do = 0.1 (referring pixel values in the range [0, 1)). Figure 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 block averaging filter, i.e., each non- overlapping block of pixels in the original image is averaged into one pixel in the low-resolution image. We use blocks of size x 2 or 4 x pixels, and assume the low-resolution image to include an additive white Gaussian noise. We showcase results on LSUN and CelebA in Figures 2 and
For compressive sensing, 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.
three random projection matrices with singular values of that
We calculate the average PSNR (peak signal-to- ratio) of each of the 8 samples in our experi- ments, as well as the PSNR of their mean, 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 dB, with the MMSE estimator having poorer perceptual quality but better PSNR.
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 LPIPS perceptual quality metric.
5.1 Assessing Faithfulness to the Measurements
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
8

---

Table 1: PSNR results for different inverse problems on 8 images from CelebA 27 times, and obtained 8 samples. The average PSNR for each of the samples is the first column, while the average PSNR for the mean of the 8 samples for each image is in the second one.
We ran SNIPS
Problem
Sample PSNR
Mean PSNR
Uniform deblurring
25.54
28.01
Super resolution (by 2)
25.58
28.03
Super resolution (by 4)
21.90
24.31
Compressive sensing (by 25%)
25.68
28.06
Compressive sensing (by 12.5%)
22.34
24.67
Original
Degraded
Samples from algorithm
Meun
Figure 6: Compressive sensing results on LSUN [56] tower images (compression by 25% and adding noise with do = 0.04).
unknown, 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 & by H, and calculate its difference from the given measurement y. obtaining y - H&. According to the problem setting, this difference should be an additive white Gaussian noise vector with standard deviation of σu- We examine this difference by calculating its empirical standard deviation, and performing the Pearson-D's Agostino [8] test of normality on it, accepting it Gaussian vector if the obtained p-value is greater than 0.05. We also calculate the Pearson correlation coefficient (denoted as p) 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 dg almost exactly, the Pearson correlation coefficient satisfies |p| < 0.1, and we obtain p-values greater than 0.05 around 95% of the samples (across all experiments). These results empirically show that our algorithm produces valid solutions to the given inverse problems.
6 Conclusion and Future Work
SNIPS, presented in this paper, is a novel stochastic algorithm for solving general noisy linear inverse problems. This method based on annealed Langevin dynamics and Newton's 's method, and relies on the availability of pre-trained Gaussian MMSE denoiser. SNIPS produces a random variety of high quality samples from the posterior distribution while guaranteeing their validity with respect to the given data. This algorithm's derivation includes an intricate choice of the injected anncaled 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.
the unknown given the measurements,
Extensions of this work should focus on SNIPS' limitations: (i) The need to deploy SVD decompo- sition of the degradation matrix requires considerable amount of memory and computations, and hinders the algorithm's 's scalability; (ii) The current version of SNIPS does not handle general content images, a fact that 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 sample images), and means for its acceleration should be explored.

---

7 Funding Transparency Statement
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 Ka Shing Fellowships and the Planning and Budgeting Committee of the Israel Council for Higher Education.
References
[1] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlicb. Solving inverse problems using data- driven models. Acta Numerica, 28:1-174, 2019.
[2] Y. Bahat and T. Michaeli. Explorable super resolution. Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 2716-2725, 2020.
[3] Besag. Markov Chain Monte Carlo for statistical inference. Center for Statistics and the Social Sciences, 9:24-25, 2001.
[4] S. A. Bigdeli, M. Jin, Favaro, and M. Zwicker. Deep mean-shift priors for image restoration. In Proceedings of the 31st International Conference Neural Information Processing Systems, pages 763-772, 2017.
[5] Y. Blau and T. Michaeli. The perception-distortion tradeoff. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition, pages 6228-6237, 2018.
[6] A. Buades, B. Coll, and J.-M. Morel. A non-local algorithm for image denoising. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'OH (05), volume 2, pages 60-65. IEEE, 2005.
[7] G. T. Buzzard, S. H. Chan, S. Srechari, and C. A. Bouman. Plug-and-play unplugged: Optimization-free reconstruction using consensus equilibrium. SIAM Journal on Imaging Sciences, 11(3):0011-2000, 2018.
[8] R. D' A'Aoosinno and E. Pearson. Tests for departure from normality. Empirical results for the distributions of 62 and Vb. Biometrika, 60(3):613-64) 622, 1973.
[9] A. Danielyan, V. Katkovnik, and K. Egiazarian. BM3D frames and variational image deblurring. IEEE Transactions on Image Processing, 21(4):1715-1728, 2011.
[10] C. Dong, C. C. Loy, and X. Tang. Accelerating the super-resolution convolutional neural network. In European Conference on Computer Vision, pages 391-407. Springer, 2016.
[11] W. Dong, L. Zhang, G. Shi, and X. Li. Nonlocally centralized sparse representation for image restoration. IEEE Transactions on Image Processing. 2224):1220-166 1630, 2012.
[12] B. Efron. Tweedie's formula and selection bias. Journal of the American Statistical Association, 1064496))600211114 2011.
[13] M. Elad and M. Aharon. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image Processing, 15(12):3736-3745 2006.
[14] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. Generative adversarial nets. In Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
[15] B. Guo, Y. Han, and J. Wen. AGEM: solving linear inverse problems via deep priors and sampling. Advances Neural Information Processing Systems, 32:547-558, 2019.
[16] H. Gupta, K. H. Jin, H. Q Nguyen, M. T. McCann, and M. Unser. CNN-based projected gradient descent for consistent CT image reconstruction. IEEE Transactions on Medical Imaging, 37(6): 1440-1453, 2018.
10

---

[17] M. Haris, G. Shakhnarovich, and N. Ukita. Deep back paoj-prooection networks for super-resolution. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1664-1673 2018.
[18] J. Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33. pages
Curran Associates, Inc., 2020.
[19] C. M. Hyun, H. P. Kim, S. M. Lee, S. Lee, and J. K. Seo. Deep learning for undersampled MRI reconstruction. Physics in Medicine & Biology, 63(13):135007, 2018.
[20] Z. Kadkhodaie and E. P. Simoncelli. Solving linear inverse problems using the prior implicit in a denoiser. arXiv preprint arXiv:2007.13640, 2020.
[21] B. Kawar, G. Vaksman, and M. Elad. Stochastic image denoising by sampling from the posterior distribution. arXiv preprint arXiv:2101.09553 2021.
[22] O. Kupyn, Т. Martyniuk, J. Wu, and Z. Wang. DeblurGAN-v2: deblurring (orders-on ofaes--ofmmaa of-magnitude) faster and better. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 8878-8887, 2019.
[23] R. Laumont, V. De Bortoli, A. Almansa, J. Delon, A. Durmus, and M. Pereyra. Bayesian imaging using plug & play priors: When Langevin meets Tweedie. arXiv preprint arXiv:2103.04715, 2021.
[24] G. Lebanon. Probability: The Analysis of Data, volume pages 104-107. CreateSpace Independent Publishing Platform, 2012.
[25] S. Lefkimmiatis. Non-local color image denoising with convolutional neural networks. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition, pages 3587-3596, 2017.
[26] H. Li, Y. Yang, M. Chang, H. Feng. Z. Xu, Li, and Y. Chen. SRDiff: single image super- resolution with diffusion probabilistic models. arXiv e-prints, pages arXiv-2104, 2021.
[27] Z. Liu, P. Luo, X. Wang, and X. Tang, Deep learning face attributes in the wild. In Proceedings of the IEEE International Conference on Computer Vision, pages 3730-3738. 2015.
[28] A. Lucas, M. Iliadis, R. Molina, and A. K. Katsaggelos. Using deep neural networks for inverse problems in imaging: beyond analytical methods. IEEE Signal Processing Magazine, 35(1): 20-36, 2018.
[29] M. T. McCann, K. H. Jin, and M. Unser. Convolutional neural networks for inverse problems in imaging: A review. IEEE Signal Processing Magazine, 34(6):85-95, 2017.
[30] T. Meinhardt, M. Moller, C. Hazirbas, and D. Cremers. Learning proximal operators: Using denoising networks for regularizing inverse imaging problems. In Proceedings of the IEEE International Conference on Computer Vision, pages 1781-1790, 2017.
[31] S. Menon, A. Damian, S. Hu, N. Ravi, and C. Rudin. PULSE: Self-supervised photo upsampling via latent space exploration of generative models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2437-2445, 2020.
[32] K. Miyasawa. An empirical Bayes estimator of the mean of a normal population. Bull. Inst. Internat. Starist., 38:181-188, 1961.
[33] G. Ohayon, T. Adrai, G. Vaksman, M. Elad, and P. Milanfar. High perceptual quality image denoising with a posterior sampling CGAN. arXiv preprint arXiv:2103.04112, 2021.
[34] Peng and K. Li. Generating unobserved alternatives: A case study through super-resolution and decompression. arXiv preprint arXiv:2011.01926 2020.
[35] A. Radford, L. Metz, and S. Chintala. Unsupervised representation learning with deep con- volutional generative adversarial networks. In 4th International Conference on Learning Representations, 2016.

---

[36] I. Ram, M. Elad, and I. Cohen. Image processing using smooth ordering of its patches. lEEE Transactions on Image Processing, 22(7):2764-277 2013.
[37] Ravishankar, C. Ye, and A. Fessler. Image reconstruction: From sparsity data addativee methods and machine learning. Proceedings of the IEEE, 108(1)886-109 2019.
[38] G. O. Roberts, R. L.. Tweedie, et al. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341 1996.
[39] Y. Romano, M. Elad, and P. Milanfar. The little engine that could: Regularization by denoising (RED). SIAM Journal Imaging Sciences, 10(4): 1844-18444, 2017.
[40] A. Rond, R. Giryes, and M. Elad. Poisson inverse problems by the plug-and and-play scheme. Journal of Visual Communication and Image Representation, 40896-108, 2016.
[41] E. Ryu, Liu, S. Wang, X. Chen, Z. Wang, and W. Yin. Plug-and-play methods provably converge with properly trained denoisers. In International Conference Machine Learning. pages 5546-5557. PMLR, 2019.
[42] C. Saharia, J. Ho, W. Chan, T. Salimans, J. Fleet, and M. Norouzi. Image super-resolution via iterative refinement. arXiv preprint wwriiv22104.0777 2021.
[43] U. Simsekli, R. Badeau, T. Cemgil, and G. Richard. Stochastic Quasi-Newton Langevin Monte Carlo. In International Conference on Machine Learning, pages 642-651. PMLR, 2016.
[44] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. In Advances Neural Information Processing Systems, pages 11918-11934 2019.
[45] Y. Song and Ermon. Improved techniques for training score-based generative models. In Advances in Neural Information Processing Systems, 33, 2020.
[46] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based generative modeling through stochastic differential equations. International Conference on Learning Representations, 2021.
[47] C. M. Stein. Estimation of the mean of a multivariate normal distribution. The annals of Statistics, pages 1135-1151, 1981
[48] M. Suin, K. Purohit, and A. Rajagopalan. Spatially-attentive patch-hierarchical network for adaptive motion deblurring. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages
2020.
[49] Y. Sun, B. Wohlberg, and U. S. Kamilov. An online plug-and-play algorithm for regularized image reconstruction. IEEE Transactions Computational Imaging, 5(3):395-40))) 2019.
[50] T. Tirer and R. Giryes. Image restoration by iterative denoising and backward projections. IEEE Transactions on Image Processing, 28(3):1 28(3):1220-122 1234, 2018
[51] G. Vaksman, M. Zibulevsky, and M. Elad. Patch ordering as a regularization for inverse problems in image processing. SIAM Journal Imaging Sciences, 9(1):287-319 2016.
[52] G. Vaksman, M. Elad, and P. Milanfar. LIDIA: Lightweight learned image denoising with instance adaptation. In Proceedings of the IEEE/CVE Conference on Computer Vision and Pattern Recognition Workshops, pages 525, 2020.
[53] S. V. Venkatakrishnan, C. Bouman, and B. Wohlberg. Plug-and-play priors for model based reconstruction. In 2013 IEEE Global Conference on Signal and Information Processing, pages 945-948. IEEE, 2013.
[54] X. Wang, K. Yu, S. Wu, J. Gu, Y. Liu, C. Dong, Y. Qiao, and C. Change Loy. ESRGAN: enhanced super-resolution generative adversarial networks. In Proceedings of the European Conference on Computer Vision (ECCV) Workshops, 2018
[55] Yang, Wright, T. S. Huang, and Ma. Image super-resolution via sparse representation. IEEE transactions image processing, 19(11):26112287 2010.
12

---

[56] F. Yu, A. Seff, Y. Zhang, S. Song, T. Funkhouser, and J. Xiao. LSUN: construction of a large-scale image dataset using deep learning with humans in the loop. arXiv preprint arXiv: 1506.03365, 2015.
[57] G. Yu, G. Sapiro, and S. Mallat. Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity. IEEE Transactions on Image Processing, 21(5):2481-284 2499, 2011.
[58] Zhang and B. Ghanem. ISTA-Net: interpretable optimization-inspinn deep network for image compressive sensing. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1828-1837, 2018.
[59] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang. Beyond Gaussian denoiser: Residual learning of deep CNN for image denoising. IEEE Transactions on Image Processing, 26(7): 3142-315 3155, 2017.
[60] K. Zhang, W. Zuo, Gu, and L.. Zhang. Learning deep CNN denoiser prior for image restoration. In Proceedings of the IEEE Conference Computer Vision and Pattern Recognition, pages 3929-3938, 2017.
[61] K. Zhang, W. Zuo, and L. Zhang. FFDNet: toward fast and flexible solution for CNN-based image denoising. IEEE Transactions Image Processing, 27(9):4608 4622, 2018.
[62] R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang. The unreasonable effectiveness of deep features as a perceptual metric. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 586-595, 2018
[63] D. Zoran and Y. Weiss. From learning models of natural image patches whole image restoration. In 201/ International Conference Computer Vision, pages 479-486. IEEE, 2011.
13


Reviews:
[{"review_id": "sDXy9rV5mYn", "reviewer": "Reviewer_VyNV", "summary": "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.\n\n\n---Edit after author responses---\n\nHaving read the other reviews and the responses posted by the authors, I am reducing my original scores from 6 --> 4. \nMy 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.", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Sections 1 and 2 are well written. The problem is clearly stated and the method is well motivated with a good discussion on prior work. However, the authors do not clearly state what the contributions of the current paper are. I suggest the authors highlight these in a subsection at the end of Section 2. It appears Section 3 and 4 are the main contributions of this work, but they fall short in delivering the clarity desired from any original piece of work.\n\n\nSection 3 is technical and intricate but lacks clarity throughout.\nIt is not clear why access to the true score function \\nabla log p(x|y) is not available. The assumptions that lead to these limitations need to be made clearer at the beginning of Section 3.\nThe derivations in Section 3.1 are not explained well. It makes very little sense to the reader as to why the specific noise levels and correlations are chosen. The authors ought to provide a lot of intuition here. In Section 3 it is very unclear what quantities are known for the problem, which of those are assumed to be known for the purpose of the analysis, and which can be assumed to be known while designing the algorithm.\n\nSection 2 makes a lot of notational assumptions which are very hard to follow. \n\nWhile the derivations in Section 3.2 are pushed to the appendix, the authors need to remind the reader what the significance of each of the equations 7-11 is by mentioning in words what these equations mean.\nIn Algorithm 1 the notation of the dagger and >, /> appear non-standard. I urge the authors to add a caption to the algorithm which clarifies this notation. \nLine 111: Is the algorithm in equation (2) called the Unadjusted Langevin algorithm (ULA)? If so perhaps the authors should mention that name.\n\n\n\nAdditional comments/typos:\nLine 47: GAN-based and Langevin sampling --> based on GANs and Langevin sampling\nLine 182: The quantity \\eta_T is not defined explicitly. Please do so\nNotation issue: \\eta_T, n_T are defined as V^T \\eta and V^T n respectively. \n\\eta_T,i and n_T,i however is a component of \\Sigma V^T\\eta_i and \\Sigma V^T n_i which is inconsistent. Please change the notation to make it easier to read.\nEquations following line 182, please remind the reader that i_1 is a function of j.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations are adequately addressed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "O1jeQZAE3Ok", "reviewer": "Reviewer_GwuL", "summary": "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.\n", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "__Strengths__:\n1. In theory, existing work cannot handle the challenges posed in the paper. The authors define their problem clearly, and they are able to derive approximations to the posterior distribution given access to a generative model that estimates annealed priors.\n\n__Weaknesses__:\n\n1. The key point of this paper is that the distribution $p(y|x)$ must be modified when you do annealing of the distribution $p(x)$, since this can introduce dependencies between the measurement noise and the annealing noise. Prior work would use Bayes rule to calculate the posterior distribution under annealing, and hence would get an incorrect functional form for the annealed posterior distribution. Is the error in your approximated distribution smaller in some helpful way? Overall, since there are no error bounds for sampling from annealed posteriors that converge to the true posterior, I fail to see the significance of an exact derivation of the true annealed posterior. In fact your derivation only gets an approximation to the annealed posterior, and hence will introduce its own error.\n\n1. The writing is extremely imprecise and the contributions are unclear.  Other than statements that prior work cannot handle ill-posed operators and noise, there isn't any empirical or theoretical evidence that prior work fails or that the proposed approach fixes this problem.\n\n1. The authors say that the rational behind sampling in the SVD domain is in the appendix, but this only contains derivations for the conditional likelihood in the SVD domain. There is no theoretical justification for sampling in the SVD domain.\n\n1. There are literally no comparisons to baselines. Reference [1] shows that a direction application of Bayes' rule gives good empirical results for inpainting and compressed sensing, while prior work cited in the paper has shown it for super-resolution, denoising, etc. Does your method have any empirical advantage?\n\n__Missing references__:\n\n[1] shows that posterior sampling is optimal for compressed sensing and needs to be cited.\n\n[1] Jalal, Ajil, Sushrut Karmalkar, Alexandros G. Dimakis, and Eric Price. \"Instance-Optimal Compressed Sensing via Posterior Sampling.\" arXiv preprint arXiv:2106.11438 (2021).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No societal concerns.", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "97iErwm_MMM", "reviewer": "Reviewer_kSnH", "summary": "The paper discusses an algorithm for sampling from the posterior distribution of a linear inverse problem with the white Gaussian noise, given the MMSE Gaussian denoiser. \n\n\n---Edit after author responses---\n\nAfter reading the reviews and responses I decided to decrease the score (from 7 to 5) mainly because of insufficient evaluation of the proposed method and issues with the derivations in Section 3.   ", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper has a nice literature overview and is overall well written. The proposed method combines the newly proposed blurring of the samples by the additional noise with dependence on the measurement noise and the Langevin dynamic algorithm performed in the SVD domain of the known degradation matrix with the correction of the speed size. The sampling relies on an available pre-trained MMSE Gaussian denoiser. The experimental results demonstrate realistic sampling on image deblurring, super-resolution, and compressive sensing tasks.\n\nThough the authors discuss the limitations of the proposed algorithm, in the numerical results the score function is estimated on the base of the same dataset, from which the linearly transformed and the blurred pictures come from. The resulting sampled pictures look really nice, but  which could be applications of the presented results? Would it be possible to assume a bit more realistic setting, when the source database is not available?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discuss the limitations of the work.  ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "67zubYOkGZ", "reviewer": "Reviewer_d1T5", "summary": "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.\n", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This work proposes a novel and interesting method for the chosen task: sampling from the posterior distribution in a linear inverse problem affected by noise. The motivation of the method and the results are quite compelling. There are, however, some issues with the description of the algorithm that should be addressed. If these changes are made, I would recommend the publication of this manuscript as part of the proceedings. More detailed comments follow below.\n\nFirst, the convention used for the SVD is a bit non-standard in that the singular value matrix Σ is rectangular (which results in odd expressions such as Σ^T and (Σ^†)^T). I would suggest using the convention that Σ is M×M and V is N×M. Some explanation for the rewrite in eq. (5), where y is written in terms of x-tilde instead of x would also be helpful. Similarly, it would be good to explain why it is desirable that z_T – Σn_T be white and independent of x_T-tilde. Throughout, explanations for various derivations would help the reader to follow the choices made and why. This is especially true for derivations and calculations that are relegated to the supplementary material. Here, it would be good to at least give an intuition for why the result should be true.\n\nThe most problematic part is the second half of Section 3.1, the derivation of the annealing noise term. First, the requirement that σ_i s_j ≠ σ_0 for all i, j is a bit worrisome. While it can be easily circumvented by an appropriate choice of variances σ_i, such an inequality constraint suggests an instability in the underlying relationships. The question is what happens when σ_i s_j is close to, but not equal to σ_0. While exact equality is very unlikely (especially given the floating-point representations used in practice), there is nothing preventing these values from being rather close to one another. What would happen if this were the case? Would convergence suffer? Do the calculations come out the same, with factors cancelling out in the limit?\n\nContinuing, the notation i_1 is a bit odd. Where does the subscript one come from? Also, the dependence on j should be made clear, so that i_1 would be come i_1(j) (or simply i_j?). The definition is also only given for s_j ≠ 0. What happens for s_j = 0? This becomes relevant later (line 185) when the special case s_j = 0 is singled out, but this is in fact undefined since i_1 is not given in this case.\n\nThe derivation of the noise term itself is quite confusing. First (on line 176), it is stated that η_i is defined as a normally distributed random vector with covariance matrix equal to a multiple of the identity matrix. Then, (on line 182), it is stated that η_{T,i} (the jth component of Σ V^T η_i) is *constructed* to satisfy certain moment conditions. As far as I can tell, what the authors are trying to say here is that η_i has a certain *marginal* distribution, which is a normal distribution with the specified covariance matrix, but that it is also jointly normally distributed with z_T and that the cross-covariance matrix (between η_{T,i} and z_T) is being specified. Nonetheless, this is quite an indirect way of stating the relationship between these random variables.\n\nThis becomes more obvious after (on line 183), where the authors state that “this implies that the layers of noise … are all portions of z_T itself”. Again, what I can tease out from this is that the η_{T,i} variables are simply scaled versions of z_T in this case. I think this section would be much better served by explicitly constructing η_i (or η_{T,i}) from z_T (plus an additional independently distributed component). Not only would this be easier to follow, but it would more clearly illustrate the purpose of the construction. The same comment applies to eq. (6). Simply stating that “it can be shown that” (and leaving the calculation to the reader) reduces the clarity of the argument. If the calculations have to be relegated to the supplementary material, at least a brief motivation would be useful here.\n\nSome minor comments follow:\n\n– On line 80, what does it mean for “the entries in the derived score switch mode”? There has been no discussion of modes prior to this.\n\n– It would be helpful for readers less familiar with Langevin dynamics to explain on line 130 why the conditional score function would not be tractable in the general case, since this is a problem that is solved by the proposed algorithm.\n\n– In the problem formulation at the beginning of Section 3.1, it would be good to mention the use of a prior to complete the description of the statistical model (since this is needed to define the posterior).\n\n– The superscript 2 on line 144 is oddly placed, coming right in the middle of some mathematical expressions. It would be better to place it at the end of the sentence.\n\n– In many places (such as lines 179 and 194), the word “white” is used when something more like “heteroscedastic” is meant. A white noise model implies that the components are uncorrelated and that variances of all the components are equal. While this may be the case for η_i, it is not true for Σ V^T η_i. In this case, the covariances are zero, but the variances are different.\n\n– In Sections 3.1 and 3.2, the authors switch between i ≥ i_1, σ_i s_j < σ_0 and s_j < σ_0 / σ_i, all of which are equivalent conditions. I suggest the authors pick one of these characterizations and use it throughout to reduce confusion.\n\n– When giving the step size heuristic in eq. (12), it would be good to provide a brief motivation for why the data term does not enter into the Hessian approximation, since the expression of the score function in eq. (10) indicates that this would be the case. As stated previously, it would also be good to provide some intuition for why these step sizes make sense.\n\n– There is a space missing between the numbers and “dB” on lines 268 and 271.\n\n– What does it mean that “SNIPS does not handle general content images”?\n\n– Several titles in the bibliography are improperly capitalized (probably missing some {} in the BibTeX file), including [3], [7], [8], [14], [21], [24], [30], [35], [36], [40], [48], [50], [53], [54], [55], [56], [57].\n\nUPDATE:\nAfter reading the other reviews and the author responses, I have decided to lower my score from 7 to 5 due to the lack of comparisons to established baselines, which make it difficult to establish the usefulness of the method.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discuss the limitations and necessary extensions in Section 6 of their work. These are reasonable and highlight areas of improvement that should be pursued. As this is a theoretical work, there is no expectation of negative societal impact and the authors have accordingly not discussed it.\n", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical issues are raised by the manuscript.\n", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper solves linear inverse problems via an approximate MMSE estimator using pre-trained Gaussian denoisers. The results look very good since the paper uses NCSNv2 which is an excellent generative model. 

The two major problems are: 
There are really no comparisons with previous experimental baselines. 
This is a major limitation for any scientific paper and the authors need to strengthen their evaluation section by discussing how prior work performs for the problems they solve. There are dozens of papers that solve deblurring, inpainting and compressed sensing problems using deep generative models. (Many papers use NCSNv2 and the more updated Score_SDE by Song et al. for various inverse problems and many others use other types of deep generators). 

Typically visual comparisons and comparisons in MSE and SSIM would help the readers understand when the current method outperforms previous work. 

Further, the technical contribution of the paper was actually not clear to the reviewers. We would encourage the authors to clarify their presentation based on the reviewer feedback and revise their manuscript. 

Decision: Accept (Poster)