Estimating High Order Gradients of the Data Distribution by Denoising
Authors: Chenlin Meng
Yang Song
Wenzhe Li
Stefano Ermon
Venue: NeurIPS 2021 Poster (2021)
Keywords: score matching
Langevin dynamics
denoising
generative modeling

Abstract: The first order derivative of a data density can be estimated efficiently by denoising score matching, and has become an important component in many applications, such as image generation and audio synthesis. Higher order derivatives provide additional local information about the data distribution and enable new applications. Although they can be estimated via automatic differentiation of a learned density model, this can amplify estimation errors and is expensive in high dimensional settings. To overcome these limitations, we propose a method to directly estimate high order derivatives (scores) of a data density from samples. We first show that denoising score matching can be interpreted as a particular case of Tweedie’s formula. By leveraging Tweedie’s formula on higher order moments, we generalize denoising score matching to estimate higher order derivatives. We demonstrate empirically that models trained with the proposed method can approximate second order derivatives more efficiently and accurately than via automatic differentiation. We show that our models can be used to quantify uncertainty in denoising and to improve the mixing speed of Langevin dynamics via Ozaki discretization for sampling synthetic data and natural images.

Estimating High Order Gradients of the Data Distribution by Denoising
Chenlin Meng Stanford University chenlin@cs. www.nnesssssss....oeooeeeeeeeeee ord. edu
Yang Song Stanford University yangsongocs.. .ttpssssssssstsffseee Sranford.co edu
Wenzhe Li
Stefano Ermon
Tsinghua University 1wz210mails ts tsinghua. edu. cn
Stanford University ermonccs. stanford... edu
Abstract
The first order derivative of a data density be estimated efficiently by denoising score matching, and has become an important component in many applications, such as image generation and audio synthesis. Higher order derivatives provide additional local information about the data distribution and enable new applications. Although they can be estimated via automatic differentiation of learned density model, this can amplify estimation errors and is expensive in high dimensional settings. overcome these limitations, we propose a method directly estimate high order derivatives (scores) of data density from samples. We first show that denoising score matching can be interpreted as particular case of Tweedie's s formula. By leveraging Tweedie's 's formula on higher order moments, we generalize denoising score matching to estimate higher order derivatives. We demonstrate empirically that models trained with the proposed method approximate second order derivatives more efficiently and accurately than via automatic differentiation. We show that our models can be used to quantify uncertainty in denoising and to improve the mixing speed of Langevin dynamics via Ozaki discretization for sampling synthetic data and natural images.
Introduction
The first order derivative of the log data density function, also known score, has found many applications including image generation [23, 24, 6], image denoising [20, 19] and audio synthesis [9]. Denoising score matching (DSM) [29] provides an efficient way estimate the score of the data density from samples and has been widely used for training score- score-based generative models [23, 24] and denoising [20, 19]. High order derivatives of the data density, which we refer to high order scores, provide a more accurate local approximation of the data density (e.g., its curvature) and enable applications. For instance, high order scores can improve the mixing speed for certain sampling methods [2, 18, 12]. similar to how high order derivatives accelerate gradient descent in optimization [1]]. In denoising problems, given noisy datapoint, high order scores be used to compute high order moments of the underlying noise-free datapoint, thus providing a way to quantify the uncertainty in denoising.
Existing methods for score estimation [8, 29, 25, 32], such as denoising score matching [29], focus estimating the first order score (i.e., the Jacobian of the log density). In principle, high order scores can be estimated from a learned first order score model (or even density model) via automatic differentiation. However, this approach is computationally expensive for high dimensional data and score models parameterized by deep neural networks. For example, given a D dimensional
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

