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 552, 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 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 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 secttttinn 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 (3x, 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 + z = 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/r,, such that
E
Tr 2T nT.
for i 2 1n
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
Xr,0,
V log (Xrlyr) = (V" Vx logr (x)),
(7)
For the case of XT,>, the expression obtained is only measurements-deeeedeen
V k logp (y)
Στ 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 [4). 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-
jj
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 ال x; V. (V'x,- c,d, (2CA32.)
(U"y -V="-V''x--> (VL s ox)=1
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-non 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- 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-supervieed 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 adaativee 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": "wEyaUkJ2YqE", "reviewer": "Reviewer_y4bh", "summary": "This paper proposes to combine Langevin dynamics and Newton's method to arrive at a posterior sampling algorithm, where the prior distribution is given by an MMSE Gaussian denoiser. The proposed method is evaluated on multiple imaging inverse problems like deblurring, super-resolution, etc.", "questions": "", "limitations": "", "rating": 7, "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: the method proposed in this paper is original and technically sound. \n\nThe authors utilize denoiser-based deep prior and blurred posterior function, which is very similar to the approach used in (Bigdeli et al. 2017). It could be helpful if similarities and differences with that work are discussed.\n\nWeakness:\n- The introduction section extensively discusses the properties of the proposed method before actually introducing the method. IMHO the introduction should be precise, and more space can be dedicated to the experimental results section.\n- The experimental results section heavily relies on subjective results. Although understandably, objective metrics like PSNR cannot reflect the perceptual quality, it is still desirable to include some objective results comparing to some baseline methods in the main text.\n- In the bibliography, venues should be properly capitalized (e.g. line 335)\n\nReferences:\nBigdeli, Siavash A., et al. \"Deep mean-shift priors for image restoration.\" Proceedings of the 31st International Conference on Neural Information Processing Systems. 2017.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Have the authors adequately addressed the limitations and potential negative social impact of their work: Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8 hours"}, "extra_text": {}}, {"review_id": "mNjQov44yjd", "reviewer": "Reviewer_KMnW", "summary": "This paper proposes a method to sample from the posterior density in the Bayesian formulation of linear (imaging) inverse problems using Langevin dynamics. The main contribution of the paper is a way to generalize the idea (which is not new) beyond diagonal observation operators (denoising and inpainting), by resorting to an SVD of the operator.\n\n ", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "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 paper proposes a method to sample from the posterior density in the Bayesian formulation of linear (imaging) inverse problems using Langevin dynamics. The main contribution of the paper is a way to generalize the idea (which is not new) beyond diagonal observation operators (denoising and inpainting), by resorting to an SVD of the operator and some \"intricate\" derivations. \n\nThis is a very nice paper and the proposed formulation is solid and constitutes an important contribution to the topic of imaging inverse problems. However, there are a few minor issues that should be addressed. \n\nA fundamental limitation of the proposed method, which should be clearly acknowledged, is that sampling from the posterior is not applicable/acceptable in many domains. Look, for example, at the samples in the compressive sensing or super-resolution experiments; in some samples, the teeth are visible in others not. If this was a medical image reconstruction/restoration scenario, this could correspond to some samples showing a tumor, others not. In such a context, what is required is THE best possible reconstruction (under some optimality criterion), not samples from the posterior, which can vary significantly. \n \nUnlike what the authors claim, it is not true that \"most image restoration techniques aim to form an MMSE estimator\". In fact, the vast majority of image restoration approaches are formulated as optimization problems, which can be interpreted as MAP estimators or, from a regularization theory point of view, as regularized inversion of the (ill-posed) observation operator. There are a few attempts at leveraging MAP/regularization approaches to yield MMSE estimates, such as https://epubs.siam.org/doi/abs/10.1137/120902276, but it is certainly not the majority. \n\nThe argument that \"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\" is not new, but it certainly requires some sort of justification, as it is far from evident that this is true. Unless the posterior is strongly multimodal, the posterior expected value could very well itself be a sample. As such, without any comment about the properties of the posterior, this is nothing but an empirical observation. \n\nIn lines 104-107, the authors write “that using MMSE Gaussian denoisers iteratively for handling general linear inverse problems has been already proposed in the context\nof the Plug-and-Play-Prior (PnP) method and RED, and their many follow-up papers. However, none of the cited methods use explicitly MMSE denoisers, but general-purpose denoisers. In fact, at the heart of the PnP is the idea that you can replace the proximity operator of some regularizer with a (Gaussian) denoiser, but this equivalence sees it as a MAP denoiser, not an MMSE one. The only work where PnP is used with an explicit MMSE denoiser is https://doi.org/10.1109/TIP.2018.2869727\n\nThe derivations and explanations in lines 149 to 166 are unnecessarily “intricate” (to use a word that the authors have used a lot). If you simply start with the equation between lines 161 and 162 and use the definitions in line 160, then everything (namely equation (5)) falls directly from there. In particular, you don’t need equation (4), which anyway is not referred to in the rest of the paper and is trivial from the formulation between lines 161 and 162. \n\nIt would be interesting to see if the formulation (and the resulting computation) is significantly simplified in the case where H is a periodic convolution (thus with M=N), in which case H is diagonalized by the discrete Fourier transform and, consequently, multiplications by U and V can be done with O(N log N) cost using FFTs.\n\nIn the experiments, the main claim is that, although having a lower PSNR (which is natural, since they are not the MMSE estimate, which maximizes the PSNR), the samples are more \"visually appealing\". This claim could be, at least in part, supported on other quality metrics, such as SSIM (structural similarity index measure), which correlates better with human quality perception.\n\nFinally, there is no point in showing the ”degraded” images in the compressive sensing experiments, since these are not degraded images, but observations that are meaningless to a human observer. \n\n\n\n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Some limitations have not been addressed, as mentioned in the main review above. Namely, that this type of approach may not make sense in some contexts, such as for medical image reconstruction. My recommendation is that this limitation should be briefly discussed in the paper.\n ", "review_text": "", "extra_scores": {"ethical_concerns": "I don't see any ethical issues in the paper. \n", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "XXSqDkwb2X", "reviewer": "Reviewer_nVFq", "summary": "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.", "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": "\nReview Summary\n--------------\n\nThis paper proposes a technique based on annealed Langevin dynamics to solve noisy linear inverse problems. IMO, the novelty of this paper is limited as it seems to boil down to applying [19] in the spectral space of the operator $H$ (**M1**). Moreover, the numerical experiments are quite limited, with mainly qualitative results and no comparison with other methods (**M2**). For these reasons, I recommend rejecting this paper. As I might be wrong for the spectral space part, I will carefully read the author response and revise my score if necessary.\n\n\nMajor comments and questions\n----------------------------\n\n* **M1:** One main concern for this paper is its limited novelty. The posterior sampling using Langevin dynamics has been described in [18] and [19] proposed a first adaptation to the inverse problem for specific operator $H$. The core contributions of the papers are thus computing the conditional score function in the spectral space of the operator and proposing to use per-coordinate step sizes:\n    - **A:** The authors states that the extension of Langevin dynamics to spectral space is far from trivial but I am not really convinced about this. Looking at the equation (10), we can see that the dynamic is independent for each coordinate in the spectral space and the propose algorithm boils down to the one in [19] but applied to all variables in spectral space, where the noise levels are different for each coordinate. Equation bellow l.182 simply states that the added noise $\\eta_{T,k}$ should sum to $z_T$ (*i.e.* the standard assumption but in the spectral space). The most complex part was to show that $\\nabla_{\\tilde x_{T, >}}\\log p(\\tilde x_{T, 0}|\\tilde x_{T, \\not 0}) = 0$ but this result is taken as an assumption (l.60 in the supplementary) while it should probably holds as $V^T$ is unitary and thus the added noise are independent for each coordinate in spectral space. Overall, it seems that the proposed derivation could be simplified a lot while retaining their efficiency.\n    - **B:** For the step-sizes adaptation, it makes sense that it should be adapted on each coordinates as the noise levels are not the same. However, the benefit of using such a scheme is not demonstrated in the paper.\n\n* **M2:** Another major concern for this paper is with the experimental results. Indeed, the paper does not propose any comparison with other methods. For inverse problem solution sampling, the authors could compare their results with the one from [30]. The mean results could also be compared with Plug-and-Play or Regularization by denoising methods, in particular with the PSNR from Table.1. Without such validation, it is hard to assess wether this method provides good results or not. Also, adding an ablation study to show the benefit of the proposed step-size compared to using the classical one would help convincing that it indeed accelerate the dynamic.\n\n\nMinor comments, nitpicks and typos\n----------------------------------\n\n\n* **m1:** (Numerical experiments) The part of the LSUN and CelebA datasets used to train the NCSNv2 models used to approximate the prior distribution are not properly described in the text. From the code, it seems that the training set is used and that the inverse problems are constructed from the test/validation sets but this should be properly described in the paper.\n\n* **m2:** (Numerical experiments) The parameters that are used in the experiments are not described in the paper. In particular, the sequence of sigma that is chosen (it seems to be geometrical from the code) and the step sizes\n\n* **m3:** (Eq.(10) & Algorithme.1 - l.9) the notation $|\\sigma_0^2 I - ... |^\\dagger$ is confusing. I think it would be clearer to use $(\\sigma_0^2 I - ... )^\\dagger$ as done in Eq.(8).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "* **L1:** One major limitation of such method is that it is limited to simple noise models as one needs to derive an expression for the score function. Moreover, the noise level must be known in advance for the method to work. These limitation should be mention in the section.6. or discuss in the experimental results, for instance showing the results when the noise is more complex or only an approximation of the noise level is given.\n\n* **L2:** Another limitation that should be mentioned is that this technique has only been tested for inverse problem in imaging. For more complex inverse problems where the ground truth is unknown and the noise more complex, such as in astrophysics or neuroscience, there is little chances to be able to apply such techniques to recover the signal of interest.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": ">12hours"}, "extra_text": {}}, {"review_id": "WJ40-G0-VQC", "reviewer": "Reviewer_YnSg", "summary": "The paper describes an apparently new approach to solving noisy inverse problems.\n\nIt shows how the approach can be applied to three different noisy inverse problems:\nwhich it calls image deblurring, super-resolution, and compressive sensing. \n\nIt validates the approach on synthetic data which is synthesized\n(eg blurred and noise added) starting from publicly available imagery.\n\nIt is carefully written and I find it readable.", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "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 describes an apparently new approach to solving noisy inverse problems.\n\nIt shows how the approach can be applied to three different noisy inverse problems:\nwhich it calls image deblurring, super-resolution, and compressive sensing. \n\nIt validates the approach on synthetic data which is synthesized\n(eg blurred and noise added) starting from publicly available imagery.\n\nIt is carefully written and I find it readable.\n\nOriginality: I have seen work on langevin dynamics in inference problems before,\nbut I think that work was theoretically focused and the computational efforts\nwere just preliminary efforts rather than mature tool construction.\nI believe the originality here is in the effort to deliver something\nthat's been carefully derived, using sophisticated identities\nrather than handwaving, and carefully tested on (stylized approximations to) \n'real-world' problems.\n \nQuality: The paper is very well written at the narrative level. \n\nClarity: I find it very clear.\n\nSignificance: I am unable to judge whether the proposed method is an improvement over\nexisting methods. Does this method do better than the state of the art on the problems\nin its purview (denoising, super-resolution, compressed sensing) where so many previous\nefforts have been made and carefully documented? I don't know.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I think what they call compressed sensing is a not at all compressed sensing.\nThey are not doing random undersampling either in Fourier space nor do they\nseem to be doing random undersampling in some other randomly orthogonally transformed domain.\nI think they are simply dealing with undersampled data in the original pixel domain.\nPossibly this is a real problem in certain settings, but it shouldn't be\ncalled compressed sensing.\n\nI think the identity (3) that they call Brilliant is well known to Mathematical Statisticians for\ndecades as \"Stein's Integration by Parts trick\" after Charles Stein, a leading mathematical\nstatistician of the 20th century. \nSee for example Stein, Charles M. (November 1981). \"Estimation of the Mean of a Multivariate Normal Distribution\". The Annals of Statistics. 9 (6): 1135–1151. doi:10.1214/aos/1176345632 (many other statistics papers use this idea, \nphysicists also know it well by other names)", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hrs"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper considers noisy linear inverse problems y = Hx +noise, and specifically, the problem of generating samples from the posterior p(x|y) using a Langevin dynamics approach. This has been considered before quite recently, but the novelty of the current paper is working with noise as well as extending ref [19] to work with more general linear operators H.

This paper generated a lot of discussion among the reviewers and between the reviewers and author(s). We all found the paper quite interesting!  There were many suggestions for improvement, even from reviewers with very high scores.  Some common themes are that some of the derivations are more complicated than necessary, and missing some prior work (e.g., from Stein), and that the innovation over [19] isn't strong enough.

Despite these weaknesses, the majority opinion is that the results are strong enough to warrant publication.  One reviewer was impressed that this work is able to make Langevin dynamic sampling so practical -- the actual experiment results shown in the figures speak for themselves.

I strongly recommend to the authors to take into account the reviewer comments in their final version of the paper. The objective is not to publish a paper in NeurIPS, but rather to publish the best possible paper that will have a large impact. We had a strong group of reviewers with good comments, so I suggest following their advice when possible.
Decision: Accept (Poster)