On the Posterior Distribution in Denoising: Application to Uncertainty Quantification
Authors: Hila Manor
Tomer Michaeli
Venue: ICLR 2024 poster (2024)
Keywords: Gaussian Denoising
Posterior Moments Estimation
Uncertainty Quantification
Uncertainty Visualization
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We derive a fundamental property of the posterior distribution in Gaussian denoising, and use it to propose a new way for uncertainty visualization, which requires no training or fine-tuning.

Abstract: Denoisers play a central role in many applications, from noise suppression in low-grade imaging sensors, to empowering score-based generative models. The latter category of methods makes use of Tweedie's formula, which links the posterior mean in Gaussian denoising (*i*.*e*., the minimum MSE denoiser) with the score of the data distribution. Here, we derive a fundamental relation between the higher-order central moments of the posterior distribution, and the higher-order derivatives of the posterior mean. We harness this result for uncertainty quantification of pre-trained denoisers. Particularly, we show how to efficiently compute the principal components of the posterior distribution for any desired region of an image, as well as to approximate the full marginal distribution along those (or any other) one-dimensional directions. Our method is fast and memory-efficient, as it does not explicitly compute or store the high-order moment tensors and it requires no training or fine tuning of the denoiser. Code and examples are available on the project [website](https://hilamanor.github.io/GaussianDenoisingPosterior/).

Published as a conference paper at ICLR 2024
ON THE POSTERIOR DISTRIBUTION IN DENOISING: APPLICATION TO UNCERTAINTY QUANTIFICATION
Hila Manor
Tomer Michaeli
Faculty of Electrical and Computer Engineering
Faculty of Electrical and Computer Engineering
Technion Israel Institute of Technology hi la manor campus. technion. ac il
Technion Israel Institute of Technology
tomer m@ee. technion. ac.il ill
ABSTRACT
Denoisers play a central role in many applications, from noise suppression in low- grade imaging sensors, to empowering score-based generative models. The latter category methods makes use of Tweedie's formula, which links the posterior mean in Gaussian denoising (i.e., the minimum MSE denoiser) with the score of the data distribution. Here, we derive a fundamental relation between the higher- order central moments of the posterior distribution, and the higher-order deriva- tives of the posterior mean. We harness this result for uncertainty quantification of pre-trained denoisers. Particularly, we show how efficiently compute the prin- cipal components of the posterior distribution for any desired region of an image, as well as approximate the full marginal distribution along those (or any other) one-dimensional directions. Our method is fast and memory-efficient, as it does not explicitly compute or store the high-order moment tensors and requires no training fine tuning of the denoiser. Code and examples are available on the project website.
INTRODUCTION
Denoisers serve as key ingredients in solving wide range of tasks. Indeed, along with their tradi- tional use for noise suppression (Aharon et al., 2006; Buades et al., 2005; Dabov al., 2007; Krull et al., 2019; Liang et al., 2021; Portilla et al., 2003; Roth & Black, 2009; Rudin et al., 1992; Zhang al., 2017a; 2021), the last decade has seen steady increase in their use for solving other tasks. For example, the plug-and plug-and-play method (Venkatakrisshan et al., 2013) demonstrated how a denoiser can be used an iterative manner to solve arbitrary inverse problems (c.g., deblurring, inapinting). This approach was extended by many, and has led to state-of-the-art results various restoration tasks (Brifman et al., 2016; Romano et al., 2017; Tirer & Giryes, 2018; Zhang et al., 2017b). Sim- ilarly, the denoising ccore-matching work (Vincent, 2011) showed how denoiser constructing a generative model. This approach was later improved (Song & Ermon, 2019), and highly related ideas (originating from (Sohl-Dickstein et al., 2015)) served as the basis for diffusion models (Ho et al., 2020), which now achieve
be used for
results on image generation.
Many of the uses of denoisers rely on Tweedie's formula (often attributed to Robbins (1956), Miya- sawa et al. (1961), Stein (1981), and Efron (2011)) which connects the MSE-optimal denoiser for white Gaussian noise, with the score function (the gradient of the log-probability w.r.t the observa- tions) of the data distribution. The MSE- clean signal conditioned on the noisy signal. Therefore, Tweedie's 's formula in fact links between the first posterior moment and the score of the data. A similar relation holds between the second posterior moment (i.e., the posterior covariance) and the second-order score (i.e., the Hessian of the prgpobbbabiltyy w.r.t the observations) (Gribonval, 2011), which is in turn associated with the deriva- tive (i.e., Jacobian) of the posterior moment. Recent works used this relation to quantify uncertainty in denoising (Meng et al., 2021), as well to improve score-basee generative models (Dockhorn et al., 2022; Lu et al., 2022; Meng et al., 2021; Mou et al., 2021; Sabanis & Zhang, 2019).
denoiser corresponds to the posterior mean of the
In this paper derive relation between higher-order posterior central moments and higher-order derivatives of the posterior mean in Gaussian denoising. Our result provides a simple mechanism that, given the MSE-optimal denoiser function and its derivatives at some input, allows determining

---