distribution, computing the (n + 1) order score value from an existing n-th order score model by automatic differentiation is on the order of D times more expensive than evaluating the latter [25]. Moreover, computing higher-order scores by automatic differentiation might suffer from large estimation error, since a small training loss for the first order score does not always lead to a small estimation error for high order scores.
To overcome these limitations, we propose a new approach which directly models and estimates high order scores of data density from samples. We draw inspiration from Tweedie's formula [4, 16], which connects the score function a denoising problem, and show that denoising score matching (DSM) with Gaussian noise perturbation can be derived from Tweedie's ' formula with the knowledge of least squares regression. We then provide generalized version of Tweedie's formula which allows us to further extend denoising score matching to estimate high order scores. In addition, we provide variance reduction techniques to improve the optimization of these newly introduced high order score estimation objectives. With our approach, we can directly parameterize high order scores and learn them efficiently, sidestepping expensive automatic differentiation.
While our theory and estimation method is applicable to scores of any order, we focus on the second order score (i.e., the Hessian of the log density) for empirical evaluation. Our experiments show that models learned with the proposed objective can approximate second order scores more accurately than applying automatic differentiation to lower order score models. Our approach is also more computationally efficient for high dimensional data, achieving up to 500 speedups for second order score estimation on MNIST. In denoising problems, there could be multiple clean datapoints consistent with noisy observation, and it is often desirable to measure the uncertainty of denoising results. As second order scores are closely related to the covaraince matrix of the noise-free data conditioned on the noisy observation, we show that our estimated second order scores can provide extra insights into the solution of denoising problems by capturing and quantifying the uncertainty of denoising. We further show that our model can be used to improve the mixing speed of Langevin dynamics for sampling synthetic data and natural images. Our empirical results on second order scores, special case of the general approach, demonstrate the potential and applications of our method for estimating high order scores.
Background
2.1 Scores of distribution
Definition 1. Given a probability density p(x) over RR,,, we define the k-th order score S (x(x). RD @TRP, where @* denotes k-folld tensor multiplications, to be a lensor with the (ix, i2, (x)- logr(x), where (in, (2, (1,.... D)*
th index given by [Sx(x)] (x)]iria...) i2... A E2.
As an example, when l = 1, the first order score is the gradient of logp(x) w.r.t. to x, defined as sı (x) + Vx log pxx). Intuitively, this is vector field of the steepest ascent directions for the log-density. Note that the definition of first order score matches the definition of (Stein) score [8]. When k 2, the second order score is the Hessian of log p(x) w.r.t. to x. It gives the curvature of density function, and with si (x) it can provide better local approximation to log p(x).
2.2 Denoising score matching
Given a data distribution Pdata (x) and a model distribution p(x: 0), the score functions of Pdata (x)
and p(x; 0) are defined as S (x) Vx log Piua (x) and S S 0) Vx log p(x: 0) respectively. Denoising score matching (DSM) [29] perturbs data sample x Pdatu (x) with pre-specified noise
distribution (x x) and then estimates the score of the perturbed data distribution go (x) = S4 (x
x) (x)dx dx which we denote (x) + Vx log qa (x)). DSM uses the following objective
Eppntt (x) (x 2
(xxx)
S (x: 0)
Vx log (x x)]
(1)
is shown that under certain regularity conditions, minimizing Eq. (1) is equivalent to minimizing the score matching [8] loss between $ (x: 0) and $ (x) [29] defined as
2 Epan (x) Ey (x S (5: 0) 5, [x]
(2)

---

When (x x) = Nxx the objective becomes
1 LDSM (0) Expm(x) Egg (8)x)
(x; 0)
1 02 (x x) 2
Optimizing Eq. (3) can, intuitively, be understood as predicting X= the added "noise" up to constant, given the noisy input x, and is thus related to denoising. Estimating the score of the noise perturbed distribution (~(X) instead of the original (clean) data distribution Pamaa Pxaa(x) allows DSM to approximate scores more efficiently than other methods [8, 25J. When is close to zero, qo (X) Pdata (x) the score of go (x) estimated by DSM will be close to that of Pdatal (x). When large, the estimated score for qo (X) plays crucial role denoising [20] and learning score-based generative models [23, 24j.
2.3 Tweedie's formula
Given a prior density Pdata (x), noise distribution (x(xx) [xxxxx o²1), and the noisy density qo (x) ( (x x)xx, Tweedie's s formula [16, 4] provides a close-form expression for the posterior expectation (the first moment) of x conditioned on x:
E[ X = (x),
(4)
where (x) Vx log fo (x). Equation implies that given a "noisy" observation X ~ 9a (x), one can compute the expectation of the "clean" datapoint x that may have produced x. a result, Equation 4 has become an important tool for denoising [19, 20J. We provide the proof in Appendix B.
A less widely known fact is that Tweedies' formula be generalized provide higher order moments of x given X, which we will leverage to derive the objective for learning higher order scores.
3 Estimating Higher Order Scores by Denoising
Below we demonstrate that DSM can be derived from Tweedie's formula [4, 16]. By leveraging the generalized Tweedie's formula on high order moments of the posterior, we extend DSM estimate higher order score functions.
3.1 DSM in the view of Tweedie's formula
The optimal solution to the least squares regression problem
min E E (x | (x:0) 0) x
(5)
well-known be the conditional expectation h 0*) E[x x]. we parameterize h 0) x + d²S (x: 0) where (x; 0) is first order score model with parameter 0, the least squares problem in Eq. (5) becomes equivalent the DSM objective:
min Exmn(x) Bg [ "o [ 0) x x| min 2 LDSM(O).
0
From Tweedie's 's formula, know the optimal 0 satisfies h(x: 0*) X ( 0*) Ex x] - (x), from which we can conclude that $ (*) (x). This proves that minimizing
the DSM objective Eq. (6) recovers the first order score.
There are other ways to derive DSM. For example, [15] provides a proof based on Bayesian least squares without relying Tweedie's 's formula. Stein's 's Unbiased Risk Estimator (SURE) [27] can also provide an alternative proof based on integration by parts. Compared to these methods, our derivation be easily extended learn high order scores, leveraging Tweedie's formula.
general version of
3.2 Second order denoising score matching
As a warm-up, we first consider the second order score, and later generalize to any desired order. Leveraging Tweedie's 's formula on x| x] and [E[x x]. obtain the following theorem.
3

---

Theorem 1. Given a D-dimensional distribution p(x) and (x) (x(x))(Xxxx)ds we have
E xx
x]
f(x,
$1)
(7)
E| xxx
xx
xx
x h(x, $1 S2),
(8)
where f(x, Si. 52) and h(i, Si 52) are polynomials of S (x). (x defined as
x + o o [
(x)x
S2 (x) T (x) +
(9)
h $1 S2)
x o ك'S ( o
(10)
Here Si (x) and $2( denote the first and second order scores of go (x).
In Theorem 1, Eq. (9) is directly given by Tweedic's formula on E[xx" x]. and Eq. (10) is derived from Tweedie's formula on both [[[x x] and E[xx" x]. Given a noisy sample x, Theorem relates the second order moment of the first order score $ (x) and second order score $2(X) of qo (X). A detailed proof of Theorem is given in Appendix B.
In the same way as how we derive DSM from Tweedie's formula in Section 3.1, we can obtain higher order score matching objectives with Eq. (9) and Eq. (10) as least squares problem.
Theorem 2. Suppose the first order score (xx) is given, we can learn a second order score model (x: 0) by optimizing the following objectives
0* arg min E
Pian (x) (x(x))
(8(x)
xx (x)
So
(11)
2
0* arg min E
piam (x) Eg (K(b))
Xxx X²i -xx- h(x (X). (2 0))
(12)
where f(-) and h(.) are polynomials defined in Eq. (9) and Eq. (10). Assuming the model has an infinite capacity, then the optimal parameter 0* satisfies (3; 0*) S2( for almost any x.
Here Eq. (11) and Eq. (12) correspond the least squares objective of Eq. (7) and Eq. (8) respectively, and have the same set of solutions assuming sufficient model capacity. In practice, we find that Eq. (12) has much simpler form than Eq. (11), and will therefore use Eq. (12) in experiments.
3.3 High order denoising score matching
Below generalize approach to even higher order scores by (i) leveraging Tweedie's 's formula to connect higher order moments of x conditioned on to higher order scores of 9o (x); ; and (ii) finding the corresponding least squares objective.
Theorem 3. Ex@"x[x]
S.), where xxx E RDD denotes n-fold tensor multiplica- tions, f, (X, Si...., S,)) is polynomial of (x, (8),, (X)) and 5x(x) represents the k-th order score of go ( = Patso (x)4o dxxx)xx.
fx, $i....,
Theorem 3 shows that there exists an equality between (high order) moments of the posterior distribution of x given x and (high order) scores with respect to x. To get some intuition, for 2 the polynomial f. (xx $1. $2) is simply the function in Eq. (9). In Appendix B, we provide recursive formula for obtaining the coefficients of f, in closed form.
Leveraging Theorem and the least squares estimation of [[@"x[x], we construct objectives for approximating the k-th order scores S (x) as in the following theorem.
Theorem 4. Given score functions s (x) Sk- (x). a k-th order score model S (ix; 0). and
0* arg min Exp (x) Ba (& x) ll @ fx (X), S. (xx 11]
We have Sa 0+) = š (x) for almost all X.
As previously discussed, when a approaches such that (x) 2 Pdlaaa (x), S (x; 0*) well- approximates the k-th order score of Pdata(X) (x).

---

Learning Second Order Score Models
Although our theory can be applied scores of any order, focus second order scores for empirical analysis. In this section, we discuss the parameterization and empirical performance of the learned second order score models.
4.1 Instantiating objectives for second order score models
In practice, we find that Eq. (12) has a much simpler expression than Eq. (11). Therefore, we propose to parameterize Sy(X) with model (x; 0), and optimize $2(x: 0) with Eq. (12), which can be simplified to the following after combining Eq. (10) and Eq. (12):
1 2z
D (0) E (x) (x(x)
(x; 0) $ (x;
5 0) (X)
(13)
where z - X-x Note that Eq. (13) requires knowing the first order score $1 (x) in order to train the second order score model $2 xx; 0). We therefore use the following hybrid objective to simultaneously
train both $ (x; 0) and 2 (x; 0):
Gomn (0) (D)M + y (
(14)
where LDSM ()) defined Eq. (3) and 7 E R=0 a tunable coefficient. The expectation for D.SM (0) and (DSM(() Eq. (14) can be estimated with samples, and we optimize the following unbiased estimator
Epoio (0)
I N
Σ
S2 (2;;;;
( 0)5; (x;;
O)TT
I 2,2Z
1
S1 (1
Z:
(15)
a²
where we define Z, 8,-x, and [ are samples from go 9x (x) Pauta (x)
which can be obtained by adding noise to samples from Paata (x). Similarly to DSM, when - 0, the optimal model 52 @*) that minimizes Eq. (15) will be close the second order score of Pauas (x) because do (x) 2 Paata (x). When is large, the learned (x; 0) can be applied to tasks such as uncertainty quantification for denoising, which will be discussed in Section 5.
For downstream tasks that require only the diagonal of S2, we can instead optimize simpler objective
Ljoint-diag (0) DDS-M-ding ( +
LDSM
where
(16)
1-20 2
DD..Midiag (0) E., (x) Byg (81x) diag $2 (x; 0)) $ (x:: 0)
s (x;
0)
(17)
Here diag(-) denotes the diagonal of matrix and c denotes element- multiplication. Optimizing Eq. (16) only requires parameterizing diag( S2 0)), which can significantly reduce the memory and computational cost for training and running the second order score model. Similar to Cjuin (0)
estimate the expectation in Eq. (17) with empirical means.
4.2 Parameterizing second order score models
In practice, the performance of learning second order scores is affected by model parameterization. As many real world data distributions (e.g., images) tend to lie on low dimensional manifolds [13, 3, 21]. we propose to parametrize S2 xx; 0) with low rank matrices defined as below
2 x 0) <x 0) + ( (); 0)
where al: 0) :R. - DxD is diagonal matrix, (): .RP R Dxr is matrix with shape D x r, and & D is a positive integer.
4.3 Antithetic sampling for variance reduction
As the standard deviation of the perturbed noise approximates zero, training score models with denoising methods could suffer from a high variance. Inspired by variance reduction method for DSM [30, 26), we propose variance reduction method for D2SM
zz!
CD.SMAR Bxxnon(x) E, B2-N(0.)) (xx,), (bxx )2
o² + (xb(x,) +(8 2x((x))
5

---

where X =x, = x- oz and d/b = $2 + Si Instead of using independent noise samples, we apply antithetic sampling and use two correlated (opposite) noise vectors centered at x. Similar to Eq. (14), we define Vjoin.VR VR + 7 LDMM VR, where LDSM.VR is proposed in [30].
We empirically study the role of variance reduction (VR) in training models with DSM and D2SM. We observe that VR is crucial for both DSM and DySM when a is approximately zero, but is optional when is large enough. To see this, we consider a 2-d Gaussian distribution (V, I) and train S1 (x; 0) and (x; 0) with DSM and DySM respectively. We plot the learning curves in Figs. la and Ib, and visualize the first dimension of the estimated scores for multiple noise scales in Figs. Ic and 1d. We observe that when o = 0.001, both DSM and DySM have trouble converging after a long period of training, while the VR counterparts converge quickly (see Fig. 1). When a gets larger, DSM and D2SM without VR can both converge quickly and provide reasonable score estimations (Figs. lc and ld). We provide extra details in Appendix C.
0,0M-R 2010
0,0M 10.33
3
CGH COT)
5,5M vi.
USM (0:11
15.5 13
CoNS (0:3)
112M0 19.5)
Ground truth
Ground trust
Opch
Eooc
(a) D2SM loss
(b) DySM-VR loss
(c) Estimated S
(d) Estimated S2
Figure I: From left to right: (a) D2SM loss without variance reduction (a = 10 (b) D2SM loss
with variance reduction (a 10-3), )) (c) Estimated $1 (d) Estimated $2, where the estimation for D2SM (0.001) is too far from the ground truth to appear on the plot.
4.4 The accuracy and efficiency of learning second order scores
We show that the proposed method can estimate second order scores more efficiently and accurately than those obtained by automatic differentiation of a first order score model trained with DSM. We observe in our experiments that $ (x; 0) jointly optimized via Lysimt or Ljoint-diag has comparable empirical performance as trained directly by DSM, so we optimize $1 (X: 0) and $2 0) jointly in later experiments. We provide additional experimental details in Appendix C.
Learning accuracy We consider three synthetic datasets whose ground truth scores are available 100-dimensional correlated multivariate normal distribution and two high dimensional mixture of logistics distributions in Table 1. We study the performance of estimating $2 and the diagonal of $2. For the baseline, we estimate second order scores by taking automatic differentiation of s 0) trained jointly with (x 0) using Eq. (15) Eq. (17). As mentioned previously, $ (X: 0) trained with the joint method has the same empirical performance as trained directly with DSM. For our method, we directly evaluate $2 (x: 0). We compute the mean squared error between estimated second order scores and the ground truth score of the clean data since we use small a and qo (x) R Pdata (x) (see Table 1). We observe that S xc; 0) achieves better performance than the gradients of S 0).
Table Mean squared error between the estimated second order scores and the ground truth on 10 test samples. Each setup is trained with three random seeds and multiple noise scales a.
Methods
d 0.01
0.05
0.1
Methods
0.01
0.05
0.1
Multivariate normal (100-d)
Mixture of logistics diagonal estimation (50-d, 20 mixtures)
6: grad (DSM)
43.80=0.012
43.7620.001
43.05:10.001
$, grad (DSM-VR) (Ours)
20.51:10.55
26.13
U.53 25.391 0.50
8, grad (DSM-VR VR)
9.4010.049
9.3910.015
9.21200.20
18.43 0.11
18.50
0.25
17.8810.
(Ours, 15)
7.12+ 0.319
6.91+00078
7.03==0039
Mixture of logistics diagonal estimation (80-d, 20 mixtures)
(Ours, 20)
5.2411.065
5.07.10.047
5.13::0.065
$ grad (DSM-VR)
32.80 0.34
32.44 0.30
31.511 0.43
(Ours, 30)
1.06:0.03
2.05:
0.06:0.0M5
S: (Ours) 21.68 0.18
22,23 0.08
22.18 0.08
Computational efficiency Computing the gradients of $) 0) via automatic differentiation can be expensive for high dimensional data and deep neural networks. To see this, we consider two models a 3-yeer MLP and U-Net [17], which is used for image experiments in the subsequent sections. We consider a 100-d data distribution for the MLP model and a 784-d data distribution for the U-Net. We parameterize $, and 82 with the same model architecture and use batch size of 10 for both settings. We report the wall-clock time averaged in 7 runs used for estimating second order scores during test time on a TITAN Xp GPU Table We observe that $2(X: 0) is 500x faster than using automatic differentiation for $ (x; 0) on the MNIST dataset.

---

Data
Noisy
Only
With
Data
Noisy
Only
With
Figure 2: Denoising 2-d synthetic data.
The incorporation of S improves uncertainty quantification.
R
G
B
R
7
7
2
Clean
= 0.3
Clean
0.3
0=
0.8
Figure 3: Visualizations of the estimated covariance matrix diagonals on MNIST and CIFAR-I0. For CIFAR-10 images, we visualize the diagonal for R, G, B channels separately. Images corrupted with more noise tend to have larger covariance values, indicating larger uncertainty in denoising. Pixels in background have smaller values than pixels near edges, indicating more confident denoising.
Uncertainty Quantification with Second Order Score Models
Our second order score model $2 (x; 0) can capture and quantify the uncertainty of denoising on synthetic and real world image datasets, based on the following result by combining Eqs. (4) and (9)
Cov(x x
E[xx"
X
E[x
E
x o S² r²1
(18)
By estimating Cov[x x via (y()), we gain insights into how pixels are correlated with each other under denoising settings, and which part the pixels has large uncertainty. To examine the uncertainty given by our (2(x;; O), we perform the following experiments (details in Appendix D).
Synthetic experiments We first consider 2-d synthetic datasets shown in Fig. 2, where train ( 0) and $2 0) jointly with Lysint- Given the trained score models, we estimate E[x x and Cov(x x] using Eq. (4) and Eq. (18). We approximate the posterior distribution ((x[x) x) with conditional normal distribution N(x [E[x X], Cov[x x]). We compare our result with that of Eq. (4), which only utilizes $ (see Fig. 2). We observe that unlike Eq, (4), which is point estimator, the incorporation of covariance matrices (estimated by ( 0)) captures uncertainty in denoising.
Covariance diagonal visualizations We visualize the diagonal of the estimated Cov[x x for MNIST and CIFAR-10 [10] in Fig. 3. We find that the diagonal values are in general larger for pixels the edges where there multiple possibilities corresponding
the same noisy pixel. The diagonal values are smaller for the background pixels where there is less uncertainty. We also observe that covariance matrices corresponding smaller noise scales tend have smaller values the diagonals, implying that the more noise an image has, the more uncertain the denoised results are.
Full convariance visualizations We visualize the eigenvectors (sorted by cigenvalues) of Cov|x *] estimated by (x 0) Fig. 4. We observe that they can correspond to different digit identities, indicating uncertainty in the identity of the denoised image. This suggests Cov[x X can capture additional information for uncertainty beyond its diagonal.
4
i
4
The first
30
200
Figure 4: Eigenvectors of the estimated covariance matrix on MNIST. The first column shows the noisy images (a 0.5) and the second column shows clean images. The remaining columns show the first 19, plus the 30, 80 and 200-th eigenvectors of the matrix. We can see digit 7 and in the cigenvectors corresponding to the noisy 7, and digit 4 and 9 in the second row, which implics that the estimated covariance matrix can capture different possibilities of the denoising results.
Sampling with Second Order Score Models
Here we show that our second order score model $2(x; 0) can used to improve the mixing speed of Langevin dynamics sampling.
7

---

Table 2: Speed analysis: direct modeling vs. autodiff.
Oata
Lengieh
tergeve
Melscn
Pipeesion
D 10c (MB.P)
D 754 (U-Nec
Astociff
Dataset
Samples (5000)
Single chein
Samples (30000)
Single chain
7.
Langevo
Oup (ra)k2200)
546
1:91
688
182
Oun (000
1840 67.3
09.4 345
Table 3: ESS on synthetic datasets. Datasets are shown Fig. 5. We use 32 chains, each with length 10000 and 1000 burn-in steps.
Dataset
Samples (3000)
Single chain
Samples (9000) Single chain
Figure Sampling with Ozaki and Langevin dynamics. We tune the optimal step size separately for both methods. The number in the parenthesis (Column 2 and 4) stands for the iterations used for sampling. We observe that Ozaki obtains more reasonable samples than Langevin dynamics using 1/6 or 1/3 iterations. Column 3 and 5 show samples within a single chain with length 31000 and 1000 burn-in steps.
Dataset
Dannsct
Methoo
Motrie
Langevin
21.81
26.33
Oraki
28,89
46.57
6.1 Background on the sampling methods
Langevin dynamics Langevin dynamics [1, 31] samples from (xat(x) using the first order score function S₁ (x). Given a prior distribution T(x),, a fixed step size >> 0 and an initial value Xy T(x).
Langevin dynamics update the samples iteratively as follows
& = x 2 (x:-1) + Vezi,
(19)
where 2, NN(0, T). As 0 and oc, X, is sample from (x)a under suitable conditions.
Ozaki sampling Langevin dynamics with Ozaki discretization tion in (x) to pre-condition Langevin dynamics:
[28] leverages second order informa-
E,/2 x; x + M1 (x + 2 ze 24 2 N)
(20)
where M,-1 - (c (8 () and Σ,-1 (c2:51(81-1) (x It is shown
that under certain conditions, this variation improve the convergence rate of Langevin dynamics [2]. In general, Eq. (20) is expensive to compute due to inversion, exponentiation and taking square root of matrices, so we simplify Eq. (20) by approximating S2 (x:-1) 1) with its diagonal practice.
In our experiments, we only consider Ozaki sampling with S₂ replaced by its diagonal in Eq. (20). small σ, S1 S₁ and
S2- We observe that diag S. (x: 0)) in Ozaki sampling can be computed in parallel with $ (5:: θ) on modern GPUs, making the wall-clock time per iteration of Ozaki sampling comparable that of Langevin dynamics. Since only use the diagonal of $2 xc 0) in sampling, we can directly learn the diagonal of $2(X) with Eq. (16).
As
6.2 Synthetic datasets
We first consider 2-d synthetic datasets in Fig. 5 to compare the mixing speed of Ozaki sampling with Langevin dynamics. We search the optimal step size for each method and observe that Ozaki sampling can use a larger step size and converge faster than Langevin dynamics (see Fig. 5). We use the optimal step size for both methods and report the smallest effective sample size (ESS) of all the dimensions [22, 5] in Table 3. We observe that Ozaki sampling has better ESS values than Langevin dynamics, implying faster mixing speed. Even when using the same step size, Ozaki sampling still converges faster than Langevin dynamics the two-model Gaussian dataset consider (see Fig. 6). In all the experiments, we use o = 0.1 and we provide more experimental details in Appendix E.
6.3 Image datasets
Ozaki discretization with learned $ 0) produces more diverse samples and improve the mixing speed of Langevin dynamics on image datasets (see Fig. 7) To see this, we select ten different digits from MNIST test set and initialize 1000 different sampling chains for each image. We update the
8

---

Data
Orei samales
Ocal camples
Onki parplees
Ozei serilees
Otali ceriple
(a) Data
(b) 50 iterations
(c) 100 iterations (d) 200 iterations (e) 300 iterations
(f) 400 iterations
Iniiiization
Langeeth samples
Langovin sambles
1
Langastn samples
Langevin samples
Langevin sampless
(g) Initialization
(h) 50 iterations (i) 100 iterations
(i) 200 iterations (k) 300 iterations ()) 400 iterations
Figure 6: Sampling two mode distribution. We the same step size 0.01 for both methods.
We observe that Ozaki sampling converges faster than Langevin sampling.
9
7
4
4
5 3
5
5 7
4
7
5
3
Langevin
4
0
7
Q 7
の 5 6 2
7
8 S 5
2 9
4
9 4
4 8 7
u 2
7
8
3 2 4
5 7 5 5
3
4
4 6
Ozaki
7
D
3 7 5
2
4
6
2 iterations
6
9
7
itcrations
iterations
iterations
200 itcrations
Figure 7: Sampling on MNIST. We observe that Ozaki sampling converges faster than Langevin
dynamics. We use step size 0.02 and initialize the chain with Gaussian noise for both methods.
4
4
4
9
9
9
4
4 7 7 9 4
4 4
9
4 4
4
Ocaid
Languen (3.03)
2
3
2
2 3 Z /
2
2 2
2
7
2
6
2
Langevn 02)
Ocaid (3.01)
Onii
Langeve n (3.
4
4 2
4
4
4
4
4
4 4 4 4 4 4
4 4
9
8
4
7 9
4
2
2 2 2
6 2 2 6 2
2
2
2
2
55a (a) Percentage of changes in class label w.r.t. iterations.
Langvin dynamics
(b) Different chains initialized with the same left panel image after 1000 iterations of update with step size 0.03.
Figure 8: Sample diversity analysis. The number in the parenthesis in Fig. 8a denotes the step size. We initialize the chain with MNIST test images and report the percentage of images that have changed class labels from the initialization w.r.t. sampling iterations. We observe that Ozaki sampling has more diverse samples.
chains with Ozaki sampling and report the percentage of images that have class label changes after fixed number of sampling iterations in Fig. 8a. We compare the results with Langevin dynamics with the same setting and observe that Ozaki sampling has more diverse samples within the same chain in a fixed amount of iterations. We provide more details in Appendix E.
7 Conclusion
We propose a method directly estimate high order scores of data density from samples. We first study the connection between Tweedie's formula and denoising score matching (DSM) through the lens of least squares regression. We then leverage Tweedie's formula on higher order moments, which allows us to generalize denoising score matching to estimate scores of any desired order. We demonstrate empirically that models trained with the proposed method can approximate second order scores more efficiently and accurately than applying automatic differentiation to a learned first order score model. addition, we show that our models can be used to quantify uncertainty in denoising and to improve the mixing speed of Langevin dynamics via Ozaki discretization for sampling synthetic data and natural images. Besides the applications studied in this paper, it would be interesting to study the application of high order scores for out of distribution detection. Due to limited

---

computational resources, we only consider low resolution image datasets in this work. However, as direct next step, we can apply our method higher-resolution image datasets and explore its application to improve the sampling speed of score-based models [23, 24, 6] with Ozaki sampling. In low rank matrix, our training general, when approximating the high-order scores with a diagonal cost is comparable to standard denoising score matching, which is scalable to higher dimensional data. larger rank typically requires more computation but could give better approximations to second-order scores. While we focused images, this approach is likely applicable to other data modalities such as speech.
Acknowledgements
The authors would like to thank Jiaming Song and Lantao Yu for constructive feedback. This research was supported by NSF (#1651565, #1522054, #1733686), ONR (N000141912149), AFOSR (TA95501910024), ARO (N91NNT-21-1-0125) and Sloan Fellowship.
References
[1] Giovanni Bussi and Michele Parrinello. Accurate sampling using langevin dynamics. Physical Review E, 75(5)0056707, 2007.
[2] Arnak Dalalyan and Avetik Karagulyan. User-friendly guarantees for the langevin monte carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):5278-53)) 5311, 2019.
[3] Sanjoy Dasgupta and Yoav Freund. Random projection trees and low dimensional manifolds. In Proceedings of the fortieth annual ACM symposium Theory of computing, pages 537-546, 2008.
[4] Bradley Efron. Tweedie's formula and selection bias. Journal of the American Statistical Association, 106(496):1602-1614, 2011.
[5] Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), (3(2):123-214, 2011.
[6] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. arXiv preprint arXiv:2006. 1239, 2020.
[7] Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in 1076, 1989.
and Computation, 1853):0059-
[8] Aapo Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(Apr):695-70) 2005.
[9] Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro. Diffwave: A versatile diffusion model for audio synthesis. arXiv preprint vrXiv:2009.09761, 2020.
[10] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[11] James Martens and Roger Grosse. Optimizing neural networks with cooneckr-ffactore approx- imate curvature. In International conference on machine learning, pages 2408-2417. PMLR, 2015.
[12] Wenlong Mou. Yi-An Ma, Martin J Wainwright, Peter L Bartlett, and Michael I Jor- dan. High-order langevin diffusion yields an accelerated meme algorithm. arXiv preprint arXiv: /908. 10859, 2019.
[13] Hariharan Narayanan and Sanjoy Mitter. Sample complexity of testing the manifold hypothesis. In Proceedings of the rdd International Conference on Neural Information Processing Systems- Volume 2, pages 1786- 1794, 2010.
10

---

[14] Tianyu Pang, Kun Xu, Chongxuan Li, Yang Song, Stefano Ermon, and Jun Zhu. Efficient learn- ing of generative models via finieediifference score matching, arXiv preprint arXiv:2007.03317 2020.
[15] Martin Raphan and Eero P Simoncelli. Least squares estimation without priors or supervision. Neural computation, 23(2):374 420, 2011.
[16] Herbert Robbins. An empirical Bayes approach to statistics. University of California Press, 2020.
[17] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-assisted intervention, pages 234-241. Springer, 2015.
[18] Sotirios Sabanis, Ying Zhang, et al. Higher order langevin monte carlo algorithm. Electronic Journal of Statistics, 13(2):3805-388 -3500, 2019.
[19] Saeed Saremi and Aapo Hyvarinen. Neural empirical bayes. Journal of Machine Learning Research, 20:1-23, 2019.
[20] Saeed Saremi, Arash Mehrjou, Bernhard Schölkopf, and Aapo Hyvärinen. Deep energy estimator networks. arXiv preprini arXiv: 08306, 2018.
[21] Lawrence K Saul and Sam Roweis. Think globally, fit locally: unsupervised learning of low dimensional manifolds. Departmental Papers (CIS), page 12, 2003.
[22] Jiaming Song, Shengjia Zhao, and Stefano Ermon. A-nice-me: Adversarial training for meme. arXiv preprint arXiv: 1 7006007616 2017.
[23] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Advances in Neural Information Processing Systems, pages 11918-11930, 2019.
[24] Yang Song and Stefano Ermon. Improved techniques for training score- score-based generative models. arXiv preprint arXiv:2006.0901/ 2020.
[25] Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: A scalable approach density and score estimation. arXiv preprint 178000,905.07088 2019.
[26] Yang Song and Diederik Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.03288 2021.
[27] Charles M Stein. Estimation of the mean of a multivariate normal distribution. The annals of Statistics, pages 1135-1151, 1981.
[28] Stramer and RL Tweedie. Langevin-type models Diffusions with given stationary distribu- tions and their discretizations. Methodology and Computing in Applied Probability, 1(3)2883- 306, 1999.
[29] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computation, 23(7):1661 1674, 2011.
[30] Ziyu Wang, Shuyu Cheng, Li Yueru, Jun Zhu, and Bo Zhang. A wasserstein minimum velocity approach to learning unnormalized models. In International Conference on Artificial Intelligence and Statistics, pages 3728-3738, PMLR, 2020.
[31] Max Welling and Yee W Tch. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference machine learning (ICML-11), pages 681-688. Citeseer, 2011.
[32] Yuhao Zhou, Jiaxin Shi, and Jun Zhu. Nonparametric score estimators. In International Conference Machine Learning, pages 11513-11522. PMLR, 2020.


Reviews:
[{"review_id": "ofwDYr-SJO", "reviewer": "Reviewer_HV7H", "summary": "Using Tweedie's formula, the paper proposes a $n$-th order denoising score matching: a parametric model learns $n$-th order gradient of the log density of perturbed data. More specifically, the paper considers a data $x \\sim p(x)$ and a perturbed random variable $\\tilde{x} \\sim p(\\tilde{x})$, which is a marginal distribution of $\\tilde{x}$ for a given a conditional distribution $p(\\tilde{x}|x)$ and $p(x)$.\n\nTweedie's formula states that when $p(\\tilde{x}|x)$ is an exponential family distribution, the expected value of the $n$-th order moment of the posterior distribution of $x$ given $\\tilde{x}$ can be represented as a polynomial function of $\\tilde{x}$ and from-first-to-$n$th order gradients of the log marginal $\\log p(\\tilde{x})$. Nothing that the optimal solution of the least-squares regression to $x$'s moments for a given $\\tilde{x}$ is the expected value of the posterior moments, the paper proposes to learn $n$-th order score models by performing the regression. The authors emphasize that for guaranteeing optimality of $n$-th order denoising score matching, ground-truth scores should be accessible up to $n$-1-th order.\n\nNevertheless, the paper proposes to jointly train the first and second-order score models for empirical studies. In particular, the authors suggest a memory-efficient parameterization of the second-order score and a variance reduction method for the training. For the application of the joint training, the authors first demonstrate that the learned second-order scores provide information about the uncertainty of denoising. Second, the paper proposes to plug in the learned second-order gradients into the Ozaki discretization of Langevin dynamics, which uses the second-order gradients to improve the mixing speed of the Langevin dynamics. The experiments on synthetic data and natural images demonstrate that the second-order score model improves sampling quality for the same discretization steps compared to the Euler-Maruyama discretization (the first order).", "questions": "", "limitations": "", "rating": 8, "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": "(Originality & Significance)  \nIn my understanding, the contributions of the paper are clear:  \n(1) It introduces new methods to learn $n$-th order gradients of log density (of perturbed data).  \n(2) The authors discuss the conditions for the theoretical optimality of the proposed method and prompt that practical implementations may not guarantee the optimal.  \n(3) The paper proposes two new training techniques for improving practices with clear motivations, followed by their ablation studies.  \n(4) Finally, the paper demonstrates the practicality of learned higher-order scores through various experiments.\n\nI also consider that the results are essential. For example, while the ML community has grown its interest in score-based generative models, many SOTA implementations rely on a large number of discretization steps. As already stated in the paper, the proposed method can improve the practicalities of such models.\n\n\n(Quality & Clarity)  \nIn general, the paper has a well-organized structure to motivate the proposed $n$-th order denoising score matchings and other practical techniques to improve training. However, I found that the description of Tweedie's formula in Section 2.3 can be improved, considering searching the standard materials for the formula seems non-trivial. In particular, it would be nice if a revised version in the section can help readers understanding how the formula can be leveraged to learn $\\mathbb{E}[xx^T | \\tilde{x}]$ and $\\mathbb{E}[x | \\tilde{x}]$ in Section 3.2.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "(Limitations)  \nPlease find comments on the (Quality & Clarity) section in the main review.\n\n(Societal Impact)  \nN/A", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": ">12hrs"}, "extra_text": {}}, {"review_id": "nKWZZzC0yvG", "reviewer": "Reviewer_FpsM", "summary": "This paper utilizes Tweedie's formula to generalize denoising score matching to higher order derivatives. The proposed method can approximate second order derivatives more efficiently and accurately than differentiating first order derivatives. Experiments on synthetic data and small image datasets also demonstrates its applications to uncertainty estimation, and Langevin dynamics sampling.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper is well organized and easy to follow. Sec. 2, 3 and 4 progressively build up the idea of directly modeling and learning the second order score, and then develop a scalable parameterization and a variance reduction method for it. Experiments on synthetic data shows clear advantages of the proposed method over automatic differentiation of first order scores. Moreover, experiments on small image datasets gives promising results on uncertainty estimation, as well as on generative modeling.\n\nThe main weakness of the paper is that the proposed method is only evaluated on small datasets, i.e., MNIST, CIFAR-10 for uncertainty estimation, and MNIST for Langevin dynamics sampling. It is not clear how well the method can scale to higher-dimensional data, and how the diagonal/low-rank approximation compares to other generative models in such cases.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "As mentioned in Appendix F, the authors consider low resolution image datasets in this work due to limited computational resources.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "gs4PZsxCSRB", "reviewer": "Reviewer_5WYk", "summary": "The paper generalizes denoising score matching, which estimates the score, i.e. the gradient of the log density, of perturbed data. In particular, it shows how higher-order derivates of perturbed data densities can be estimated and learnt. Much of this is derived via interesting connections of denoising score matching to Tweedie's formula, a technique for denoising noisy data. While the paper theoretically shows how to learn derivatives of the log data density of arbitrary order, its experiments focus on second-order derivatives, corresponding to the Hessian of the log density. The learnt second-order score is leveraged for improved image (and toy data) denoising, where the second-order score allows to parametrize a covariance matrix, that captures the uncertainty in the denoising operation. Furthermore, it uses the second-order score for sampling applications where it is leveraged to design more efficient samplers via preconditioning. It is also shown that directly learning the second-order score as proposed is more accurate and faster than computing it via auto-differentiation from first-order scores at test time.", "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": "__Overall Impression and Significance:__ I enjoyed reading the paper and think that it addresses an important problem. Regular denoising score matching has recently found various applications, for example in generative modeling and image denoising. I can imagine that many of these applications will benefit from access to higher order scores, which can be estimated and learnt with the methods described in this paper. Therefore I consider the contributions as quite significant. However, some concerns about the method's scalability remain. All experiments are run either on small toy distributions or small images such as MNIST and the proposed method is potentially memory and compute intensive for high dimensional data.\n\n__Originality:__ To the best of my knowledge the proposed ideas, i.e. how to estimate and learn higher-order generalizations of the score, are original and novel. I found many of the derivations quite elegant and insightful. I also found the experiments insightful and well designed.\n\n__Clarity:__ The paper is well written and easy to follow.\n\n__Questions, Concerns:__\n1. My main concern is the method's scalability: Parametrizing second-order scores seems to be generally memory intensive as we need to form very high-dimensional tensors (unless we only do the diagonal version of course). Hence, does this imply practical limitations of the method? How easily could one estimate second-order scores of high resolution images, say 1024x1024 pixels, for example (this would be very useful in denoising applications, potentially)? How can we scale the method towards that?\n2. Also, how much slower is training of these models compared to regular first-order denoising score matching only?\n3. Do you envision any practical applications of scores of even higher order, like third-order scores?\n\n__Minor Feedback:__\n1. Eq. (13): I think the $z$'s should be properly introduced in the text. There is also no expectation over them. It becomes only clear below Eq. (15). \n2. In line 161, the authors refer to an expectation in Eq. (14). However, Eq. (14) has no expectations. I assume this is a typo and it should be Eq. (13).\n3. In Eq. (14), the authors use $\\gamma$, in Eq. (15) $\\lambda$ for the weight of the regular DSM objective term. I would suggest a consistent notation.\n\n__Conclusions:__ I think this is a very nice paper. The method essentially gives practitioners access to second-order scores in a similar manner like regular scores are learnt with standard denoising score matching. I think this will probably lead to further interesting applications, for example in generative modeling. Some concerns regarding the scalability of the method remain, considering the small experiments. Nevertheless, I recommend this paper for acceptance.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "__(Negative) Societal Impact:__ Has been discussed adequately in the Appendix.\n\n__Limitations:__ Have been discussed very briefly in Appendix and could potentially be discussed a bit more thoroughly. In general, I think we pay a price for all the advantages we get by having a second order score model. As was briefly mentioned, parametrizing a full second order score model can be memory intensive, due to the high dimensionality of this object, and I assume learning and convergence also becomes a bit slower than regular score matching. I could well imagine that this somewhat limits the applicability of the method on very high dimensional data, such as large high resolution images (denoising of large high resolution images would be very practically relevant). I think that these aspects could be discussed in more detail.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "17oMJjRlzv1", "reviewer": "Reviewer_yoPv", "summary": "- The author show that denoising score matching with Gaussian noise can be derived from Tweedie's formula through the lens of least square regressions.\n- This provides a new interpretation of the first score (gradient of log density).\n- They use the generalization of Tweedie's formula to higher orders to derive a similar estimation of higher order scores.\n- This is accompanied by a thorough experimental validation for the second order score.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The good:\n- The paper is well written, and easy to understand, and the task is well motivated. Being able to compute higher order scores more efficiently is very important, as those can be used to improve a lot of machine learning tasks like denoising and data generation.\n- The theoretical results are sound, and the connection of DSM to Tweedie's formula is interesting. \n- They include a thorough experimental analysis  They focus on the second order score in the experimental validation, and show that their method can learn the second order score more efficiently and accurately. Moreover, they apply the second order score to quantify the uncertainty of  denoising, and for Ozaki sampling.\n- Overall a solid paper in terms of presentation, but also contribution and content.\n\n\nWhat can be improved:\n- Discussion and limitations of the paper, as well as the societal impact were relegated to the Supplemental Material.\n- The tensor notation was used in section 2.1, but was only explained in section 3.3.\n- Similarly, the notation $\\mathbf{z}=\\frac{\\tilde{\\mathbf{x}} - \\mathbf{x}}{\\sigma}$ was introduced few lines after it ws first used in equation (13).\n- The notations can be improved: the authors use $\\mathbf{s}_1(\\tilde{\\mathbf{x}})$ for the score of the perturbed density. I think this is a bad notation, as the letter used to denote the argument of a function is mutable. It can be improved by using $\\tilde{\\mathbf{s}}_1$ for e.g. In fact, in Theorem 3 (and more generally in Theorems 1 and 2) they use $\\mathbf{s}_k$ as an argument of $\\mathbf{f}$, and it is not clear if they mean $\\mathbf{s}_1(\\tilde{\\mathbf{x}})$ or $\\mathbf{s}_1(\\mathbf{x})$; I also believe there is a mistake in line 138, where it should've been $\\mathbf{s}_1(\\tilde{\\mathbf{x}})$ and not $\\mathbf{s}_1(\\mathbf{x})$.\n- The word \"score\" is used without specifying the density. It can sometimes be confusing to know which score is referred to, for e,g, line 141.\n- Do equation (11) and (12) have the same set of solutions? This should probably be mentioned (unless I somehow missed it) to justify using (12) over (11) in practice.\nFigures 1c and 1d are a bit hard to read, as most of the lines are stacked on top of each other. may be a change of scale can help make it clearer?\n- In line 178, do you mean $D \\times r$ instead of $D \\times (D \\times r)$?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors included a discussion of the limitations and impact of the paper, but it was relegated to the Supplemental Material. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers are in unanimous agreement for acceptance. Please incorporate the reviewers' feedback for the camera-ready version. Congratulations on nice work!


Decision: Accept (Poster)