Published as a conference paper at ICLR 2024
the entire posterior distribution of clean signals for that particular noisy input (under mild condi- tions). Additionally, we prove that a similar result holds for the posterior distribution of the projec- tion of the denoised output onto a one-dimensiona. direction.
We leverage our results for uncertainty quantification in Gaussian denoising by employing a pre- trained denoiser. Specifically, we show how results allow computing the top eigenvectors of the posterior covariance (i.e,, the posterior principal components) for any desired region of the image. We further use our results for approximating the entire posterior distribution along each posterior principal direction. As we show, this provides valuable information on the uncertainty in the restoration. Our approach uses only forward passes through the pre-trained denoiser and is thus advantageous over previous uncertainty quantification methods. In particular, it is training- free, fast, memory-efficient, and applicable to high-resolution images. We illustrate our approach with several pre-trained denoisers on multiple domains, showing its practical benefit in uncertainty visualization.
2 RELATED WORK
Many works studied theoretical properties of MSE-optimal denoisers for signals contaminated by additive white Gaussian noise. Perhaps the most well-known result is Tweedie's 's formula (Efron, 2011; Miyasawa et al., 1961; Robbins, 1956; Stein, 1981), which connects the MSE-optimal de- noiser with the score function of noisy signals. Another interesting property, shown by Gribonval (2011), is that the MSE-optimal denoiser can be interpreted as a maximum-a -a-nuuum---ttttioriiii (MAP) es- timator, but with a possibly different prior. The work most closely related to ours is that of Meng et al. (2021), who studied the estimation of high-order scores. Specifically, they derived a rela- tion between the high-order posterior non-central moments in Gaussian denoising task, and the high-order scores of the distribution of noisy signals. They discussed how these relations can be used for learning high-order scores of the data distribution. But due to the large memory cost of storing high-order moment tensors, and the associated computational cost during training and in- ference, they trained only second-order score models and only on small images (up to 32 x 32). They used these models for predicting the posterior covariance in denoising tasks, as well as for improving the mixing speed of Langevin dynamics sampling. Their result is based on a recursive relation, which they derived, between the high-order derivatives of the posterior mean and the high- order non-central moments of the posterior distribution in Gaussian denoising. Specifically, they showed that the non-central posterior moments mi, m2; m3, .. admit a recursion of the form m = f(ms, Vm,, m1).
In many settings, central moments are rather preferred over their non-central counterparts. Indeed, they are more numerically stable and relate more intuitively to uncertainty quantification (being di- rectly linked to variance, skewness, kurtosis, etc.). Unfortunately, the result of (Meng et al., 2021) does not trivially translate into a useful relation for central moments. Specifically, one could use the fact that the ktth central moment, M4x can be expressed in terms of {m, } and that each m; can be written in terms of { But naively substituting these relations into the recursion of Meng et al. (2021) leads to an expression for 14 that includes all lower-order central-moments and their high-order derivatives. Here, we manage to prove a very simple recursive form for the central moments, which takes the form μ [[4,, Vuk, (22). Another key contribution, which we present beyond the framework studied by Meng et al. (2021), relates to marginal posterior distribu- tions along arbitrary cross-sections. Specifically, we prove that the central posterior moments of any low-dimensional projection of the signal, also satisfy similar recursion. Importantly, we show how these relations can serve as very powerful tools for uncertainty quantification denoising tasks.
Uncertainty quantification has drawn significant attention in the context of image restoration. Many works focused on per-pixel uncertainty prediction (Angelopoulos et al., 2022; Gal Ghahramani, 2016; Horwitz & Hoshen, 2022; Meng et al., 2021; Oala et al., 2020), which neglects correlations between the uncertainties of different pixels in the restored image. Recently, several works forayed into more meaningful notions of uncertainty, which allow to reason about semantic variations (Kutiel et al., 2023; Sankaranarayanan et al., 2022). For example, a concurren work by Nehme et al. (2023) presented method for learning the posterior principal components of arbitrary inverse problems. However, all existing methods either require a pre-trained generative model with a disentangled latent space (e.g., StyleGAN (Karras al., 2020) for face images) or, like many of their per-pixel
2

---

Published as a conference paper at ICLR 2024
counterparts, require training. Here present training- computationally efficient, method that only requires access to a pre-trained denoiser.
MAIN THEORETICAL RESULT
We now present our main theoretical result, starting with scalar denoising and then extending the discussion to the multivariate setting. The scalar case serves two purposes. First, it provides intu- ition. But more importantly, the formulae for moments of orders higher than three are different for the univariate and multivariate settings, and therefore the two cases require separate treatment.
3.1 THE UNIVARIATE CASE
Consider the univariate denoising problem corresponding the observation model
y= x n, where x is a scalar random variable with probability density function Ps and the noise n ~ ~N(0, a²) statistically independent of The goal in denoising is provide a prediction of x, which is function of the measurement y. It is well known that the predictor minimizing the MSE, E[(x x)
the posterior mean of given y. Specifically. given particular measurement y = y, the MSE- optimal estimate is the first moment of the posterior density Pxly (-lv), which we denote by
μ1 (y) = E = yl.
(2)
While optimal in the MSE sense, the posterior mean provides very partial knowledge on the possible values that x could take given that y y. More information is encoded in higher-order moments of the posterior. For example, the posterior variance provides measure of uncertainty about the MSE-optimal prediction, the posterior third moment provides knowledge about the skewness of the posterior distribution, and the posterior fourth moment can already reveal bimodal behavior.
Let us denote the higher-order posterior central moments by
(3) Our key result is that knowing the posterior mean function μt() (.) and its derivatives at can be used recursively compute all higher-order posterior central moments at y (see proof in App. A). Theorem (Posterior moments in univariate denoising). In the scalar denoising setting of (1), the high-order posterior central moments of x given y satisfy the recursion
(x(y)= = E [(x 1ik--ua (y))))) yy =y 2
p2(y) a² p1(())
pa(y) = o² pl)()).
(x+1(y) = o² p4 (v) + kpp-
3
(4)
Thus, (3+ is uniquely determined by pı p4 (v), (i) (z) (v)-
Figure illustrates this result via a simple example. Here, the distribution of x is mixture of two Gaussians. The left pane depicts the posterior density (-)y(-1-) as well as the posterior mean function (.).). We focus on the measurement = y shown as a vertical dashed line, for which the posterior Pxly(-ly"") bimodal (right pane). This property cannot be deduced by merely examining the MSE- optimal estimate μn (y"). However, this information does exist in the derivatives of (2(-) at y". To demonstrate this, we numerically differentiated μz (-) at y". used the first three derivatives to extract the first four posterior moments using Theorem and computed the maximum entropy distribution that matches those moments (Botev & Kroese, 201 1). As can be seen, this already provides good approximation of the general shape of the posterior (dashed red line).
Theorem I has several immediate implications. First, it is well known that if the moments do not grow too fast, then they uniquely determine the underlying distribution (Lin, 2017). This is the case e.g., for distributions with a compact support and is thus relevant to images, whose pixel values typically lie in [0, 1]. For such settings, Theorem 1 implies that knowing the posterior mean the neighborhood of some point y. allows determining the entire posterior distribution for that point. second interesting observation, is that Theorem can be evoked to show that the posterior is Gaussian whenever all high-order derivatives of fn (.) vanish (see proof in App. F).
Corollary 1. Assume that pl
(k) (y") for all k 1. Then the posterior Psly (ly")) is Gaussian.
3

---

Published as a conference paper at ICLR 2024
Pelp
Prtp unl//))
0.3
GT
Retimated
4
0
Au (9)
-10
10 -5
10 15
20
Figure I: Recovering posteriors in univariate denoising. The left pane shows the posterior distri- bution (-1) and the posterior mean function (1(-) for the scalar Gaussian denoising task (I). On the right we plot the posterior distribution of x given that y = y², along with an estimate of that dis- tribution, which we obtain by analyzing the denoiser function (11 ))) the vicinity of y². Specifically, this estimate corresponds to the maximum entropy distribution that matches the first four moments, which are obtained from Theorem 1 by numerically approximating μd{ (y²) p A"( A.pp" (y)).
3.2 THE MULTIVARIATE CASE
We move to treat the multivariate denoising problem. Here is random vector taking values in Rd,, the noise n ~ (0, Ia) is white multivariate Gaussian vector that is statistically independent of x, and the noisy observation is
= x n.
(5)
As in the scalar setting, given noisy measurement y. are interested in the posterior distri- bution P (-ly). The MSE-optimal denoiser is, again, the first-order moment this distribution,
i (y) E[x y y),
which is a d dimensional vector. The recond-order central moment is the posterior covariance 22 = Cov(x |y y):
which is d x d matrix whose (i1, i2) entry is given by
(pe(1)) = E (xx,, (y)))), (x,, (11(y))s2)l y)
(8)
d array with k indices (a kth order tensor), whose component at multi-index (f1....,iz. ..) is given by
For any 1 3, the posterior keth-rrder central moment
(u)(y))........
[(x;, [(r(y)) (x1, [xz(r))i y)
As we now show, similarly to the scalar case, having access to the MSE-optimal denoiser and its derivatives, allows to recursively compute all higher order posterior moments (see proof in App. B). Theorem 2 (Posterior moments in multivariate denoising). Consider the multivariate denoising setting of with dimension 2 2. For any 2 and any k k+1 indices i....,,i441 E {1,,-.. (1....,d) d), the high-order posterior central moments of x given y satisfy the recursion
[H2(y))i = (1)r(
dy
[us())liidsiss = 01p2 (y)lew
Oy
Σ [ (x))e (12(1))
1 3,
(10)
Oy
where& e A (di.....dy 11-1, 1741 ..., ix). Thus, Pa+ (g) is uniquely determined by μ (3) and by the derivatives up to order k of its elements with respect to the elements of the vector y.
Note that the first line (10) can compactly written as
3pp (
(2(1) o²
(11)
dy
4

---

Published as a conference paper at ICLR 2024
Algorithm Efficient computation of posterior principal components
Input: N (number of PCs), K (number iterations), (+) (MSE-optimal denoiser), (noisy
input), o² (noise variance), < 1 (linear approx. constant)
Initialize (m)) ( N
o2)
1i=1
2: for k <- 1 to K do
3: for f- 1 to N do
U + Mily + (2) -2((1)
(6)
4:
5:
Q,
R. QR OR.DPOOOPOSITIOX(
1k 1
(M)
6:
k
QQ
Q
7:
(0)
V
K
(i)
8:
Il
+
cu1- (1))l
where (1) denotes the Jacobian of μi at y. This suggests that, in principle, the posterior co- By
variance of an MSE-optimal denoiser could be extracted by computing the Jacobian of the model using e.g., automatic differentiation. However, in settings involving high-resolution images, storing this Jacobian is impractical. In Sec. 4.1, we show how the top eigenvectors of (2 (iee. the posterior principal components) can be computed without having to ever store (2 in memory.
Moments of order greater than two pose an even bigger challenge, as they correspond to higher- order tensors. In fact, even they could somehow be computed, is not clear how they would be visualized in order to communicate the uncertainty of the prediction to a user. A practical solution could be to visualize the posterior distribution of the projection of x onto some meaningful one- dimensional space. For example, one might be interested in the posterior distribution of projected onto one of the principal components of the posterior covariance. The question, however, how to obtain the posterior moments of the projection of onto deterministic d-dimensional vector v.
Let us denote the first posterior moment of o x (i.e., its posterior mean) by μu (y). This moment is given by the projection of the denoiser's output onto v.
# (y) [ x y y)
E fxly y) pa (3).
(12)
Similarly, let denote the kth order posterior central moment of v T
xy
# ( E
(v [x T T (v)) = 2.
(13)
As we show next, the ccalar-valued functions (nE satisfy recursion similar to (4) (see proof in App. C). In Sec. 5, we use this result for uncertainty visualization.
Theorem 3 (Directional posterior moments in multivariate denoising). Let u be deterministic d-dimensional vector. Then the posterior central moments of v are given by the recursion
5 (yy Dup' (3).
16 (y) Duji (y),
p)+1(y) o² Duu (g) + kpb k: 3.
(14)
Here ( denotes the directional derivative of function f RRd - RR in direction v al
4 APPLICATION TO UNCERTAINTY VISUALIZATION
We now discuss the applicability of our results in the context of uncertainty visualization. We start with efficient computation of posterior principal components (PCs), and then illustrate the approximation marginal densities along those directions.
4.1 EFFICIENT COMPUTATION OF POSTERIOR PRINCIPAL COMPONENTS
The top eigenvectors of the posterior covariance, (2( capture the main modes of variation around the MSE-optimal prediction. Thus, as we illustrate below, they reveal meaningful information re- garding the uncertainty of the restoration. Had we had access the matrix 222 computing these
5

---

Published as a conference paper at ICLR 2024
Noisy Image
43 33,70
Figure 2: Computing posterior principal components for pre-trained face demoising model. For each noisy image y. we depict one of the posterior PCs obtained with Alg. 1. To the right that PC, we show the denoiser's output, P (y), and its perturbation along that PC. As can be seen, this visualization captures the denoiser's uncertainty along semantically meaningful directions, such as the color of the moustache, the thickness of the lips, and the extent which the mouth is open.
top eigenvectors could be done using the subspace iteration method (Arbenz, 2016; Saad, 201 1). This technique maintains set of N vectors, which repeatedly multiplied by 22( and orthonor- malized using the QR decomposition. Unfortunately, storing the full covariance matrix is commonly impractical. To circumvent the need for doing so, recall from (11) that (1y corresponds the Jacobian of the denoiser μu (y). Thus, every iteration of the subspace method corresponds to rectti product. For neural denoisers, such products can be calculated using automatic dif- ferentiation (Dockhorn et al., 2022). However, this requires computing backward pass through the model each iteration, which can become computationally demanding for large images!. Instead, propose to use the linear approximation
0u)(y) By
u vo mn + cv) m1(
(15)
which holds for any Id when R is sufficiently small. This allows applying the subspace iteration using only forward passes through the denoiser, as summarized Alg. 1. As we show in App. H, this approximation has negligible effect on the calculated eigenvectors, but leads e,g., to 6x reduction in memory footprint for a 80 92 patch with the SwinIR denoiser (Liang et al., 2021). We note that to compute the PCs for user-chosen region of interest, all that is required is mask out all entries of a outside that region in each iteration.
Figure illustrates this technique in the context of denoising of face images contaminated by white Gaussian noise with standard deviation a = 122. We use the denoiser from (Baranchuk et al., 2022), which was trained part of DDPM model (Ho al., 2020) on the FFHQ dataset (Karras et al., 2019), Note that here we use it as plain denoiser (as used within a single timestep of the
Note that backward passes for whole images also often avoided during training of neural denoisers. Indeed, typical training procedures limited-sized crops.

---

Published as a conference paper at ICLR 2024
12.5
Posterise alnng 15
.25
Poscecir slong
Eximatal
1014
19/4
0.06
7.5
7.5
50-
0004
-2V3
0v7
2v3
2v3
pnly ony)
25
4ol9 ony)
0,01
0.02
-25
0.00
-2/5
0v/7
2v7
-2./
Figure 3: Computing marginals along principal components. On the left, we show the prior Px as a heatmap. a noisy sample y (red), the corresponding MSE-optimal estimate μ4 (7) (black), and the two principal axes, computed using Alg. Here, we used the closed form for (y)) The second pane shows the marginal posterior distribution along the first principal component, computed both using our proposed procedure (dashed red), and by using the closed-form solution (solid blue). On the right we show the same experiment, but with simple neural network trained on data samples.
DDPM). We showcase examples from the CelebAMask-HQ dataset (Lee et al., 2020). As can be seen, different posterior principal components typically capture uncertainty in different localized regions of the image. Note that this approach can be applied any oogite-of-interest within the image, chosen by the user at test time. This is in contrast to a model that is trained to predict a low- rank approximation of the covariance, to the specific input size on which it was trained, and cannot be manipulated at test time to produce eigenvectors corresponding some user-chosen region (cropping a patch from eigenvector is not equivalent to computing the eigenvector of the corresponding patch in the image). In App. K we report quantitative comparisons to the naive baseline of estimating the PCs using a posterior sampler, and quantitatively evaluate the accuracy of the eigenvalues predicted by our method.
in (Meng et al., 2021). Such a model inherently limited
4.2 ESTIMATION OF MARGINAL DISTRIBUTIONS ALONG CHOSEN DIRECTIONS
more fine- fine-graine grained characterization of the posterior can be achieved by using higher-order moments along the principal directions. These can be calculated using Theorem 3, through (high-order) nu- merical differentiation of the one-dimensional function f(a) v. T n1 (y + av)) at a = 0. Once we obtain all moments up to some order, we compute the probability distribution with maximum entropy that fits those moments. In practice, we compute derivatives up to third order, which allows us to obtain all moments up to order four.
Figure illustrates this approach on a two-dimensional Gaussian mixture example with a noise level of a 2. On the left, show heatmap corresponding to Px (-), as well as noisy input y (red point) and its corresponding MSE-optimal estimate (black point). The two axes of the ellipse are the posterior principal components computed using Alg. 1 using numerical differentiation of the closed-form expression of the denoiser (see App. E). The bottom plot on the second pane shows the the largest eigenvector. We numerically computed its derivatives up to order three at a 0 (black point), from which we estimated the moments up to order four according Theorem 3. The top plot on that pane shows the ground-truth posterior distribution of v. x, along with the maximum entropy distribution computed from the moments. The right half of the figure shows the same experiment only with neural network that was trained on pairs of noisy (pink) samples and their clean (blue) counterparts. This denoiser comprises 5 layers with (100, 200, 200, 100) hidden features and SiLU (Hendrycks & Gimpel, 2016) activation units. We trained the network using Adam (Kingma & Ba, 2015) for 300 epochs, with learning rate of 0.005.
function fa) corresponding
Figure 4 illustrates the approach on a handwritten digit from the MNIST (LeCun, 1998) dataset. Here, we train and use simple CNN with 10 layers of 64 channels, separated by ReLU activation layers followed by batch normalization layers. As can be seen, fitting the maximum entropy distri- bution reveals more than just the main modes of variation, as it also reveals the likelihood of each reconstruction along that direction. It is instructive to note that although the two extreme reconstruc- tions, l (y) + VXzt2, look realistic, they are not probable given the noisy observation. This is the
7

---

Published as a conference paper at ICLR 2024
Noisy Image P;
D/2
P3
-3,2,p
9 9
9
9
9
4
4
3
- 2
1 7/
1 Ag
2 2
3 Ag
Figure 4: Uncertainty quantification for denoising a handwritten digit. The first three PCs corresponding to the noisy image are shown on the left. On the right, images along the third PC, marked in blue, are shown, together with the marginal posterior distribution we estimated for that direction. The two modes the possible restoration, corresponding to the digits and 9, can clearly be seen as peaks in the marginal posterior distribution, whereas the MSE-optimal restoration in the middle is obviously less likely.
reason their corresponding estimated posterior density is nearly zero. In App. K, we quantitatively validate the advantage of using higher-order moments for estimating the marginal distribution.
Our theoretical analysis applies to non-blind denoising, in which is known. However, we empiri- cally show in Sec. 5 and Fig. 5 that using an estimated a is also sufficient for obtaining qualitatively plausible results. This can either be obtained from noise estimation method (Chen al., 2015) or even from the naive estimate d² = ( yll? where μr (y) is the output of a blind denoiser. Here we use the latter. We further discuss the impact of using an estimated in App.
5 EXPERIMENTS
We conduct experiments with our proposed approach for uncertainty visualization and marginal posterior distribution estimation on additional real data in multiple domains using different models.
We showcase our method on the MNIST dataset, natural images, human faces, and on images from the microscopy domain. For natural images, we use SwinIR (Liang et al., 2021) that was pre-trained on 800 DIV2K (Agustsson & Timofte, 2017) images, 2650 Flickr2k (Lim et al., 2017) images, 400 BSD500 (Arbelaez et al., 2010) images and 4,744 WED (Ma al., 2016) images, with patch sizes 128 128 and window size 8 x 8. We experiment with two SwinIR models, trained separately for noise levels o = {25. 50}. and showcase examples on test images from the CBSD68 (Martin et al., 2001) and Kodak (Franzen, 1999) datasets. For the medical and microscopy domain we use Noise2Void (Krull et al., 2019), trained and tested for blind-denoising on the FMD dataset (Zhang al., 2019) in the unsupervised manner described by Krull et al. (2020). The FMD dataset was collected using real microscopy imaging, and as such its noise is most probably not precisely white Gaussian, and the noise level is unknown in essence (the ground truth images are considered as the average of 50 burst images). Accordingly, N2V is a blind-denoiser, and we have no access to the "real" o. therefore, for this dataset we used an estimated in our method, described in Sec. 4.2.
Examples for the different domains can be seen in Figs. 2, 4, and 5. As can be seen, in all cases, our approach captures interesting uncertainty directions. For natural images, those include cracks, wrinkles, eye colors, stripe shapes, etc. In the biological domain, visualizations reveal uncertainty in the size and morphology of cells, as well as in the (in)existence of septums. Those constitute important geometric features in cellular analysis. More examples can be found in App. I.
One limitation of the proposed method is that it relies on high-order numerical differentiation. As this approximation can be unstable with low precision computation, we use double precision during the forward pass the networks. Another method that can be used to mitigate this is to fit a low 0, and then
degree polynomial to f(a) v uu + av) around the point of derivation, use the smooth polynomial fit for the high-order derivatives calculation. Empirically we found the on the choice of the polynomial degree and the fitted range. This caused bad fits even for the simple two-component GMM example, whereas the numerical derivatives approximations worked better.
polynomial fitting also be sensitive,
8

---

Published as a conference paper at ICLR 2024
Noisy Image V
V2
3/3
-3,20
+3/2DD
Natural Images (SwinIR)
Miczoso00 Images (N2V)
3,7 -2 1
1,7
2/44
3/74
Figure 5: Uncertainty quantification for natural image denoising using SwinIR (top) and mi- croscopy image denoising using N2V (bottom). In each row, the first three PCs corresponding to the noisy image are shown on the left, and one is marked in blue. On the right, images along marked PC are shown above marginal posterior distribution estimated for this direction. The PCs show the uncertainty along meaningful directions, such as the existence of cracks on an old vase and changes in the tiger's stripes, as well as the sizes of cells and the existence of septum, which constitute important geometric features in cellular analysis.
CONCLUSION
Denoisers constitute fundamental ingredients in variety of problems. this paper we derived relation in the denoising problem between higher-order derivatives of the posterior mean to higher- order posterior central moments. These results were then used in the application of uncertainty visualisation of pre-trained denoisers. Specifically, we proposed method for efficiently computing the principal components of the posterior distribution, in any chosen region of an image. Addition- ally, we presented scheme to use higher-order moments to estimate the full marginal distribution along any one-dimensional direction. Finally, we demonstrated our method on multiple denoisers across different domains. Our method allows examining semantic directions of uncertainty by using only pre-trained denoisers. in fast and memory-efficin way. While the theoretical basis of our method applies only to additive white Gaussian noise, we show empirically that our method provides qualitatively satisfactory results also in blind denoising real-world microscopy data.

---

Published as a conference paper at ICLR 2024
REPRODUCIBILITY STATEMENT
As part of the ongoing effort to make the field of deep learning more reproducible and open, we pub- lish our code https: / / /hi Lamanor... github...........g.. /oo//aussae anDenoia wwwisssinnooeiiiooooooooooo The repository includes scripts to regenerate all figures. Researchers that want to re- re-implement the code from scratch can use Alg. 1 and our published code as guidelines. and detailed proofs for all claims in the paper in Appendices A, B, C, E, and F of the supplementary material. Finally, we provide Appendix a translation from our notation to the notation of Meng et al. (2021) to allow future researchers to use both methods conveniently.
addition, we provide full
ETHICS STATEMENT
In many scientific and medical domains, signals contaminated by noise, and deep learning based denoising models have emerged as popular tools for restoring such low-fidelity data. However, de- noising problems inherently ill-posed. Therefore, a system that presents users with only a single restored signal, may mislead the data-analyst, researcher, or physician into making flawed decisions. To avoid such situations, is of utmost importance for systems to also report and conveniently vi- sualize the uncertainties in their predictions. Such systems would be much more trustworthy and interpretable, and will thus support making credible deductions and decisions. The method we pre- sented in this paper, can help visualize the uncertainty in denoiser's 's prediction, by allowing users to explore the dominant modes of possible variations around that prediction, accompanied by their likelihood (given the noisy measurements). Such interactive denoising systems, would allow users to take into consideration other, and sometimes even more likely possibilities than e.g., the minimum MSE reconstruction that is often reported as single solution.
ACKNOWLEDGEMENTS
The Miriam and Aaron Gutwirth Memorial Fellowship supported the research of HM. The research of TM was supported by the Israel Science Foundation (grant no. 2318/22), by the Ollendorff Miverva Center, ECE faculty, Technion, and by gift from Elbit. The authors are grateful to Elias Nehme, Rotem Mulayoff, and Matan Kleiner for their insightful discussions and input throughout this work, and Noa Cohen for her invaluable help with the algorithm.
REFERENCES
Eirikur Agustsson and Radu Timofte. NTIRE 2017 challenge on single image super-resolution: Dataset and study. In Proceedings of the IEEE conference computer vision and pattern recog- nition workshops, pp. 126- 135, 2017. URL https: http://daaaaaavi.eeeeee. DTV2K/.
vision. wewionn.eecee ww. ch/ cvl/
Michal Aharon, Michael Elad, and Alfred Bruckstein. K-SVD: An algorithm for designing over- complete dictionaries for sparse representation. IEEE Transactions signal processing, 54(11): 4311-4322, 2006.
Anastasios N Angelopoulos, Amit Pal Kohli, Stephen Bates, Michael Jordan, Jitendra Malik, Thayer Alshaabi, Srigokul Upadhyayula, and Yaniv Romano. Image-to-image regression with distribution-free uncertainty quantification and applications in imaging. In International Confer- ence on Machine Learning, pp. 717-730. PMLR, 2022.
Pablo Arbelaez, Michael Maire. Charless Fowlkes, and Jitendra Malik. Contour detection and hierarchical image segmentation. IEEE transactions on pattern analysis and machine intelli- gence, (3(5):898-916 2010. URL https: http://ww//////c/////////w eecs berkeeereeeee...e edu Reseeache Projects CS/viss/issss s sn ggrogrroppoooooo ng/ resourcesssseeee html.
Peter Arbenz. Lecture notes on solving large scale eigenvalue problems. 2016.
Dmitry Baranchuk, Andrey Voynov, Ivan Rubachev, Valentin Khrulkov, and Artem Babenko. Label- efficient semantic segmentation with diffusion models. In International Conference on Learning Representations, 2022.
10

---

Published as a conference paper at ICLR 2024
Zdravko I Botev and Dirk P Kroese. The generalized cross entropy method, with applications to probability density estimation. Methodology and Computing in Applied Probability, 13(1):1-27, 2011.
Alon Brifman, Yaniv Romano, and Michael Elad. Turning denoiser into super-resolver using plug and play priors. In 2016 IEEE International Conference on Image Processing (ICIP), pp. 1404-1408. IEEE, 2016.
Antoni Buades, Bartomeu Coll, and J-M Morel. non-local algorithm for image denoising. In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR'OH 005, volume 2, pp. 60-65. IEEE, 2005.
Guangyong Chen, Fengyuan Zhu, and Pheng Ann Heng. An efficient statistical method for image noise level estimation. In Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV), pp. 477-485, 2015.
Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian. Image denoising by sparse 3-d transform-domain collaborative filtering. IEEE Transactions on image processing, 16 (8):2080-2095, 2007.
Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large-scale hi- erarchical image database. In 2009 IEEE conference on computer vision and pattern recognition, pp. 248-255. leee, 2009.
Tim Dockhorn, Arash Vahdat, and Karsten Kreis. GENIE: Higher-order denoising diffusion solvers. In Alice H. Oh, Alckh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022.
Bradley Efron. Tweedie's formula and selection bias. Journal of the American Statistical Associa- fion, 106(496):160226111)) 2011.
Rich Franzen. Kodak lossless true color image suite. source: http://rOk: us/graphics/kodak 1999.
Yarin Gal and Zoubin Ghahramani. Dropout as bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning, pp. 1050-1059.
PMLR, 2016.
Rémi Gribonval. Should penalized least squares regression be interpreted maximum posteriori estimation? IEEE Transactions on Signal Processing, 59(5):2405-2410 2011.
Dan Hendrycks and Kevin Gimpel. Gaussian error linear units (gelus). arXiv preprint
arXiv:1605.084155 2016.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6855 2020.
Eliahu Horwitz and Yedid Hoshen. Conffusion: Confidence intervals for diffusion models. arXiv preprint arXiv:2211 1.099.9 2022.
Tero Karras, Samuli Laine, and Timo Aila. A style-based generator architecture for genera- tive adversarial networks. In Proceedings of the IEEE/CVF conference computer vision and pattern recognition, pp. 4401-4410, 2019. URL htttps://githu..... cot.uull.co////llbb...... ffhq dataset.
Tero Karras, Miika Aittala, Janne Hellsten, Samuli Laine, Jaakko Lehtinen, and Timo Aila. Training generative adversarial networks with limited data. Advances in neural information processing systems,
2020.
Asem Khmag, Abd Rahman Ramli, SA Al-Haddad, and Noraziahtulhidayu Kamarudin. Natural im- age noise level estimation based on local statistics for blind noise reduction. The Visual Computer: International Journal of Computer Graphics, 34(4):575-587, 2018.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International
Conference Learning Representations, 3, 2015.

---

Published as a conference paper at ICLR 2024
Priyanka Kokil and Turimerla Pratap. Additive white gaussian noise level estimation for natural images using linear scale-space features. Circuits, Systems, and Signal Processing. 40(1):353- 374, 2021.
Alexander Krull, Tim-Oliver Buchholz, and Florian Jug. Saase2ooid-eeaninng denoising from sin- gle noisy images. In Proceedings of the IEEE/CVF conference computer vision and pattern recognition, pp. 2129-2137, 2019.
Alexander Krull, Tomáš Vičar, Mangal Prakash, Manan Lalit, and Florian Jug. Probabilistic noise2void: Unsupervised content-aware denoising. Frontiers in Computer Science, 2:5, 2020.
Gilad Kutiel, Regev Cohen, Michael Elad, Daniel Freedman, and Ehud Rivlin. Conformal prediction masks: Visualizing uncertainty in medical imaging. In ICLR 2023 Workshop on Trustworthy Machine Learning for Healthcare, 2023.
Yan LeCun. The MNIST database handwritten digits. 1998. URL http:://apat Lecun.
com/ hon//xxdb/mag ticb//mnist/
Cheng-Han Lee, Ziwei Liu, Lingyun Wu, and Ping Luo. Maskgan: Towards diverse and interactive facial image manipulation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5549-5558, 2020.
Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte. Swinir: Image restoration using swin transformer. In Proceedings of the IEEE/CVF international confer- ence on computer vision, pp. 1833-1844, 2021.
Bee Lim, Sanghyun Son, Heewon Kim, Seungjun Nah, and Kyoung Mu Lee. Enhanced deep resid- ual networks for single image super-resolution. In The IEEE Conference Computer Vision and Pattern Recognition (CVPR) Workshops, July 2017. URL https: //tptgggtttttt///t////g/ com/ limbee/ NTIRE2017. 7.
Gwo Dong Lin. Recent developments on the moment problem. Journal of Statistical Distributions and Applications, 4(1):5, 2017.
Wei Liu and Weisi Lin. Additive white gaussian noise level estimation in svd domain for images. IEEE Transactions Image processing, 22(3):872 883, 2012.
Xinhao Liu, Masayuki Tanaka, and Masatoshi Okutomi. Single-image noise level estimation for blind denoising. IEEE transactions image processing, 22(12):5226-5237
Cheng Lu, Kaiwen Zheng, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Maximum likelihood training for score-based diffusion odes by high order denoising score matching. In International Conference on Machine Learning, pp. 14429-14460, PMLR, 2022.
Kede Ma, Zhengfang Duanmu, Qingbo Wu, Zhou Wang, Hongwei Yong, Hongliang Li, and Lei Zhang. Waterloo exploration database: New challenges for image quality assessment mod- els. IEEE Transactions on Image Processing, 1604)1000 2016. URL https://ece. uwaterloo. ca/caa//// aa/-k2/ k2 9ma/ exploration/.
David Martin, Charless Fowlkes, Tal D., and Jitendra Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecologi- cal statisties. In Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 2, pp. 413-423. IEEE, 2001. URL httpt:/ http::::://ttttb.ooo// c ausmi chele/ -ast6eddatase....
Chenlin Meng, Yang Song, Wenzhe Li, and Stefano Ermon. Estimating high order gradients of the data distribution by demoising. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp- 25359-25366 Curran Associates, Inc., 2021.
Koichi Miyasawa et al. An empirical bayes estimator of the mean of a normal population. Bull. Inst. Internat. Statist, 38(181-188))1-2 1961.
12

---

Published as a conference paper at ICLR 2024
Wenlong Mou, Yi-An Ma, Martin Wainwright, Peter L Bartlett, and Michael Jordan, High- order langevin diffusion yields an accelerated MCMC algorithm. Journal of Machine Learning Research, 22(42):1-41, 2021.
2021.
Elias Nehme, Omer Yair, and Tomer Michaeli. Uncertainty quantification via neural posterior princi- pal components. In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
Luis Oala, Cosmas HeiB, Jan Macdonald, Maximilian März, Wojciech Samek, and Gitta Kutyniok. Interval neural networks: Uncertainty scores, arXiv preprint arXiv:2003. 1166, 2020.
Javier Portilla, Vasily Strela, Martin J Wainwright, and Eero P Simoncelli. Image denoising using scale mixtures of gaussians in the wavelet domain. IEEE Transactions Image processing, 12
(11):1338-1351, 2003.
Herbert Robbins. An empirical bayes approach to statistics. In Proceedings of the Third Berke- ley Symposium on Mathematical Statistics and Probability, 1954- 1955, volume 1, pp. 157-163. Berkeley and Los Angeles: University of California Press, 1956.
Yaniv Romano, Michael Elad, and Peyman Milanfar. The little engine that could: Regularization by denoising (RED). SIAM Journal on Imaging Sciences, (0(4):1804-18 2017.
Stefan Roth and Michael J Black. Fields of experts. International Journal of Computer Vision, 82: 205-229, 2009.
Leonid Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259-263 1992.
Yousef Saad. Numerical methods for large eigenvalue problems: revised edition. SIAM, 2011.
Sotirios Sabanis and Ying Zhang. Higher order langevin monte carlo algorithm. Electronic Journal of Statistics, 13(2):3805-3850, 2019.
Swami Sankaranarayanan, Anastasios Nikolas Angelopoulos, Stephen Bates, Yaniv Romano, and Phillip Isola. Semantic uncertainty intervals for disentangled latent spaces. Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learn- ing. pp. 2256-2265. PMLR, 2015.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Charles M Stein. Estimation of the mean of a multivariate normal distribution. The Annals of Statistics, 96)):1135-1151 1981.
Tom Tirer and Raja Giryes. Image restoration by iterative denoising and backward projections. IEEE Transactions on Image Processing, 28(3):1220-1233 2018.
Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In 2013 IEEE Global Conference on Signal and Information Processing, pp. 945-948. IEEE, 2013.
Pascal Vincent. connection between score matching and denoising autoencoders. Neural compu- tation, 23(7):1661-166 2011.
Yinhuai Wang, Jiwen Yu, and Jian Zhang, Zero-shot image restoration using denoising diffusion null-space model. In The Eleventh International Conference on Learning Representations, 2023.
URL httph:tp////c....t hwtp//peeeeeeeeeeeeeeeeeeeeeeeee review net / forum? i d=mR: eQgMtNTQ. 1
Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang. Beyond gaussian denoiser: Residual learning of deep cnn for image denoising. IEEE transactions on image processing, 26 (7):1142-311 20.
13

---

Published as a conference paper at ICLR 2024
Kai Zhang, Wangmeng Zuo, Shuhang Gu, and Lei Zhang. Learning deep CNN denoiser prior for image restoration. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3929-3938, 2017b.
Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte. Plug-and-play image restoration with deep denoiser prior. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):6300633 6376, 2021.
Lei Zhang. Xiaolin Wu, Antoni Buades, and Xin Li. Color demosaicking by local directional in- terpolation and nonlocal adaptive thresholding. Journal of Electronic Imaging, 20(2):023016, 2011.
Yide Zhang, Yinhao Zhu, Evan Nichols, Qingfei Wang, Siyuan Zhang, Cody Smith, and Scott Howard. A poisson-gaussian denoising dataset with real fluorescence microscopy images. In Pro- ceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 11710- 11718, 2019. URL https: // //gtttut...g//gg/////////g...h.t wwwm//ytraaao/eeeoiiiiiii fluorescence.
14

---

Published as a conference paper at ICLR 2024
SUPPLEMENTARY MATERIAL
A PROOF OF THEOREM
We start with the case & 2 2 (bottom two lines in (4)). In this case, the conditional moment (4) can be expressed using Bayes' formula as
(x(() E [x (1(y))* y vl
/ (o (1(1)) Paly Paly)))PPly( (2ly)dr
[(x (1)) Pyxx))
Py (y(
f(x (2To²)-3 3/2 (1())* expl 200 (v - x)2)px(x)d
(S1)
Py(1)
Denoting the numerator by q(y) & (2 s(x - #1 (v)) 22z (y x)) we can write the derivative of (x) as
14(y)) (())) a(v))(())
PV(())
(((v)
((v) (1(
Py(y)
Py(v) (y())
d'(w)
(M)
py(1)
(T)
d'(v)
dlogp,(v)
(v)
Py(1)
dy
('(v)
1
σ2 (u(l)((aa() -y),
(S2)
(v)
where we used the fact that dy (1)) y) (see e.g., (Efron, 2011; Miyasawa et al., dy
1961: Stein, 1981)). The first term in this expression is given by
('(2) (2mo?) + f dy [(x (1)))* expl ( - m) (x())
Py(v)
(v(
(2 1 / - (1))=)(())))~-1111 (1( - (x - (1(()) ( - x)) exp( 2 (u x))ppx(>) (x))da
Py(y)
I (-k(e -n())) - (x - (1())* da (y x)) pyla (2)
(v)
/ - k(x - (1()) - (x un (x))** 1
σ² (y - x²) (oy(al/))de
= E ( - (x - (y)) & 0 (y x) y = v
(S3)
To allow unified treatment of the cases l = 2 and l: > 2, let us denote
th (v) e E - (y))) y = m] = 00
k= 1,
(ua(v)
(S4)
k 2.
15

---

Published as a conference paper at ICLR 2024
We therefore have
q' ('(y) ko 11)=)1())1(0)
1
1
x² Ux(uly
+
E [(x E (x x1)=)))>//=>>))>> yy= vl
Py(y)
0²
1
σ² vx(u)v+ + = -kbb-1 (v)
-kph - (v)
El(x - (1(1)) - x2(x)) (n()))ly vl
+
(Uk+1(1) - 2x (1)
0²
=-) (v)
x² vx+1(() 02 (x(()) (MI(8) v)
(S5)
Substituting this back into (S2), we obtain that
(v) = kee ( (3) + x2 1x+1
+ 02 (x( (MI(1) - y) (u)( (41(1) - y)
1
=- (y) +
x² (k+1
(S6)
where we used the fact that (y) (x(y) for all 2. Now, for k 2 this equation reads
1 v2((1))
12(1)) =
(S7)
and for k 2 3, it reads
1 vk(1)) = Kun (3))-100) + (e))
(S8)
We thus have that
pa(y) = 0² (d/(1);
(x+1 o² (1) + ( 2 3.
(S9)
Note that an equivalent expression for the last line is obtained by replacing (1) with u2(y)) as we prove below. This completes the proof for k 2 2.
The case h = can be treated similarly. Here,
μi (v) = Fx y y)
(2TO?)-
1/2 ʃ I expl 2t x)² Pxx
(S10)
Py(1)
so that we define ((v) A (2 - f
x 1xeppf 1
202 (y - (ux))²]ppl0 (x)))p((x)dd. We thus have
(2m0²)-3 1 [xexp[ - (u 2 px (x()ddr
q (v) (y(y)
Py (y
(2TT²?) 1 2/ (x y) exp( 2
1 (y 2))ppx(0 (x((x)dz
(y())
E[x(x y) ly = vl
[E|x2 yy = y) - --(())))
(S11)
σ²
Therefore,
('(y) 1
(v(2)=
(1)())(110)) v)
(v(
1
ly = u] - _ (o) - -))
(B[x²
1
(B|x2 ly = vl - (z()))
1
(B(x2ly 1 = v) E(x v)
1
#22(w)
(S12)
16

---

Published as a conference paper at ICLR 2024
which demonstrates that
Nz(v) = (0).
(S13)
This completes the proof for k = 1.
B PROOF OF THEOREM 2
We begin with the case A = 1 (first line in (10)), by directly deriving the matrix form (11), Using Bayes' formula, the posterior mean μa (3) can be expressed
41 (v) = E(xly = yl
xpx y (x
RR
Jaa mp Ix ( px
Py
(2x02) $ JR. x expl 20 l1y 2=12 Jpx (æ)d
(S14)
Py (y
Therefore, denoting the numerator by ((v)
(2Te2) & Jae 2 exppp 255 11 2112) 1 (2)) we can
write the Jacobian of l at y as
0pl(y)
Ogly) Dy
By Py( (19) (Vpy (2))
dy
0q(y) By
((v)
(VDy(t))
Py (y)
(y)
Py(1)
8g(y)
(Vpy (y))
0y
(1)
Py (y)
Py(
8g(y)
(y) (V log (3))
Py (1)
Dg(y)
0y
2 (1) (1) -y
(S15)
Py
Here, Ogly) E dRd denotes the Jacobian of Rd + Rd at y. and we used the fact that v 1l 1 (x) 3) (Efron, 2011: Miyasawa et al., 1961: Stein, 1981). The first term in (S15) can be further simplified as
0g(w) By
(2mo2) Sgr x excpp 1 11y x 112 3 (x 8) px(xx)
Py (
(y(
Jo # xx y) (x))
Py (v)
x(x x - y) Pxly (ar (
Rd 2
1
EE] Xxx T y y) Elxly =y]
2
[E[xx] ly y) nu (V) y
(S16)
22
17

---

Published as a conference paper at ICLR 2024
Substituting (S16) back into (S15), we obtain
0p1 (v)
1
[ Exxx y y) 1---1 (1) y) ( (x -
v))
dy
1
(E(xx' y = yl (1(1))1(()
x²
1
(E(xx)) ly y) Eixly B/E y) ))
0²
1
Cov(xly y)
1
d2 H2(()
(S17)
This completes the proof for k 1.
We now move on to the cases = and 2 3 (second and third lines in (10)). Element (f1,...,44) of the posterior kth order central moment can be expressed
(1))
[xx;, (y))s,) ... (XE, [p) y = y (2mo2) Si (00,, - [M1 (1)),, ... (Di - (un(3)lle)) exxp 23 11 - x]1²}pp (xx((x))e
(y(
((y)
(S18)
Py Py()
where (3))
1 (2mo2) 1 SRS (æi, - INn (B))s,) ... (a,, - [M (y))li))) expl- 2nh2 11y - x } (xxx(()) ddx.
Therefore, for any i&t1 fl,....dd d}. the derivative of [xc: (3) (i......... with respect yiats given
by
Ml (1 i
89
0g(v)
Py
(3))
9(())
BMik
Bp.
Oyin
(²(1))
0g(v)
Opy(w)
a(y)
SMiz+1
Py (y)
Py Py(
0q(w)
Ö Ö (3)
[ee (Mx())) Ak
Py (y)
Dytzt
ag(y)
G2 [ua(t))
(1I( Mixii)
(S19)
(y)
where in the last line we used the fact that V log Py (y) (1r (y) y) (Efron, 2011; Miyasawa et al., 1961; Stein, 1981). The first term here can be written as
0g(y)
J Byii4s [(xi, [H)(())),, -.. (x., [Mr (w)) )xpp expl 20 11l x=|?]] px px(x)) dx
(2
Py
(y
alui(Wilis (xi1 [Ke (1))aa)))) (wo [Mn (M)) expl 20 l1y x1121pxx (2)
f
OVicii
(2TO2)&pyll (1)
almi(u)lie expl 2/2 11 x 1px (2) (x)xx Oehe
(xii, (3))le)) ... (zit-1 [M4] (1)) l
(2TTO2) 2) py(())
+ (ax,, -14(()))s) (ar - (M)) kz) 33 (I - Miemn) expl 1 lly x 1| P da (2 (y
(S20)
18

---

This demonstrates that Published as a conference paper at ICLR 2024
Let treat the cases k and k > separately. When A = 2, the above expres- sion contains precisely three terms, but the first two vanish. Indeed, the first term reduces to (w)lu ([ (y)),,) 0 and the
(w)l Byig E[x;, (yı(y) y yl (f E[x,, - [pl (y)),, ly = yl (لw ([x1(v]) ( = 0.
B9is
ond term to 8yis Therefore, when
8yy; we left only with the last term, which simplifies
On(v)
By's
[xx, (xi(y));,) (x1, - (x)()))))) (Xi, -va) y y)
Py Py(
o²
1
E
[xxi, [mi(r));,) (x1² (x)(r));, X 1y y) o² (22(
x²
1 [xx - (y))))) (xi² (x) (xi, (11(y)) + (y)))s ly y)
0²
0a (22
[ss (43))))i,, 0² (x 0² (22(
(S21)
σ² Substituting this back into (S19), we obtain 01ux (y)lenne
1 न ( 2 (22(1))i
1
σ (M3()))i..2n +
o² (22(v))en (1ur - Vis))
o² | [1)
(S22)
[ma(y)linnais o² dyints
(S23)
which completes the proof for k 2.
When l 2 3, none of the terms in (S20) vanish, and the expression reads Bo(y)
09
0(1z(1)l (2))
0lpe (1)lie)
(14x-(())
+... (1)=1())),, (
(y
dy
dyknts
(1)
Vi--1 = [(x;; (y))))) (Xi, -(m(()))i)) X ly= y)
σ²
( 1
(1x-1(3)1e
σ² (px()))) Bieti
dy
j=1
1
E [[xi, - (m (y))())) (xix - -un (Xie1s - [mI(y)) + (()())) ly y]
Σ ( (1 1 j=1
k
1 Leax ( Dints + oo (ME4a
Syan 02
1 σ² |e4 (y)): (
0lui
1 + oz (1
(m)-1
dym-.
(y)) (()(())) - Virrz):
(S24)
where we used the definition E, (b1...., 15-1, ..... 4x). Substituting this expression back into (S19), we obtain
b [ (v)le, Olmi(y)ls + o² [ (1)))
(S25)
1=1
19

---

Published as a conference paper at ICLR 2024
This demonstrates that
= o² Tlpe (y)]esmmne in
A
alur
[ (y)las
o²
Σ
| (y)le,
dyliut
dyin-r
/=1
3 Jee (y) limmii o²
[ (y)lee (
(S26)
j=1
where we used (S17). This completes the proof for k 2 3.
C PROOF OF THEOREM 3
We will use the fact that for any k 2 the posterior kth order central moment of oTx can be written explicitly by expanding brackets as
K Σ vs [x (n(e))l, y =
E
(v ( (x - (H)((()) y w) =E
Un ... Ui 18 ((xi, (z)))))))..... (x. (mi()))
1=1
d
Σ
Ui
Ui [ (y))
(S27)
111
Let us start with the second moment. From (S27), it is given by
(3( Σ
U1 2/ (
1=11==
= v (z(Y)
(1)
~²o
dy
,²Vn (v (1(1))
0² Du (TT1((y))
DD (1)
(S28)
This proves the first line of (14). Next, we derive the third moment. From (S27), it is given by
(g( =
Vi Ui Vis
2 1l y 1n
2
0i Uig 1
ay
11=1 1==
i (v T (2())o)
Vis
By
20 o Vy (v (2(2)
o² D. (v (2(2)))
(S29)
where in the last line used (S28). This proves the second line of (14).
20

---

Published as a conference paper at ICLR 2024
Finally, we derive the (k: + 1)th moment for any A 2 3. From (S27), it is given by
d
(y)
M
Vi Vinti (
#++1=1
o² (y)plyy)ll..... + Σ
[a (wi)le, (eee())) (1))l
Ui
Ui2t1
Oyliet
Die+i dyikii
Vi Via [uN (u(Y)lle..ng +
$
M1=1
==1
Σ
Vi ... Vi,- Ubzts ... Vie [ (y)le, Σ Σ vyus ( inyiies
1=12
(v)
0²
Uists 0y
Σ
1h (
1++1=1
o'v' (y) kp - (1)
D kpp (1) 1² (z).
(S30)
where in the second line we used (10). This completes the proof of the third line of (14).
D RELATED WORK: ESTIMATION OF HIGHER ORDER SCORES BY DENOISING
The work most related ours is that of Meng et al. (2021). Here, present their results while translating to our notation. Given probability density Py over Rd they define the kth order score syy as the tensor whose entry multi-index (i1, -2..--.. ik) is
(sx(
log py (y),
(S31)
dy: dy
By
for every i,...., in E [I....,dj*. Using our notation, and under the assumption (5) that y is a noisy version of x Px. the denoising score matching method estimates the first-order score (y), which is simply the gradient the log probability, Vy log py P3). This is done by using Tweedie's formula, which links s with the first posterior moment (the MSE-optimal denoiser)
i (y) E[xly y] y o² s1))
(S32)
As noted by Meng et al. (2021), a similar relation links the second-order score with the second posterior moment (i.e,, the posterior covariance) as
(2( Cov(xly y) = ( + I
(S33)
Note from (S31) that S2 is the Hessian of the log-probability. Vý log (yg)y(y), equivalently the Jacobian of the gradient of the log probability, By Vy log (y). And since we have from (S32) that V log (y) 3 (M) y). Eq. (833) be equivalently written
à (12(1)) 0 oy,, o²
(mll) - y
+ + 2 3 1 (v))
(S34)
dy
This illustrates that the second-order formula of Meng et al. (2021) equivalent to (10).
Moving on to higher-order moments, following our notations, Lemma of Meng et al. (2021) states that
E
xly =y] o0 Oy
E E = y] " = v ()())((()) y )
vk 2 1, (S35)
21

---

Published as a conference paper at ICLR 2024
where @ x Rd
denotes k-fold tensor multiplication. This lemma is used in Theorem 3 of Meng et al. (2021), to derive a recursion relating higher-order moments and scores. Substitut- ing (S32), this relation be written
a ay
B x = y o²
E
x y) E y) (1
4k 2
(S36)
Denoting the non-central posterior moment of order l by mx Eq. (S36) can be written com- pactly
mx+ mx+1 (1) 0² dy m. (y) + mx(y) (1))): Vk 2
(S37)
Writing out the elements of m2++1 (y) explicitly, this relation reads
0 mm ( (y))imm + (m))))))s [m=+1 (y)le 0²
(y))
Vk 2
(S38)
dy
It is interesting compare this expression with the recursion for the central moments in Theorem 2. We see that the non-central moments satisfy a sort of one-step recursion (if we disregard the depen- dence on (1) ). in the sense that m depends only on myg In contrast, as be seen in Theorem 2, the central moments satisfy sort of two-step recursion (if we disregard the dependence on (42). in the sense that Mk+1 (z) depends on both μu (y) and (
E POSTERIOR DISTRIBUTION FOR GAUSSIAN MIXTURE PRIOR
In Fig. and Fig. demonstrated our approach one-dimensional and two-dimensional Gaus- sian mixtures, respectively. In both cases, we showed plots of the marginal posterior distribution in the direction of the first posterior principal component, as well the posterior mean for particular noisy input sample. Those simulations relied on the closed-form expressions of the posterior dis- tribution and the marginal posterior distribution along some direction for Gaussian mixture prior. In addition, Fig. and Fig. 3 also contain the maximum entropy distribution estimated using our method, which uses the numerical derivatives of the posterior mean. Here as well we used the nu- merical derivatives of the posterior mean function, which we computed in closed-from. We now present these closed-form expressions for completeness.
Suppose Px is a mixture of L Gaussians,
L px(x) = Σ TWN TNN Ex).
(S39)
Let be random variable taking values in (1...., L} with probabilities T1,.... TL. Then we can think of x as drawn from the ěth Gaussian conditioned on the event that e = l. Therefore,
Pxly (xly)
Pxly. ( (Ply)()))
Σ (xly. 8) Pylc(y) (l)
Py
Σ Nix me, Ee)
Pe
(S40)
where we denoted
A. Nyymm Σ, Ei+0* o2)),
m, Σ, + +010)=100000 m:) mi
Σ Σ, EE +021)
(S41)
22

---

Published as a conference paper at ICLR 2024
As for the marginal posterior distribution along some direction v, it is easy to show that
Po XoTxly((u (ay(())))
Σ Po (0 (f
1=1
Σ Po (aly. e) py Pyic( (()
Py (
Σ va me, [ Eçe)
PL
(S42)
Σ 6'=1 Pe n
F PROOF OF COROLLARY 1
We start by reminding the reader of (4)
(2(y) = 0² (1))
ps(u) = o² (á(();
M (v) = 0² p/ (v) + (v)
k k 3.
We will prove by complete induction that
(m) 0 for all k 2 2 and m 2 1. Pk
(S43)
Note that since for any m 2 2 we have (mm) (v") = 0. for any m 2 we have (m)
Base
11 (v") = 0²H (m++1) 1) (v²)
0
(m) (y*) = o² pi! (mmmt (v²)
0²Mn ui (11+2) (v¹)
0
(1m) (5¹) 0² Uy
(y²)) + fm+t1)(
8ym (ng(w))
8mh
44
N=y"
(1) o Ug
2 (m+1)
(y) + 2
m
(m
(1)))²
10 (2')
U12
1=0
02/16mm (gr) 3 112 (z)))))))) + .... 112 222 (mm) (1'))
(m)
a² P/3 (m+ (z²)
0,
(S44)
where (1) results from the general Leibniz rule.
(m) Induction Assume that pnn (y/) = 0 for all 4 S n
k + and m 2 1. Then,
mm
( (v²) = Tym (o² (4(()) + (1)=)19)
y-
gm
2² u
(m+ 1) (8) k Jy
(p=-1())
y=y"
11) 0 Hk
(m+1) 1) (y") + k
(m) ( -
(m-
(v
(v*) (v")
1=0
2 Us
(m+1) 1)
(y/ + k
(m)
(y )22(²) +... + kp 1 (1) puz (3)")
0,
(S45)
23

---

Published as a conference paper at ICLR 2024
where for (1) the general Leibniz rule was used again, and in (2) used our induction assumption.
This concludes the induction.
Using (S43) therefore obtain for all lo 2 3 that
Mx+1(y") = kpp- p2()",
k-k: (1 "(x")(")))
1-lk 2] 4)p2( (u==5"
A!llyz
111 (V*)
& is odd,
(S46)
k is even.
Since My(y") = o2 uz(")) = 0 as well, the posterior moments are the same as those of a Gaussian distribution. Indeed, the central moments of random variable N ([[2z) o²) are given by
[[(2 - -|)) fo (d - 1)!! d is even, { 0
(S47)
d is odd.
To conclude the proof, all that remains show is moment-deterninacy (i.e., that the sequence of moments uniquely determines the distribution). This is the case, since the moments of a Gaussian distribution are trivially verified satisfy e.g., Condition (h6) of (Lin, 2017). This implies that the posterior is moment-determinet and is Gaussian.
G EXPERIMENTAL DETAILS
Algorithm requires three hyper- hyper-paranene parameters input. The first is the small constant which is used for the linear approximation (15). The second is N, which is the number of principal components we seek. The last is K, which is the number of iterations to preform. In all our experiments we 100 while for the rest we used
used c = 10-5 and = For the N2V experiments we used K K= 50.
Figure S1 depicts the convergence of the subspace iteration method for two different domains. For each noisy image and patch for which find the principal components (marked in red), the plot to the right shows the convergence of the first N = principal components. Specifically, for each principal component U; calculate its inner product with the same principal component in the previous iteration. As the graph shows, K 50 iterations suffice for convergence.
H THE IMPACT OF THE VACOBIAN-VECTIE DOT-PRODUCT LINEAR APPROXIMATION
As described in Sec 4.1, Alg. calls for calculating the Jacobian-vector dot-product of the denoiser. While for neural denoisers this calculation can be done via automatic differentiation, we propose using linear approximation instead (See Eq. (15)). This can reduce the computational burden, while retaining high-accuraee in the computed eigenvectors. For example, in an experiment using SwinIR and a 50, the cosine similarity between the principal components computed with the approximation and those computed with automatic differentiation typically reaches around 0.97 at the 50th iteration. However, terms of computational burden, the differences can sometimes be dramatic. For example, with the SwinIR model, when calculating one eigenvector for a patch of size 80 92, the memory footprint using automatic differentiation reaches 12GB, while using the linear approximation method it only reaches 2GB. These differences will increase for running on larger images. visual comparison of the resulting principal component can be found in Fig. $2.
I THE IMPACT OF ESTIMATING a
Our theoretical analysis is developed for non-blind denoising, and accordingly, most of our exper- iments conform this setting. These include the experiments on faces (Fig. 2 and Fig. S7),
24

---

Published as a conference paper at ICLR 2024
Noisy image
Subspace iterations convergence
1.0
0.8
0.4
3
0.2
10
20
30
40
50
1/2
0.8
03
0.6
0.4
10
20
30
40
50
Figure S1: Convergence of the subspace iteration method. In each row noisy image is shown with a red patch marking the region for which the posterior principal components are calculated. To the right, plot for each of the first 3 principal components the inner product between the principal component in consecutive iterations. As the graph shows, K = 50 iterations suffice to guarantee convergence in those domains.
MNIST digits (Fig. 4), on natural images (top part of Fig. 5, S5 and S6), and the toy problem of Fig. Namely, all those experiments the noise level a was assumed known.
Nevertheless, we show empirically that our method can also work well in the blind setting. This the case in the real microscopy images (bottom part of Fig. 5). this experiment, we estimated using the naive formula 8² à Ell (1) - ylla. where ja (1) is the (blind) N2V denoiser. It is certainly possible to employ more advanced noise-level estimation methods order to obtain an more accurate estimate for Indeed, noise-level estimation, particularly for white Gaussian noise, has been heavily researched, and as of today state-of-the-arr methods reach very-high pre- cision (Chen et al., 2015; Khmaa al., 2018; Kokil & Pratap, 2021; Liu & Lin, 2012: Liu al., 2013). For example, when the real a equals 10, the error in estimating sigma is around 0.05 (see e.g., Chen et al. (2015)). However, we find that even with the naive method described above, we get quite accurate results. Particularly, the impact of small inaccuracies in a on our uncertainty estimation turn out to be very small. To illustrate this, we applied our method with SwinIR model that was trained for a = 50, on images with noise levels of a 47.5, 52.5. This accounts for 5% errors in a, that are significantly higher than typical 0.5% errors of good noise level estimation techniques. Despite the inaccuracies in the eigenvectors produced using our method are quite similar, as can be seen in Fig. $3.
J USE IN NON-GAUSSIAL SETTINGS
In Sec. 5 empirically show our method provides sensible results also real microscopy images (bottom part Fig. 5), where the noise model is not known. In this setting, the noise distribution in each pixel is likely close to Poisson-Gaussiar the noise level is unknown, and it is not clear if the noise is completely white. However, the theory developed in this paper holds only for the non- blind Gaussian denoising case. We therefore aim to provide here intuition to why method can still find meaningful results for blind Gaussian denoising.
Suppose that the observation model is
y + on,
(S48)
25

---

Published as a conference paper at ICLR 2024
Using Lincar Approximation
Using Backpropogation
Figure S2: The impact the linear approximation on the calculated principal component. The first principal component calculated with SwinIR and 50, using the linear approximation in Eq. (15), and using automatic differentiation (Backpropegation). Both methods achieve similar re- sults, with a cosine similarity of 0.96 over 50 iterations. However, the linear approximation methods uses drastically less memory.
Noisy Image
D
-3,7,01-
(1
+3 7:D:
= 52.5
0 50
0 = 47.5
Figure S3: The effect of small inaccuracies in uncertainty estimation. The first principal component calculated using SwinIR, for an assumed a 50, for three different actual noise levels in the image.
26

---

Published as a conference paper at ICLR 2024
where x is random vector taking values in R the noise n NN(0, Ia) is multivariate Gaussian vector that is statistically independent of x, and the noise level or is a random variable sampled from some distribution Pa- The noise level is unknown to the denoiser (all that is known is the distribution of noise levels pa).
In this case, the Jacobian of the MSE-optimal denoiser is given by
01)
aE
x
= y
Oy
dy
E [Exx|y y, olly = =y]
dy
E
E[xly = y o) y y- y =
dy
= E Cov(xly = y: σ) σ)
y
(S49)
6²
where we used the law of total expectation in the second line, and Theorem 2 in the last line. Namely, instead of Cov(xlu) which we had in the Gaussian setting, here the Jacobian reveals the mean of the o²
posterior covariance divided by o². where the mean is taken over all possible noise levels a. This matrix is linear combination of the posterior covariances corresponding to different noise levels, so that it captures some notion of spread about the posterior mean, similarly to the regular posterior covariance that arises in the non-blind setting. Thus, intuitively, we expect that the top eigenvectors of this matrix capture meaningful uncertainty directions, similarly to the non-blind setting.
K VALIDATION OF THE PREDICTED PRINCIPAL COMPONENTS
It is impossible directly measure the quality of our estimated posterior PCs, since denoising datasets contain only one clean image æ for each noisy image y. This single x is just one sample from the posterior Pxly and therefore cannot be used to extract a ground-truth posterior covariance matrix or ground-trute truth PCs to compare against. To validate our method beyond the controlled toy- experiment of Fig. 3, in which the ground-truth posterior distribution was known analytically and thus so were the PCs, here we provide the two following experiments.
First, we employ the use of a diffusion based posterior sampler for evvess-prrblemm. to generate many posterior samples for a noisy image y. The posterior principal components can then by ex- tracted by performing PCA on those samples. We note that this approach is impractical for real- world applications because of its very high computational cost, and is brought here only for eval- uating our method against some baseline. Indeed, when using a diffusion- each sample requires many neural function evaluations (NFEs) to generate, and many samples are needed for obtaining accurate PCs. This is while our method can faithfully extract each posterior PC with only 10 NFEs, shown in the convergence graphs in Fig. S1.
posterior sampler,
For each noisy image, we generated many posterior samples using DDNM (Wang et al., 2023) and used them calculate the PCs of the posterior. As can be seen in Fig. S4, the number of posterior samples increases, the PCs estimated using this baseline become cleaner and more similar to our PCs. However, even with 500 samples, the PCs of this baseline do not seem to have fully converged, and generating 500 posterior samples using DDNM requires 50,000 NFEs. Therefore, for extracting e.g., PCs. our method roughly 1000 faster than this naive approach.
We further supply quantitative results in Tab. 1, over 100 randomly selected images from CelebA- 19 Baranchuk al. (2022), a subset of CelebAMask- mean of the samples generated by the posterior samples should theoretically approximate the pos- terior mean, which is the MSE-optimal restoration. As we verify, this estimate is indeed very close to our denoiser's ' output, and they both achieve practically the same RMSE to the ground-truth im- ages. Second, we compare the PCs our method to those generated by the suggested baseline by measuring the norm of the error after projecting it onto these PCs. The larger this norm, the larger the portion of the error that these PCs account for. Mathematically. this measure defined as | (ax nu [1, where V is matrix containing the PCs as columns. We compute the ratio of this norm relative to the measured MSE, and find that the mean of this measure for both methods is also very close.
Lee et (2020). First, the empirical
27

---

Published as a conference paper at ICLR 2024
500 Samples
50 Samples 100 Samples 200 Samples 500 Samples
Ours
Rotated
V1
V2
V3
V4
Figure $4: Comparison to PCs computed by applying SVD on different numbers of posterior samples. The posterior samples were produced using DDNM Wang et al. (2023). In the second column to the right, rotate the produced PCs best match our estimated PCs, while constraining them to remain orthonormal (using Procrustes analysis). The fact that the rotated PCs are quite similar our PCs, shows that both sets of PCs span similar subspaces. However, as can be seen, our PCs are more disentangled within that subspace.
Table 1: Quantitative evaluation of the estimated PCs.
RMSF(as, /n (v)) 4
|VT (v 14 (v)) /13 (1 1
Baseline Ours
4.026 10
7.767 10-3
4.022 10-2
7.638 10-3
Finally, we verify the predicted eigenvalues by comparing the projected test error over the first PC, of (x - (y)), to the predicted 1st eigenvalue Å₁. The average the ratio between those two quantities should theoretically be 1. For the same 100 randomly sampled face images, we found that the average of this ratio is 1.03. We also verified the predicted eigenvalues by calculating the ratio for the natural images domain. For this, randomly selected 100 natural images from the CBSD (Martin et al., 2001), Kodak (Franzen, 1999) and McMaster (Zhang et al., 201 1) datasets, and applied our algorithm using SwinIR (Liang al., 2021). For each image calculated the PCs
28

---

Published as a conference paper at ICLR 2024
Table 2: Quantitative evaluation of the estimated marginal posterior distributions.
Face images NLL +
Natural Images NLL &
Moments 1 &
1.83
0.05
Moments 1 4
1.81
0.03
Noisy Image
M4 (v)
+3, d,v
D2
V23
-3,7 -2
-1
1/41
2//7
3
Figure S5: Additional examples natural images using SwinIR (Liang et al., 2021). In each row, one of the first three PCs corresponding to the noisy image is shown on the left. On the right, images along the PC are shown above the marginal posterior distribution estimated for this direction. The principal components reveal uncertainty delicate parts of the wall-painting such as the thin rays of the sun, or the existence of mullions in the windows.
100 100 sized patch, located randomly within the image. For these images, the ratio computed was 0.93.
In addition, to quantitatively verify that the marginal posterior distributions estimate along the PCs are accurate, we measure the negative log likelihood (NLL) of the ground-truth images projected onto those directions (lower is better). We compared this the NLL of Gaussian distribution defined by only the first two estimated moments. Tab. 2 provides the results for the same 100 randomly selected face images and natural images. In both cases, the NLL of estimation is lower.
L ADDITIONAL RESULTS
Figures S5 and S6 provide additional results on test images from the McMaster (Zhang et al., 2011) dataset and images from ImageNet (Deng et al., 2009). In the supplementary material attach video showing more examples on face images, demonstrating different semantic principal compo- nents.
29

---

Published as a conference paper at ICLR 2024
Noisy Image
V/
33
-3,7
k2(+)
++3 +3377
3/4 -2 - - 1
1/4
2
3/24
Figure S6: Additional examples on natural images using SwinIR (Liang al., 2021). In each row, the first three PCs corresponding to the noisy image are shown on the left, and one is marked in blue. On the right, images along the marked PC are shown above the marginal posterior distribution estimated for this direction. The principal components catch semantic directions such the pattern on the owl's feathers, the embroidery pattern, or the length of the Axolotl's gills.
L.I POLYNOMIAL FITTING EXAMPLES
As discussed briefly in Sec. 5, we experimented with fitting polynomial to the function f(a) v 14 (y + av), and using the derivatives of the polynomial at 0 instead of using numerical derivatives of (a) itself at a = Here, we provide the results of an experiment where we polynomial degree six over the range [-VX, v for the ith principal component. As can be seen in Fig. $7, the marginal distribution estimates are quite smooth. Presumably, these posterior estimates are smoother than the true posterior, as the low degree polynomial smooths the directional posterior mean function.
30

---

Published as a conference paper at ICLR 2024
Noisy Image
1/2
V3
3
(2
+3/2ro
3/ -2 -1 A
0
1/74
2//4
3
Figure S7: Additional examples on face images, using polynomial fit marginal distribution estimate. In each row, the first three PCs corresponding to the noisy image are shown on the left, and one marked in blue. On the right, images along the marked PC are shown above the marginal posterior distribution estimated for this direction. The principal components highlight meaningful uncertainty, such as eyes shape or the existence of wrinkles. Note as an example in the first row how the optimal-MSE restoration is the mean of the more probable mode, depicting no hair on the forehead, and the distribution'ss tail, yielding the less- probable transtueenl hair.
31


Reviews:
[{"review_id": "kK4KjtlUPY", "reviewer": "Reviewer_HsMN", "summary": "This paper proposes a new method to compute posterior central moments of minimum mean squared error (MMSE) estimates from observations contaminated by additive white Gaussian noise. Then, this theoretical development is used to obtain uncertainty estimates in the context of image denoising. The method is cleverly adapted to the case where the denoiser is a black box (e.g., a neural network) only requiring forward passes through it, rather than back-propagation to compute derivatives exactly.", "questions": "I have no questions.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The method proposed in the paper is a relevant contribution to the important problem of estimating uncertainty in high-dimensional estimation problems, in this case, in image denoising from additive white Gaussian noise. The paper is very well written, with very precise and clear notation, and is a pleasure to read. The graphical illustrations of the results are very clear and useful and the experiments are also well presented, although not as impactful, as the variations depicted are very subtle.", "weaknesses": "The paper does not have, in my opinion, any major weaknesses, although a few minor things could be improved. \n\nWhen the authors first mention denoising (second line of the introduction), they cite a couple of papers, the oldest of which is from 2017. This may give the wrong idea that image denoising started in 2017, when in fact it is arguably the oldest and longest-standing problem in image processing, going back at least to the 1960s. More classical references should be mentioned here, rather than just a couple of recent deep-learning-based methods. \n\nThere are two common meanings for \"score\": gradient of the log-likelihood w.r.t. the parameters (more common in statistics) or w.r.t. the observations (more common in machine learning). It would be nice to make clear that you're using the second one, to make sure that some statistician reading the paper doesn't get confused. \n\nEquation (1) is not a denoising problem; it is the observation model underlying a denoising problem. \n\nAccording to Efron (2011), Tweedie's formula was derived by Robbins in 1956, which predates Miyasawa (1961). \n\nNotice that the first equality in Equation (4) is basically equivalent to Equation (2.8) in the paper by Efron (2011).\n\nMinor typos: \"memory efficient\" -> \"memory-efficient\"; \"...which connects between the MSE-optimal denoiser and the score function of noisy signals\" -> \"...which connects the MSE-optimal denoiser with the score function of noisy signals\"; \"...the most well known ...\" -> \"... the best known ...\" or \"...the most well-known ...\".", "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": "", "extra_scores": {}}, {"review_id": "OGo8IwWaWV", "reviewer": "Reviewer_a3nf", "summary": "Motivated by the recent success of using deep denoisers as priors (e.g, PnP and score matching),  this paper conduct connections between posterior central moments and the derivatives of posterior mean in a g higher order and recursive manor within the context of AWGN denoising. Then such conduction is applied for denoising  uncertainty visualization on multiple image datasets.", "questions": "1), While the authors claim the ability to compute the principal components (PCs) of the posterior distribution for any specified image region, the implementation details remain unclear to me. Seems there are no clear connection between the PCs and image features. At least it is not controllable.\n\n2), Figure 4, and Figure 5, the authors claims that the PCs show the uncertainty along meaningful directions. However, interpreting these findings remains challenging. Given that images often follow complex distributions, it's not evident how the curves presented by the authors establish a direct connection with changes in ground-truth geometric features.\n\n3), The computational complexity should be reported in the revision.\n\n4), More importantly, absolute error to the ground-truth should be also presented in Fig. 4 and Fig.5 to show that the uncertainty indeed can better reflect the restoration error.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1)\tThe paper is well organized and easy to follow. The demonstration proof about posterior central moments and directional projection in a recursion manner is well presented and informative.\n\n2)\tThe intuitive of using higher order moments to approximate posterior distribution is overall interesting. \n\n3)\tThe linear approximation in Eq.15 seems to be effective to handle differentiation complexity for higher-order tensor.\n\n4)\tFinally, the paper has shown its potential for uncertainty visualization within the context of image Gaussian denoising.", "weaknesses": "1), This paper only considers AWGN removal, making its practical relevance to real-world denoising (e.g., Poisson, Laplacian, Speckle, etc) very limited. Moreover, it would be more interesting to show how this method can be applied to other imaging inverse problems like inpainting and super-resolution both empirically and theoretically.\n\n2), At the same time, since the analysis is built on Gaussian distribution,  it is difficult to evaluate its theoretical contributions since the recurrence relation for the central moments of the Normal Distribution is somehow standard.\n\n3), In this paper, no baseline methods are compared against with. No quantitative results about the uncertainty calibration such as expected calibration error etc. are reported.", "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": "", "extra_scores": {}}, {"review_id": "aVWzhzTRXO", "reviewer": "Reviewer_EK8x", "summary": "Authors show that given a pretrained score-based model, they are able to approximate the higher-order moments of the posterior distribution of the unknown noise-free image given noisy image.", "questions": "See weaknesses.", "limitations": "", "rating": 3, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 1, "strengths": "**Originality.** While authors provide theoretical results for higher-order moments, as they mention, the second order moment has already been used in the context of denoising uncertainty quantification.\n\n**Quality and clarity.** The paper was not too hard to follow. There are a few typos and some sentences are not clear.\n\n**Significance.** It is unclear what this approach brings to the table that was not possible with existing uncertainty quantification methods.", "weaknesses": "* $\\sigma$^2 not defined in Eq. 4.\n\n* This particular method for image denoising uncertainty quantification is not motivated. Why would one use this approach? Why not just write down the posterior distribution using a likelihood model and the pretrained score-based model as (log) prior and use MCMC? Why not use amortized Bayesian inference and just simply sample from the posterior distribution?\n\n* Related to the previous point, a lot of emphasize has been given to higher-order moments of the posterior. Why do we care about computing these? Can't we compute these quantities using samples from the posterior distribution?\n\n* How would one use this method when the noise level is unknown?", "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": "", "extra_scores": {}}, {"review_id": "WeaV4guEDi", "reviewer": "Reviewer_N9TK", "summary": "A method is proposed to quantify uncertainties when solving denoising problems with data-driven pre-trained denoisers.  The high-order central moments of the target posterior distribution are related to higher-order derivatives of its mean.  For high-dimensional problems the top eigenvectors of the posterior covariance are computed to capture the main modes of variation for uncertainty visualization.  Experiments are presented to demonstrate application of the method on a number of data-sets.", "questions": "Is there any way to create modified images that more clearly show uncertainties?  I suppose increasing the scaling of the eigenvalue could be considered, although this is not very well motivated.  It might also be interesing to consider adding multiple eigenvectors at once, to show the mutliple modes of variation at once, rather than just one mode at a time.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "For denoising problems, denoisers implicitly encode a posterior for the data on which they are trained, which is related to the denoiser through Tweedie's formula.  Tweedie's formula explicitly relates the first posterior moment to the score of the data.  Similar relations hold for higher order momements of the posterior and higher order scores.  These expressions are presented in Theorem 1 and 2, which are proved in the appendices.  For high dimensional problems (e.g. imaging), computing these terms directly is not computationally feasible.  Instead the top eigenvectors are computed by the subspace iteration method.  These eigenvectors are then used to visualize the main modes of variation in the posterior.  The method requires only a pre-training denoiser and forward evaluations of the denoiser, avoiding any need for additional training.  It is fast and memory efficient and so suitable for high-resolution images.", "weaknesses": "While the methodology introduced in is very nice, the uncertainty quantification visualizations presented are a little underwhelming. The MNIST example (Fig. 4) clearly highlights the merit of the approach, where different potential modes of a 9 or 4 are apparent.  However, for the experiments performed on natural images (e.g. Fig. 5), it is difficult to distinguish much variation in the images.  I had to really zoom into these images to see the changes that are highlighed by the annotated arrows.  However, I suppose the eigenvectors themselves may be a more useful visualization aid since it is more clear what structure of the image changes.  Nevertheless, the difficulting in distringuishing changes in the modified images may limit the practical application of the method proposed.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)