Learning normalized image densities via dual score matching
Authors: Florentin Guth
Zahra Kadkhodaie
Eero P Simoncelli
Venue: NeurIPS 2025 poster ()
Keywords: energy-based models
image dimensionality
score matching
diffusion models
manifold hypothesis
negative log-likelihood
Primary area: probabilistic_methods
TL;DR: We train an energy-based model on image datasets through a dual score matching objective and analyze the local geometry of the learned energy landscape.

Abstract: Learning probability models from data is at the heart of many machine learning endeavors, but is notoriously difficult due to the curse of dimensionality. We introduce a new framework for learning \emph{normalized} energy (log probability) models that is inspired by diffusion generative models, which rely on networks optimized to estimate the score. We modify a score network architecture to compute an energy while preserving its inductive biases. The gradient of this energy network with respect to its input image is the score of the learned density, which can be optimized using a denoising objective. Importantly, the gradient with respect to the noise level provides an additional score that can be optimized with a novel secondary objective, ensuring consistent and normalized energies across noise levels. We train an energy network with this \emph{dual} score matching objective on the ImageNet64 dataset, and obtain a cross-entropy (negative log likelihood) value comparable to the state of the art. We further validate our approach by showing that our energy model \emph{strongly generalizes}: log probabilities estimated with two networks trained on non-overlapping data subsets are nearly identical. Finally, we demonstrate that both image probability and dimensionality of local neighborhoods vary substantially depending on image content, in contrast with conventional assumptions such as concentration of measure or support on a low-dimensional manifold.

Learning normalized image densities via dual score matching
Florentin Guth
Zahra Kadkhodaie
Center for Data Science, New York University Flatiron Institute, Simons Foundation florentin, guth@nyu.edu florentin.guth@nyuuee
Flatiron Institute, Simons Foundation zk388@nyu. edu edu
Eero Simoncelli New York University
Flatiron Institute, Simons Foundation
eero. wwm.ncellmmmynn edu
Abstract
Learning probability models from data is at the heart of many machine learning endeavors, but notoriously difficult due the curse of dimensionality. We introduce new framework for learning normalized energy (log probability) models that is inspired by diffusion generative models, which rely on networks optimized to estimate the score. We modify score network architecture to compute an energy while preserving inductive biases. The gradient this energy network with respect to its input image is the score of the learned density, which can be optimized using denoising objective. Importantly, the gradient with respect to the noise level provides an additional score that can be optimized with novel secondary objective, ensuring consistent and normalized energies across noise levels. We train an energy network with this dual score matching objective on the ImageNet64 dataset, and obtain a cross-entropy (negative log likelihood) value comparable to the state of the art. We further validate our approach by showing that our energy model strongly generalizes: log probabilities estimated with two networks trained on non- overlapping data subsets nearly identical. Finally, demonstrate that both image probability and dimensionality of local neighborhoods vary substantially depending image content, in contrast with conventional assumptions such concentration of measure or support on a low-dimensional manifold.
Introduction
Many problems in image processing and computer vision rely, explicitly or implicitly, on prior probability models. However, learning such models by maximizing the likelihood of a set of training images difficult. The dimensionality of the space (i.e., the number of image pixels) large, and worst-case data requirements for estimation grow exponentially (the "curse of dimensionality" "). The machine learning community has developed a variety of methods to train a parametric network to estimate log p(x), known as an "energy model" (Hinton et al., 1986; LeCun et al., 2006), relying on the inductive biases of the network to alleviate the data requirements. For all but the simplest of models, this approach is frustrated by the intractability of estimating the normalization constant (Hinton, 2002; LeCun and Huang, 2005; Yedidia et al., 2005; Gutmann and Hyvärinen, 2010; Dinh et al., 2014; Rezende and Mohamed, 2015; Dinh et al., 2017; Song and Kingma, 2021).
clever means of escaping this conundrum is to estimate the gradient of the energy with respect to the image (known as the "score"), which eliminates the nommalization constant, and can be learned
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

from data with a "score- mattrr-matchin"" objective (Hyvärinen and Dayan, 2005). The recent development of "diffusion" generative models (Sohl-Dicksterd et al., 2015; Song and Ermon, 2019; Ho et al., 2020; Kadkhodaie and Simoncelli, 2020) builds on this concept, by estimating a family of score functions for images corrupted by Gaussian white noise at different amplitudes. These scores may then be used to draw samples from the corresponding estimated density, using an iterative reverse diffusion procedure. These methods have enabled dramatic improvements in both the quality and diversity of generated image samples, but the learned density is implicit. An explicit and normalized energy model (and density) can be obtained through integration of the divergence of the score vector field along trajectory (Song et al., 2021a), at tractable but substantial computational cost.
Here, we leverage the power of diffusion models to develop robust and efficient framework for directly learning a normalized energy model from image data. We approximate the energy with a deep network that takes as input both a noisy image and the corresponding noise variance, and derive two separate objectives. The first, obtained by differentiating with respect to the noisy image, is a denoising objective (as used in diffusion models). The second, obtained by differentiating with respect to the noise variance, ensures consistency of the energy estimates across noise levels, which we show to be critical for obtaining accurate and normalized energies. We optimize the sum of the two, in a procedure that we refer to as "dual score matching". We also propose novel architecture for energy. computed as the inner product of the input and output images of a score network. This preserves the inductive biases of the base score network, leading to equal superior denoising performance (and thus, sample quality).
We train our energy model on ImageNet64 (Russakovsky et al., 2015; Chrabaszez et al., 2017), 7), and show that the estimated energies lead to negative log likelihood (or cross-entropy) values comparable to the state of the art. We further demonstrate that the energy model strongly generalizes in the sense of Kadkhodaie et al. (2024): two separate models trained on non-overlapping subsets of the training data assign essentially the same probabilities to each image. This convergence is observed for training set sizes far smaller than the worst case prediction of the curse of dimensionality, We find that the distribution of log probabilities over ImageNet images covers a broad range, with densely textured images at the lower end, and sparse images at the upper end. The probability is relatively insensitive to changes in image luminance, but decreases with dynamic range. Finally, we highlight two geometrical properties of the learned image distribution. The first is an extremely tight inverse relationship between volume and density that leads to absence of concentration of energy values. They furthermore follow a Gumbel distribution, revealing an unexpected form of statistical regularity. The second is that the dimensionality of the energy landscape in the neighborhood of image varies greatly depending on image content and neighborhood size. We observe images with full-dimensional neighborhoods of non-negligible size and images with lower-dimensional neighborhoods even at sub-quantization scales. These results challenge traditional presuppositions regarding high-dimensional distributions, such as the concentration of measure phenomenon (Vershynin, 2018; Wainwright, 2019) and the manifold hypothesis (Tenenbaum et al., 2000; Bengio et al., 2013). Pre-trained models and software for running all experiments are available at https: //gtttu.....//////// com/ com//Florentii hww.aaaaa/aaaaaaoeaatttttttteeia
2 Learning normalized energy models with dual score matching
Estimating high-dimensional probability density from samples faces two significant challenges as a result of the curse of dimensionality. The first is statistical: realistically-sizer datasets do not contain enough information to learn an unknown high- high-dimensional distribution. One thus needs powerful inductive biases, typically in the form of a parametric network architecture, to hope to recover the data distribution. The second challenge is computational: the traditional objective aims to maximize the likelihood of the model over the data, which is intractable due to the need to compute the normalization constant. Nevertheless, recently developed generative models seem to succeed despite these problems. Drawing inspiration from diffusion models (Section 2.1), we develop a more direct solution by deriving a novel objective (Section 2.2) and architecture (Section 2.3) for learning normalized log probabilities (energies) from data. We validate our approach in Section 2.4.
2.1 Motivation
Traditionally, energy models are defined in terms of a parametric function (z(2) that approximates the unnormalized log density over RR:: pa (x) 2, ----uo -Uo(x) with normalizing constant

---

Energy U
Ground trith
Log produblity ((8)
Stugle SM
Single SkI (aljaate))
SM
-10-
Space
-15
Figure 1: Comparison of single and dual score matching recovering the energy of a scale mixture of two Gaussians in d 1000 dimensions. Experimental details are provided in Appendix C.3. Left: Radial slices of the log probability. The single score matching estimate (green dashed curve) fails to recover the true energy (blue solid curve), even after global normalization (green dotted curve), while dual score matching (red dashed curve) succeeds. Middle: Radial components of the scores. Single score matching learns an accurate score over the support of the data (blue bar plot) but not outside of it. Right: Energy landscape across space and time (noise level) for mixture of two Gaussians in one dimension. The direct path between the modes at = 0 crosses a large energy barrier (green curve), which is alleviated path that is not restricted to t 0 (red curve).
Zo = Jo (-Uo)))) The parameters 9 are estimated by minimizing the expected negative log likelihood (NLL), E., log po(r)] Er, [Uo[]]] + log Ze. which is equivalent to minimizing the KL divergence between (r) and the data distribution. The normalizing constant Zo plays a critical role in learning, representing the total energy over H² which trades off against the energy of the data. Unfortunately, direct estimation (i.e., computing the integral) is typically intractable.
The normalization constant be eliminated from the NLL by differentiating w.r.t. x, yielding quantity known as the (negative) score: V. log po = V (z)). As a result, a score model can be efficiently fitted to data via "score matching" (Hyvärinen and Dayan, 2005), which minimizes the Fisher (as opposed to KL) divergence between the model and the data. For data corrupted by Gaussian white noise, this amounts to solving denoising problem (Vincent, 201 Raphan and Simoncelli, 201 1). But this computational advantage comes at a statistical cost: good approximation of the score does not always lead to a good model of the energy (Kochler al., 2023), as we now illustrate. Consider the case of an equal mixture of two multivariate Gaussian distributions with zero mean and variances of and o² In high dimensions, this distribution concentrates near the two hyper-spheres with radii a v and 02 v d, leading to data scarcity in the rest of the space. We show in Figure energies estimated from samples by optimizing their gradient via (single) score matching. This fails to recover the true energy, even after adjusting a normalization constant. Indeed, evaluating the energy difference between the two Gaussians requires good estimation of the score along integration path between them. If there is no data in-between modes to constrain the score, due to an energy barrier concentration phenomena, this leads to inconsistent energy values across modes. In other words, single score matching estimates energy values up to a mode-dependent additive constant.
Diffusion models expand on score matching by learning scores of data corrupted by white Gaussian noise for a range of different noise amplitudes, -V Uy, t), with y + (o, tId). The evolution of the density p(yl() "time" (noise variance) l. increases is a diffusion process, and thus the corresponding energies increase in smoothness with This multiscale family of scores can be used to draw high-quality samples from p(z) using reverse diffusion algorithm, which follows a trajectory of partial denoising steps (Song and Ermon, 2019; Song et al., 202 Ho et al., 2020; Kadkhodaic and Simoncelli, 2020). In fact, the diffusion scores implicitly capture density model of the data (Song et al., 202 Kingma et al., 2021): the relative energy levels between modes, to which score matching at = 0 is blind (Zhang et al., 2022), are encoded in the score at the time when they merge (Raya and Ambrogioni, 2023; Biroli et al., 2024). This implicit density model can be evaluated in various ways. which all involve tractable but costly integration of score divergences (or denoising errors) as a function of time (Song et al., 202 la Kong et al., 2023; Skreta et al., 2024; Karczewski al., 2025). See Appendix A.2 for more in-depth exposition.
Intuitively, diffusion models deal with multimodal distributions by providing a high-probability path between modes through space and time, visualized in the right panel of Figure 1. In other words,
3

---

the joint distribution ply. t) qualitatively has a connected support. This suggests that we may learn the joint "space-time" energy ty, by score matching on (y. t). Specifically, the "space score" -V (y, 1) can be learned using denoising objective, and the "time score" -0,06(v 1) (Choi et al., 2022) can be learned using an analogous score matching objective. Matching both space and time scores correctly constrains the energy levels across modes: the energy difference between two modes can be explicitly recovered by integrating the energy derivative along the red path in Figure (although we will not need to do explicitly).
2.2 Dual score matching: Objective function
Space and time score matching. We aim to learn time-dependent energy function, Uo l), to approximate the NLL of the noisy image distribution at all noise levels l:
Uty t) = log (f p((z)ee de fall-~al w/l? & log(2m)) dz
(1)
Differentiating this NLL with respect gives the (negative) "space score", which can be expressed using the M Tweedde identity (Robbins, 1956; Miyasawa, 1961) as
y/
Vy,(y.t) = Ex
y
(2)
t
This leads the denoising score matching objective (Vincent, 2011; Raphan and Simoncelli, 2011; Saremi et al., 2018) used to train diffusion models:
fDSM (0, 1) = E.y.y VoUo(y, 1/
(3)
Differentiating the energy with respect to t yields (negative) "time score" (Choi al., 2022), for which we can derive a similar identity (see Appendix B.1):
d
11y
(,U 1)
F.
y
(4)
2t
2+
This leads an analogous "time score matching" objective:
êTSM (0, = Ex [
(acitte
11
x
+
(5)
26
21
Combining objectives across noise levels. The two objectives in eqs. (3) and (5) are defined for a given noise level f. We form an overall objective by integrating these over f, with an appropriate weighting. For the denoising score matching objective, natural choice is the so-called maximum- likelihood weighting (Song et al., 2021b; Kingma et al., 2021), which provides a bound on the KL divergence with the data distribution:
/ (0, (),
1
KL (1))
S
(6)
2
where 69 is the implicit density of the generative diffusion model. This integral can be approximated by Monte-Carlo sampling from a distribution of noise levels. In practice, have found it best to sample log uniformly over a finite interval (specifically, p(() x 1/1, E [kminy /max)), which corresponds to the following implementation of the integral:
DDSM
(0, (OSM(0,C)
dt [t (O, t)).
(7)
Note that the resulting term in the expected value is unit-less: it is invariant simultaneous rescaling of the data and noise level. We thus choose to weight the time score matching objective in eq. (5) by 1², so that it is also unit-less, and to evaluate it over the same distribution p()). Finally, after appropriate normalization by the dimensionality d that ensures that the two objectives have comparable orders of magnitude, we simply add them (we found no significant improvement from tuning tradeoff hyperparameter). In summary, our dual score matching objective is:
{(0)= = Eg (D.M(0,t) + (2) (ISM((,))
2
(8)

---

Normalization. Since we match only the partial derivatives of (v(v, t) to those of the true energy Uly t). we only recover the energy up global constant: the end of training, (y, t) s U(v, () + const. Note again that this crucially relies on p(y, 1) having a connected support. An important aspect of our framework is that enables estimation of this constant, which determines the normalization of -Ua(9.4), Indeed, the time score objective ensures that this normalizing constant does not depend on time: mass is conserved through the diffusion. Since the true distribution is approximately Gaussian at large noise levels, p(y|fmax 2 (0, Id)x we can set this constant to the entropy of this Gaussian distribution:
1o(y,1) - () (y, () [max] + d log((2relmax)-
Figure provides a high-dimenniooal numerical example verifying that dual score matching indeed provides an accurate estimate of normalized energy values (in particular, for t = 0), in contrast to single score matching.
Related approaches. Similar combinations of space and time scores were considered by Choi et al. (2022) (who referred to it as the "pathwise" method) and Kobler and Pock (2023). These time score objectives, however, relied on second derivative in time instead of regression objective. Our time score objective is also a special case of the conditional time score matching objective developed concurrently in Yu et al. (2025). Finally, Yadin et al. (2024) train an energy model with objective combining score matching with a classification cross-entropy loss for estimating a discretized version of the noise level f.
2.3 Dual score matching: Architecture
How should one choose an architecture to compute the energy Ug(y, t)? Rather than designing one from scratch (Salimans and Ho, 2021; Cohen et al., 2021; Thiry and Guth, 2024), construct one by modifying a score based denoising architecture that is known to have appropriate inductive biases. Let so: Rd IR - Rd be such architecture (e.g., a UNet). We wish to define new energy architecture Ua: IR x IR + R such that V VUU0 = se. preserving the inductive biases of sg. To achieve this, we set
Uolu, t) (y, so(u,)))
(10)
This inner product parameterization has been concurrently proposed by Thornton et al. (2025). We show in Appendix B.2 that (y. t) = (y. t) if the score network sg is both conservative and homogeneous (Romano al., 2017; Reehorst and Schniter, 2018). Homogeneity has been shown to hold approximately in the related setting of blind denoisers (Mohan* et al., 2020; Herbreteau al., 2023), and can be enforced architecturally. Note that the seemingly similar choice of a squared norm architecture Us (y, t)- 1 +|se@(u (y, t)l used in some previous energy models (Salimans and Ho, 2021; Hurault et al., 2021; Du et al., 2023; Yadin et al., 2024; Thiry and Guth, 2024) leads to the same desirable homogeneity properties y but not in sa, and thus fails to preserve the optimization properties of the original score network. Further architectural details are provided in Appendix C.I.
2.4 Performance evaluation
Denoising performance. We first verify that the gradient of the proposed energy network, VoU provides as good a denoiser as the score network 8g on which it is based. We train the two networks separately, Ug with the dual score matching objective (eq. (8) using double back- paek-propagation), and se with the standard (denoising) score matching objective (eq. (3)). Further details are provided in Appendix C.2. Table compares denoising performance across noise levels. For all but the smallest noise levels, the energy-based model achieves (slightly) better denoising performance than the score model. Thus, there no penalty in modeling the energy rather than the score. This is contrary to the results of Salimans and Ho (2021), and is likely due to our use of an architecture that is homogeneous, and for which non-conservativenes. of the score is not a large source of denoising error (Mohan* et al., 2020; Chao et al., 2023). Finally, these results demonstrate that the two components of our dual score matching objective are not trading off against each other; rather, they complement and even reinforce each other (i.c. their minima coincide).
5

---

Table Denoising performance several noise levels for corresponding score and energy networks, averaged over images in ImageNet64, with intensities scaled to lie in [0. 1]. All quantities are expressed as peak signal-to-noise ratio (PSNR) in dB: PSNR - - 1o
Noise variance
90
75
60
45
30
15
-15
-30
Score network
90.20
75.47
60.43
47.19
35.58
25.92
18.84
13.31
-0.11
Energy network
90.09
75.17
60.45
47.25
35.67
26.01
18.88
13.48
2.53
Table 2: Negative log likelihood (in bits/dimension) on ImageNer64 test set. See Appendix A. for more details.
Method
Anti aliasing
Augmentation
Discreteness
Type
Single NFE
NLL
Glow (Kingma and Dhariwal,
None
Contimuous
Normalized
PixcICNN (Van den Oord al., 2016)
X
None
Diserete
Normalized
X
3.57
I-DDPM (Nichol and Dhariwal, 2021)
X
None
Continoous
Upper bound
X
3.54
VDM (Kingma al., 2021)
X
Nonr
Discrete
Upper bound
X
3.40
FM (1 (pmman et al, 2023)
None
Uniform
Normalized
X
3.31
NEDM (Bartosh al., 2024)
Horizontal llips
Uniform*
Normalized
TarFlow (Zhai et al, 2024)
X
Horizental flips
Uniform
Normalized
2.99
Ours
Hotizontal Rlpps
Continoous
Estimate
3.36
Negative log likelihood. How accurate are our energy estimates? A standard evaluation proba- bilistic models consists of estimating the cross-entropy [ log po (2:) between the data distribution and the model, also known as negative log likelihood (NLL). In practice, NLL is computed by com- puting the probability the model assigns held-out test set. There are however several subtleties that make direct model comparisons challenging. Specifically, there are 3 factors that can cause variations in NLL of =0.5 bits/dimension more: details of data pre-processing (e.g., downsampling or data augmentation). conversion method from continuous to discrete probability, and the estimator type (exactly normalized, variational bound, approximately normalized). We expand on these issues in Appendix A.
We compare NLLs of our method and a variety of recent energy models in Table 2. This evaluation demonstrates that our model is comparable to the best-performing models the literature, within the variability arising from the three factors of variation mentioned in the previous paragraph. Two unique advantages of our method are that it provides (1) direct (one-shot) estimates of energy (as opposed to other density estimation approaches used with diffusion models, which we review in Appendix A.2), and (2) access to energy across all noise levels, providing a window into larger-scale features of the energy landscape. For instance, our energy network compute the probability of 50k images in 12s on an A100 GPU, whereas the score network requires upwards of 3h20, even with as few as 100 times steps and 10 noise samples compute the energy (see Appendix A.2). This justifies the longer training time due to the cost of double- buuble-backroopggtton (training for 1M steps on ImageNet64 a single A 100 GPU respectively took 120 hours, compared 32 hours for the score network), which may be alleviated by the use of sliced score matching (Song et al., 2020).
3 Analysis of the learned energy-based model
3.1
Generalization
The previous section demonstrated that our energy-based model achieves near state-of-the-art NLL on ImageNet64. That is, the model on average assigns high probability to set of held-out test images. Next, we establish that the energies the individual images are reliable. In particular, we verify that the model's energy assignment is stable under change of the training data.
To this end, we use the strong generalization test developed Kadkhodaie et al. (2024). We partition the training data into two non-overlapping sets, train a separate energy- based model on each set, and then compare the energies computed by these two models on images from both training subsets. We gradually increase the size of each training set until the two models assign approximately equal log probability across all images. Figure 2 shows the results of this experiment. The two models assign very different probabilities to the same image when the training set size, N. is small. But

---

100
1,000
10, 10,000
100,000
Thiased A
Tatsset
log
(15
(xg
Figure 2: Convergence of energy estimates. The data set is split into two halves (denoted A and B), and separate energy models are trained N samples drawn from each half. Each scatterplot compares the energy estimates of the two models at / = O, over all 2N training images. As N increases, the energy estimates of the models converges for all images.
10-1
Data
Gumbel fit
Prequency
10-2
10
10
40
60
(a)
(dB//im)
Figure 3: Histogram of log probabilities of test images in the ImageNet dataset. Color-coded arrows indicate values for the example images on the right, and the leftmost (brown) and rightmost (green) arrows indicate values for a uniform noise image [0, 1] and a constant image of intensity 0.5, respectively. The distribution is well-fit fit by Gumbel distribution (red line). Additional example images organized by probability are provided in Figures to 8 (Appendix D.1).
they converge gradually and compute nearly the same values at = 10 Note that the rate of convergence depends on image probability: more data needed before the models agree on the high-probability images. is also worth noting that the transition from memorization to generalization of energy models is marked by large increase of variance in energy the training set (starting N 100).
This result establishes that the model variance vanishes with feasible training set size. However, does not guarantee that the values are accurate- the models could be biased. Direct calculation of model bias requires access the true density, which is only available for synthetic data (as in Figure l). However, the small NLL values over test data (Table 2) do suggest a small model bias.
3.2 Distribution of energies and relationship to image content
We now study propertics of the learned energy model. What is the entropy (average energy) of the image dataset? Do all images have the same probability? If not, what determines the probability of an image? To investigate these questions, we compute the log probability of all 50, 000 images in the ImageNet64 test set. We express log probability in units of decibels per dimension (dB/dim), computed as 2/2 10 log Po (x), with d = 64 x 64 4096. In this scale, an increase of 10 dB/dim corresponds to a multiplicative rescaling of the probability density by factor of 10d.
Entropy. The average value of - loo log))(z) provides an estimate of 11..4 dB/dim for the differ- ential entropy of ImageNet. Note that the uniform distribution over the hypercube [0, 1]d has an entropy of This indicates that natural images occupy a fractional volume about 1.14. This can be converted to an estimate of discrete entropy by assuming that log probability is constant within quantization bins. For 8-bit images (with 256 possible intensity values), this corresponds to an entropy of 4.20 bits/dimension (the deviation from the value in Table 2 is due to the use of grayscale images here). In other words, there are 105.180 quantized ImageNet images out of 109. possible images this resolution.

---

1.0
30
0.8
40
Log prolubillgy (4B3 dimal
20
0.6
30
Imng (az
300
20
20
0.4
10
00
0.2
0.0
10-
5
02
0.6
0.8
1.0
0.2
Putteipation natio
Figure 4: Influence of image statistics on probability. Left. log (a as a function of and b. Middle. Horizontal slice (b- = 2) of the left panel for the example images of Figure 3. Right. Log
probability function of sparsity, measured as the participation ratio of wavelet coefficients.
Lack of concentration. Many high-dimensional probability distributions exhibit a concentration phenomenon: typical realizations have nearly equal energy. For instance, in simple image probability models such as a Gaussian model or a sparse wavelet model (where wavelet coefficients are indepen- dent), the energy of independent components add, and the law of large numbers imply concentration of the energy around its mean. contrast, our energy network reveals enormous diversity in the log probabilities of individual ImageNet images (Figure 3, left). Over the entire test set, they span a range of 34.4 dB/dim, corresponding to a probability ratio of x 10]4,000, Given this enormous ratio, it is surprising that any low probability image appears in finite dataset. To make sense of this apparent contradiction, it is important to realize that the total probability mass of images with given value of log p dataset is equal to this probability value multiplied by the volume of the corresponding probability level set. Observing a range of probability values of 1024,000 thus reveals that these enormous variations in probability density must be nearly compensated for by inverse variations in the volume of their corresponding level sets.
Shape of the distribution of log probabilities. The distribution of log probability values is highly skewed (Figure 3, left): there is heavy tail of high-probability images and a much lighter tail of low- probability images. We also observed these qualitative features on the CelebA and CIFAR-10 datasets, shown Appendix D.2. Surprisingly, this distribution is well approximated by Gumbel distribution (we report parameter values in Appendix C.3), which arises as a limit distribution for the maximum of many i.i.d. random variables. This surprising observation is, to the best of our knowledge, novel, and calls for an explanation. In independent component models, the log probabilities are sums of i.i.d. random variables, and are therefore caussian-distribtter model that leads to Gumbel distribution is high-dimensional sphericallyssymmedt distribution with an exponentially-distiibee
high dimensions. In contrast, a simple
radial marginal (Lyu and Simoncelli, 2009).
Image content. We also show in Figure (right) example images with a variety of log probabilities. High-probability images contain small objects, a blank (often white) background. Conversely, low-probability images are filled with dense detailed texture Since the energy of an image, when expressed in bits, corresponds to the size of its optimally compressed representation, this implies that nigher-probability images have shorter codes and thus less "information content". We show in Appendix D.4 images that have been obtained by locally maximizing minimizing probability.
Intensity range and sparsity. Following these visual observations, further examine the influence of intensity range and sparsity on image probability (Figure 4). First, we evaluate the log probability of reference images as we manipulate their brightness and contrast through an affine operation. We find that the brightness has minimal effect on the probability, while higher-contrast images have lower probability. This reveals that the distribution support is star-shaped: any pair of images are connected by high-probabblity path passing through a constant image. In contrast, we show in Appendix D.3 that linear interpolation between pairs of images typically increase probability: the distribution support is not convex. Next, we verify that log probability is correlated with simple measure of sparsity of images. We compute a multi-scale wavelet decomposition the images, and measure the (--norm divided by the (--oorm of the coefficients. The square of this quantity. ||>11/// / E (0, 1] known as the participation ratio, with smaller values indicating higher sparsity. The right panel of Figure shows that this simple measure captures significant portion of the variance of log po(x).
8

---

3.3 Effective dimensionality of the energy landscape
Beyond estimating log probability for a given image, we aim to characterize the local behavior of the density in the vicinity of that image. For instance, we would like to assess whether the probability is locally concentrated near a low-dimensional "tangent" subspace, and if so, to estimate its dimensionality. In this section, we explain how these quantities may be computed from our energy model. We then show that the local dimensionality around an image is often much lower than the ambient dimensionality of the space, but that this dimensionality depends on both the particular image and the scale that is used to define the local neighborhood.
Multi-scale dimensionality. Consider hypothetical distributions supported on the blue regions in the left two panels of Figure 5. For the left example, the effective dimensionality of a local neighborhood decreases with the size of that neighborhood. The opposite behavior is also possible, as shown in the right example. These examples demonstrate that dimensionality measures which aim to describe these geometrical structures need to depend on both the location and the scale of the neighborhood (Mohan* * et al., 2020; Tempezyk et al., 2022).
Effective dimensionality. We now introduce an effective dimensionality measure which can be equivalently defined from the optimal denoiser or the evolution of the probability landscape as diffuses. Given noisy observation of clean image x with noise variance I, the optimal denoiser estimates z with the conditional expectation R[xx| y]. Its average deviation from æ capture the local support of the data distribution around a: at scale t. Intuitively, from observing y, the optimal denoiser identifies that x is located on a dett-dimeniooal "tangent" space and projects y onto it. This preserves the components of the noise that lie along the subspace, incurring a denoising error tdott. Thus, we define the local effective dimensionality around x at scale t as
1 dot () , Ey || Ex wlll
(1))
Effective dimensionality can be equivalently defined directly from the energy. Consider the forward diffusion which progressively adds more noise to an image x, blurring the probability landscape. At time t. we can define an effective energy [y, t) æ]. Its rate of change with captures the effective dimensionality around x. Intuitively, if this landscape is locally a dest-dimensional subspace, then adding more noise causes probability to diffuse in the d dest normal "off- "off-manifold" directions, spreading over a volume ~ /(d-dr)/2 Taking the logarithm, we obtain Egy [U(y xx] ~ (d - dett) 1/ log L. Thus, we can equivalently define
des (zr, 1) = d - 0 log r FBg (y, 1) [=]]
(12)
We show in Appendix B.3 that the definitions in eqs. (11) and (12) are equivalent, unifying two seemingly distinct points of view underlying these dimensionality measures (Wu and Verdú, 201 I; Mohan et al., 2020; Tempezyk et al., 2022; Stanczuk et al., 2022; Kamkari et al., 2024; Horvat and Pfister, 2024).
Numerical results. We show in Figure 5 the behavior of the effective energy and dimensionality as function of the noise level. As the noise level increases, the log probability initially remains constant, indicating that probability is uniformly spread in a ball around x. The probability eventually decreases as the diffusion radius becomes larger than the local support size, and all lines converge to the mean (the negative entropy), consistent with the asymptotic Gaussian behavior of the energy. The effective dimensionality (fig. 5, right) vanishes at large noise levels (images have compact support, which behaves like a point at sufficiently large noise levels) but increases as the noise level reduced, eventually reaching the ambient dimensionality d (our model has continuous non-zero density everywhere by construction). In-between these two extremes, there exists sizable range of scales E [10-9, 10 where almost the entire range of dimensionalities coexist. In particular, images have diveediimeiionall neighborhoods, even at relatively small scales, while lower- probbbiibabilityy images have nearly dimenniioii neighborhoods, even at relatively large scales. This result empirically confirms the point raised in Section 3.2: the local probability mass around high and low probability images are of the same order, despite the enormous gap between their probability values. The local high dimensionality of the low-probability images corresponds to a significantly larger volume occupied by their local neighborhood. We note two potential limits to these estimates. First, the quantization of pixel values into bins of size 1/256 limits the resolution
9

---

Elllggeis: (e(s:)) [+]
(at(s.))
1000
Los probeaiiyyy (dl2)dmn))
2400)
2000
-10
10 10
Nidoo
Figure 5: Left: Two hypothetical examples illustrating how local effective dimensionality depends on the scale of the neighborhood. For both examples, support of density corresponds to the blue regions. In the left example, the dimensionality around the red point decreases with scale (from 2 down to 0), while the opposite true for the right example. Right: Log probability and effective dimensionality as function of noise level. Colored lines correspond to different example images I (shown in right panel of Figure 3), while the dashed black line shows the average over the ImageNet test set. The vertical gray line indicates the minimum noise level presented during training (t = 10 10-9), and the horizontal gray line the ambient dimensionality of the dataset (d = 4096).
to 10 possibly explaining the momentary decrease in dimensionality decreases. Second, our energy model was only trained on noise levels down to mmin 10-9. Thus, energies at smaller values of are solely determined by model inductive biases, which favor a constant energy.
Discussion
We have developed novel framework for estimating the log probability (energy) of a distribution of images from observed samples. The estimation method is simple and robust, and converges to stable solution with relatively small amounts of data (in our examples, 100, 000 images of size 64 x 64 suffices). The framework leverages the tremendous power of generative diffusion models, relying on networks trained to estimate the score of the data distribution by minimizing denoising objective. We augment this with secondary objective that ensures consistency of the model across noise levels. and an architecture constructed from an existing denoiser so to preserve its inductive biases. We validate our model by verifying that it achieves denoising performance equal to or surpassing the original denoiser, and NLL values comparable to the state of the art. We note that it is straightforward to extend our approach to conditional models: given conditioning information, the network needs only to compute the conditional energy function. The expression of losses (DSM and fTSM would be unchanged, as well as the energy architecture.
We mention two important limitations of our work. First, our energy model has a roughly quadrupled training time compared to the base score network, due to the use of double back pack-propagation. We believe this can be improved through architecture improvements, specialized auto-differenniaatin functions, or through the use of sliced score matching (Song et al., 2020). But in any case, the additional training time is offset by the dramatic efficiency gains in computing energy values. Another limitation of our approach is that we cannot provide any theoretical guarantees that minimizing our objective leads to good approximation of the true energy. We expand on this question in Appendix A.3.
Our learned energy model provides new insights into the geometrical properties of the space of photographic images. Notably, we observe lack of concentration of measure, with log probability values varying over a wide range and following a Gumbel distribution. As a limit distribution for the maximum, rather than the sum, of i.i.d. random variables, we expect that this might arise in broad class of digh-diiensional distributions. We also introduced a novel image- and scale- dependent measure of effective dimensionality which unifies two separate points of view developed in the literature. We demonstrate that different image neighborhoods display both high and low dimensionality over a wide range of scales. These results challenge simple interpretations of the manifold hypothesis. Moreover, the estimated dimensionalities remain too high to explain how the curse of dimensionality is avoided. This suggests that the energy landscape of natural images has additional undocumented geometrical regularity.
10

---

Acknowledgments
FG thanks Louis Thiry for starting this research direction (Thiry and Guth, 2024). We also thank Joan Bruna and Beerre-Etienne Fiquet for inspiring discussions. We gratefully acknowledge the support and computing resources of the Flatiron Institute (a research division of the Simons Foundation).
References
Grigory Bartosh, Dmitry Vetrov, and Christian Andersson Naesseth. Neural flow diffusion models: Learnable forward process for improved diffusion modelling, Advances in Neural Information Processing Systems, 37:73952- 2024.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. IEEE transactions pattern analysis and machine intelligence, 35(8):1798-188) 18288 2013.
Sagnik Bhattacharya, Abhiram R Gorle, Ahmed Mohsin, Ahsan Bilal, Connor Ding, Amit Ku- mar Singh Yadav, and Tsachy Weissman. ItDPDM: Information- theorrttc discrete poisson diffusion model. arXiv preprint arXiv:2505.0500 2025.
Giulio Biroli, Tony Bonnaire, Valentin De Bortoli, and Marc Mézard. Dynamical regimes of diffusion models. Nature Communications, 15(1):9557, 2024.
Joan Bruna and Jiequn Han. Provable posterior sampling with denoising oracles via tilted transport. Advances in Neural Information Processing Systems, 37:82863 0882663-22894, 2024.
Chen-Hao Chao, Wei-Fang Sun, Bo-Wun Cheng, and Chun-Yi Yi Lee. On investigating the conservative property of score-based generative models. International Conference on Machine Learning, 2023.
Kristy Choi, Chenlin Meng, Yang Song, and Stefano Ermon. Density ratio estimation via infinitesimal classification. In International Conference on Artificial Intelligence and Statistics, pages 2552- 2573. PMLR, 2022.
Patryk Chrabaszez, Ilya Loshchilov, and Frank Hutter. downsampled variant of imagenet as an alternative to the CIFAR datasets. arXiv preprint urXiv:1707.0885 2017.
R Cohen, Y Blau, D Freedman, and E Rivlin. It has potential: Gradient-drives denoisers for convergent solutions to inverse problems. Ady Neural Information Processing Systems (NeurIPS), 34, 2021.
Laurent Dinh, David Krueger, and Yoshua Bengio. NICE: Non-linear independent components estimation. arXiv preprint 111.v:1410.8516, 2014.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In International Conference on Learning Representations, 2017.
Yilun Du, Conor Durkan, Robin Strudel, Joshua B Tenenhaum, Sander Dieleman, Rob Fergus, Jascha Sohl-Dickstein, Arnaud Doucet, and Will Sussman Grathwohl. Reduce, reuse, recycle: Composi- tional generation with energy-based diffusion models and meme. In International conference on machine learning, pages 8489-8510. PMLR, 2023.
Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Society, 1993.
Dongning Guo, Shlomo Shamai, Sergio Verdú, et al. The interplay between information and estimation measures. Foundations and Trends@ in Signal Processing, 6(4):243-422) 2013.
Michael Gutmann and Aapo Hyvärinen. Noise-contrastive estimation: new estimation principle for unnormalized statistical models. In Proceedings of the thirteenth international conference artificial intelligence and statistics, pages 297-304. JMLR Workshop and Conference Proceedings, 2010.
Sébastien Herbreteau, Emmanuel Moebel, and Charles Kervrano. Normalization-equivari neural networks with application to image denoising. Advances Neural Information Processing Systems, 36:5706-5728 2023.

---

Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8)1171-1800,, 2002.
Geoffrey E Hinton, Terrence Sejnowski, et al. Learning and relearing in Boltzmann machines.
Parallel distributed processing: Explorations in the microstructure of cognition, 1(282-317):2 1986.
J Ho, A Jain, and P Abbeel. Denoising diffusion probabilistic models. Adv Neural Information Processing Systems (NeurlPS), 33, 2020.
Jonathan Ho. Xi Chen, Aravind Srinivas, Yan Duan, and Pieter Abbeel. Flow++: Improving flow- based generative models with variational dequantization and architecture design. In International conference on machine learning, pages
PMLR, 2019.
Christian Horvat and Jean-Pascaa Pfister. On gauge freedom, conservativity and intrinsic dimen- sionality estimation in diffusion models. The Twelfth International Conference on Learning Representations, 2024.
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Gradient step denoiser for convergent plug-and-play. arXiv preprint arXiv:2110.03220, 2021.
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statisties-Siullation and Computation, 1856) 1989.
Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Z Kadkhodaie and E P Simoncelli. Solving linear inverse problems using the prior implicit in denoiser. arXiv preprini arXiv:2007. 13440, Jul 2020.
Zahra Kadkhodaie, Florentin Guth, Eero P Simoncelli, and Stéphane Mallat. Generalization in diffu- sion models arises from geometry-adaptin harmonic representations. In The Twelfth International Conference on Learning Representations, 2024.
Hamid Kamkari, Brendan Ross, Rasa Hosseinzadeh, Jesse Cresswell, and Gabriel Loaiza-Ganem. geometric view of data complexity: Efficient local intrinsic dimension estimation with diffusion models. Advances in Neural Information Processing Systems, 37:38307-38354, 2024.
Rafat Karczewski, Markus Heinonen, and Vikas Garg. Diffusion models as cartoonists! The curious case of high density regions. In The Thirleenth International Conference on Learning Representations, 2025.
Diederik Kingma, Tim Salimans, Ben Poole, and Jonathan Ho. Variational diffusion models. Advances in neural information processing systems, 34:21696-21707, 2021.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions. Advances in neural information processing systems, 31, 2018.
Erich Kobler and Thomas Pock. Learning gradually non-convex image priors using score matching. arXiv preprint arXiv:2302. 10502, 2023.
Frederic Koehler, Alexander Heckett, and Andrej Risteski. Statistical efficiency of score matching: The view from isoperimetry. ICLR, 2023.
Xianghao Kong, Rob Brekelmans, and Greg Ver Steeg. information-theoretii diffusion. In The Eleventh International Conference Learning Representations, 2023.
Chieh-Hsin Lai, Yuhta Takida, Naoki Murata, Toshimitsu Uesaka, Yuki Mitsufuji, and Stefano Ermon. FP-diffusion: Improving score-based diffusion models by enforeing the underlying score fokker-planck equation. In International Conference Machine Learning, pages 18365-18398. PMLR, 2023.
Yann LeCun and Fu Jie Huang. Loss functions for discriminative training of energy-based models. In International workshop artificial intelligence and statistics, pages 206-213. PMLR, 2005.
12

---

Yann LeCun, Sumit Chopra. Raia Hadsell, M Ranzato, Fujie Huang, et al. A tutorial on energy-based learning, Predicting structured data, 100) 2006.
Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In The Eleventh International Conference Learning Representations, 2023.
S Lyu and E E Simoncelli. Nonlinear extraction of 'independent components" of natural images using radial Gaussianization, Neural Computation, 2166):1485-1519, Jun 2009, doi: 10.1162/meco.0009 04-08-773.
K Miyasawa. An empirical Bayes estimator of the mean of normal population. Bull, Inst. Internat. Statist., 38:18 181-188 1961.
S Mohan*, 2 Kadkhodaie*, E P Simoncelli, and C Fernandez-Granda. Robust and interpretable blind image denoising via bias-free convolutional neural networks. In Int'l Conf on Learning Representations (ICLR), Addis Ababa, Ethiopia, Apr 2020.
Eric Nalisnick, Akihiro Matsukawa, Yee Whye Teh, Dilan Gorur, and Balaji Lakshminarayanan. Do deep generative models know what they don't know? In International Conference on Learning Representations, 2019.
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International conference on machine learning, pages 8162-8171. PMLR, 2021.
M Raphan and E P Simoncelli. Least squares estimation without priors or supervision. Neural Computation, 23( Feb 2011. doi: 10. 162/NCCO 00076.
Gabriel Raya and Luca Ambrogioni. Spontaneous symmetry breaking in generative diffusion models.
Advances in Neural Information Processing Systems, 36:66377-66389, 2023.
Edward T Reehorst and Philip Schniter. Regularization by denoising: Clarifications and new interpretations. IEEE transactions on computational imaging, 5(1):52-67, 2018.
D Rezende and S Mohamed. Variational inference with normalizing flows. In Ins'l Conf on Machine Learning (ICML), pages 1530-1538. PMLR, 2015.
H Robbins. An empirical bayes approach to statistics. In Proc Third Berkeley Symposium on Mathematical Statistics and Probability, volume I, pages 157-163. University of CA Press, 1956.
Yaniv Romano, Michael Elad, and Peyman Milanfar. The little engine that could: Regularization by denoising (red). SIAM Journal on Imaging Sciences, 10(4): 184-)88444 2017.
Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, al. Imagenet large scale visual recognition challenge. International journal of computer vision, 115:211-2522 2015.
Tim Salimans and Jonathan Ho. Should EBMs model the energy or the score? In Energy Based Models Workshop-ICLR 2021, 2021.
Saced Saremi, Arash Mehrjou, Bernhard Schölkopf, and Aapo Hyvärinen. Deep energy estimator networks. arXiv preprint arXiv: /80.008306 2018
Marta Skreta, Lazar Atanackovic, Jocy Bose, Alexander Tong, and Kirill Neklyudov. The superposi- tion of diffusion models using the Itô density estimator. In The Thirteenth International Conference on Learning Representations, 2024.
J Sohl-Dickstein, E Weiss, N Maheswaranathan, and S Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. Francis Bach and David Blei, editors, Proc 32nd Int'l Conf on Machine Learning (ICML), volume 37 of Proceedings of Machine Learning Research, pages 2256-2265 Lille, France, 07-09 Jul 2015. PMLR.
Song and S Ermon. Generative modeling by estimating gradients of the data distribution. Adv Neural Information Processing Systems (NeurIPS), 32, 2019.
13

---

Y Song, Sohl-Dickstein, P Kingma, A Kumar, S Ermon, and B Poole. Score-based generative modeling through stochastic differential equations. In Int'l Conf on Learning Representations (ICI.), 2021a.
Yang Song and Diederik P Kingma. How to train your anergy-based models. arXiv preprint arXiv:2101.03288, 2021.
Yang Song. Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: scalable approach to density and score estimation. In Uncertainty in artificial intelligence, pages 574-584. PMLR, 2020.
Yang Song, Conor Durkan, lain Murray, and Stefano Ermon. Maximum likelihood training of score-based diffusion models. Advances in Neural Information Processing Systems, 34:1415-1428, 2021b.
Stanczuk, Georgios Batzolis, Teo Deveney, and Carola-Bibiane Schönlieb. Your diffusion model secretly knows the dimension of the data manifold. arXiv preprint arXiv:2212.1200 2022.
Piotr Tempezyk, Rafat Michaluk, Lukasz Garncarek, Przemystaw Spurek, Jacek Tabor, and Adam Golinski. LIDL: Local intrinsic dimension estimation using approximate likelihood. In Interna- lional Conference on Machine Learning, pages 21205-21231. PMLR, 2022.
Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319 2323, 2000.
Theis, A van den Oord, and Bethge. A note on the evaluation of generative models. In International Conference Learning Representations (ICLR 2016), pages 1-10, 2016.
Louis Thiry and Florentin Guth. Clasiffcatioon-deniiiing networks. arXiv preprint arXiv:2410.03505. 2024.
James Thornton, Louis Béthune, Ruixiang ZHANG, Arwen Bradley, Preetum Nakkiran, and Shuangfei Zhai. Composition and control with distilled energy diffusion models and sequen- tial monte carlo. International Conference on Artificial Intelligence and Statistics, pages 3259-3267. PMLR, 2025.
Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Instance normalization: The missing ingredient for fast stylization. arXiv preprint arXiv:1607,0800 2016.
Aaron Van den Oord, Nal Kalchbrenner, Lasse Espeholt, Oriol Vinyals, Alex Graves, et al. Conditional image generation with PixeICNN decoders. Advances in neural information processing systems, 29, 2016.
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
P Vincent. A connection between score matching and denoising autoencoders. Neural Computation, 23(7):1661-1674 2011.
Martin Wainwright. High-dimensional statistics: non-asymptotic viewpoint, volume 48. Cam- bridge university press, 2019.
Yihong Wu and Sergio Verdú. MMSE dimension. IEEE Transactions on Information Theory, 57(8):
4857-4879, 2011.
Shahar Yadin, Noam Elata, and Tomer Michaeli. Classification diffusion models. arXiv preprint 200.0 2024.
Jonathan S Yedidia, William T Freeman, and Yair Weiss. Constructing free-energy approximations and generalized belief propagation algorithms. IEEE Transactions on information theory, 51(7): 2282-2318 12, 2005.
Hanlin Yu, Arto Klami, Aapo Hyvarinen, Anna Korba, and Omar Chehab. Density ratio estimation with conditional probability paths. In Forty-second International Conference on Machine Learning, 2025.
14

---

Shuangfei Zhai, Ruixiang Zhang, Preetum Nakkiran, David Berthelot, Jiatao Gu, Huangjie Zheng, Tianrong Chen, Miguel Angel Bautista, Navdeep Jaitly, and Josh Susskind. Normalizing flows are capable generative models. arXiv preprint arXiv:24/2.06322, 2024.
Mingtian Zhang, Oscar Key, Peter Hayes, David Barber, Brooks Paige, and Prancois-Xavits Briol. Towards healing the blindness of score matching. In NeurIPS 2022 Workshop on Score-Based Methods, 2022.
Kaiwen Zheng, Cheng Lu, Jianfei Chen, and Jun Zhu. Improved techniques for maximum likelihood estimation for diffusion ODEs, In International Conference on Machine Learning, pages 42363- 42389. PMLR, 2023.
15

---

A Discussions on density estimation
A.1 On comparison of NLL values
We identify three issues that arise in estimation of NLL values, each of which limits the precision with which they can be compared.
A first issue concerns the selection and pre-processing of the training data. There are two different versions of ImageNet at 64 x 64 resolution: one introduced in Van den Oord al. (2016), no longer available, and a second which uses anti- anti-aliasing during downsampling (reducing entropy and thus NLL), introduced in Chrabaszez al. (2017). This change in the dataset accounts for variations of about 0.3 bits/dimension (Zheng et al., 2023). Similarly, the data augmentations used can both increase entropy (e.g., random flips) or reduce it (e.g., center crops, or resizing operations with anti-aliasing).
A second point is that NLL estimates must be computed on discrete probability model, and NLL estimates for continuous models are dependent on the method used to discretize the distribution. The simplest option is assume that the probability is uniform within quantization bins. Discrete probabilities are then computed as a product of the continuous probability with the volume of the quantization bin, which corresponds to an additive shift of 8 bits in the NLL for image data. This can be enforced by adding uniform noise to the data, a technique known as "uniform dequantization" (Ho al., 2019). which leads an upper bound on the NLL of the corresponding discrete model (Theis et al., 2016). These differences can account for variations of 0.1 to 0.5 bits/dimension ( (Kong et al., 2023). A more sophisticated variational dequantization procedure can improve results by an additional 0.1 bits/dimension (Ho et al., 2019; Song et al., 2021b). In contrast, directly modeling the discrete data with appropriate methods can lead reductions in NLL of more than 2 bits/dimension (Bhattacharya et al., 2025).
Finally, a third point is that depending on the method, the NLL should be interpreted differently. Classically, Pe is a normalized probability distribution. typically obtained through a change-of- variable formula (as in flow-based models), so that up to an unknown additive constant (the entropy of the data), the NLL directly evaluates the KL divergence KL(p || pe). For some other models (such as VAEs), the NLL is not directly tied to probabilistic model and is instead evaluated through a variational lower bound, leading to upper bounds on p(z) for each z. The NLL. however remains an upper bound on the entropy of the data. The MSE-based formulas of Kong et al. (2023) also fall in this category when using (necessarily suboptimal) network denoiser. Lastly, approximately normalized energy-based models such as ours or CDM (Yadin et al., 2024) compute estimates of - log p(z) that are neither lower nor upper bounds, so that NLL can only be interpreted as an estimate of the entropy of the data.
As result of these differences, NLL values of different methods are often not directly comparable. Here, we aim to learn an accurate continuous probability model of image distributions rather than obtaining the best upper bound on the discrete entropy of the data.
When evaluating previously published NLL results, we found that these details were often not provided, and we thus had to make assumptions in compiling Table 2. The "anti-aliasing" column refers to the version of ImageNet64 used, we assumed that articles that cite Van den Oord al. (2016) for the dataset made use of the aliased version. We assume no data augmentation if none is mentioned (for TarFlow (Zhai et al., 2024), the authors mention performing center crops of the data, but an inspection of their code indicates that it has been disabled for ImageNet64). The "discreteness" column refers to the nature of the probability model and the potential conversion from continuous to discrete probabilities ("discrete" for discrete probability models, "continuous" for an additive shift of bits, and "uniform" for uniform dequantization by adding uniform noise). We assume a continuous model if no dequantization is mentioned (FM (Lipman et al., 2023) and NFDM (Bartosh et al. 2024) use slightly different notion of uniform dequantization which improves NLL by 0.04 bits/dimension). The "type" column refers to the nature of the LLM estimate ("normalized" for an exactly normalized probability model coming from a flow- model, "upper bound" for variational lower bounds on log probability, and "estimate" for other approaches). Finally, we indicate those methods that use single neural function evaluation (NFE) to compute log probability the "single NFE" column.
16

---

A.2 Background on density computation in diffusion models
We review several methods of estimating log probabilities from diffusion models that have appeared in the literature.
In early publications on diffusion methods, probabilities are computed with the so-called probability flow ODE (Song et al., 202 This corresponds to the distribution of samples generated by the backward ODE. Given large noise level maax the backward ODE solves the equation
dæ dt
VoUo(EE,))
(13)
backwards in time from 2 tmax (0, mmax and produces an approximate sample x 20 PODE. The log probability of this sample can be calculated as (Song et al., 2021a)
||0ttmaxx| 11 + m2g(2x))mx)
d
mmax
log PODE(L)
ApUo(11, t)dt.
(14)
2fmax
Note that eq. (14) is also valid for arbitrary test points x; in which case the ODE (13) needs to be solved forward forward in time from Ho = x at = O to = mmax-- Equation (14) requires estimating the divergence of the score (the Laplacian of the energy), which is typically approximated with the Hutchinson trace estimator (Hutchinson, 1989).
Another approach is to use variational bound as in Song et al. (2021b); Kingma al. (2021); Kong et al. (2023). This variational bound arises from the exact identity (Kong et al., 2023)
(td (td-Eyy Ey || - llog = EylU [T(( (max) | 20| f.
dt
E, lvll x
262 (15)
The first term is the effective energy at mmxx. which is equivalent The second term features the optimal mean squared error over noise levels / E [0, mmxx] Note that the integrand can be rewritten as (d deu (x., ())d log t). and eq. (15) can thus be derived by integrating the equivalent definitions of effective dimensionality (11) and (12) (see Appendix B.3). Equation (15) naturally leads to an upper-bor souur-bound on the negative log probability of the data when replacing the optimal denoiser Ex [xr y, d] with the denoiser derived from the Tiwese expression, vV, U,,,,, t):
log((2re/ as mmxx o0.
log PMSE (x) log max)) -J l
d
x
dt 212
(16)
2
This framework can be generalized to other noise distributions (Guo et al., 2013) such as Poisson noise, which more adapted to discrete distributions (Bhattacharya et al., 2025).
Finally, it was recently observed in Skreta et al. (2024); Karczewski et al. (2025) that cheaper unbiased stochastic estimator of this bound can be obtained with the Itô formula. Given realization (24))ER, of the forward SDE (Brownian motion) dr, dw, started at
p(x). then
Fanas (V)U(((,))).d))) -
1 [ U(U d)((at) 2
- log p(x) = U (1 (mise tmax)
(17)
I 112 + 2 log( (0g(22/mmx)) Replacing the 24maae
When maxx + one can again use U ( (max) unknown true energy U in the integrand with model Up leads to biased stochastic estimate of the log probability:
| 11? d
(fmax //m (IV,Uo(zn,)))da)) t), dz,) - || - (V, fI]
log (2DE(x)
log log(2rtmax))-
2bmax 2
(18)
Karczewski al. (2025) show that averaging SDE trajectories started at the same xe = x leads to an upper bound on the NLL, in fact recovering the one in eq. (16).
In summary, the ODE gives deterministic probabilities exactly corresponding to the corresponding generative model (eq. (14)), while for the SDE one can choose between a denoising- upper bound (eq. (16)) cheap stochastic estimator of it (eq. (18)). These evidently similar approaches
17

---

are equivalent if and only if the model is consistent across noise levels (i.e., satisfies the diffusion equa- tion). In practice. the choice of estimation method can lead to variations of up to 0.5 bits/dimension (Kong et al., 2023), see in particular Karczewski et al. (2025) for careful study). These approaches also apply to our model, in addition to the direct evaluation of Uo( I = 0). While this latter estimate does not correspond to generative model (like the ODE) or a variational bound (like the SDE), its main advantage is its computational efficiency, as it does not require integrating over noise levels and gives deterministic values without averaging over noise realizations (for divergence term or mean squared error).
A.3 On the theoretical justification of dual score matching
A limitation of our approach is that we offer no theoretical guarantees that minimizing our objective leads to a good approximation of the energy. It would be desirable to show that our training objective (assuming infinite training data) quantitatively controls the distance between the learned energy and the data energy. This could be achieved by showing that the joint distribution ply, t) (with log t uniformly distributed in [log mini log (max) satisfies Poincaré inequality, i.e.,
Var(Us(y.t) Uv, t) S CE CI
VDU t)ll (0,Uolu,t) t))) (19)
for some constant C' that is not too large. While we conjecture that such result holds (or a variant, such as when considering the distribution of (y/ /f, log to match our noise level weighting), note that this weaker than a control in the Kullback-Leible divergence. Replacing Var Var t) Uy, t)] with (L Il Po) in eq. (19) would require that the model distribution (v t) satisfies a Sobolev inequality, which could only be enforced with specific architectures. As result, there is no control of the learned energy outside the support of the data distribution, so that out-of-distriuttion detection may be unreliable. It also implies our model may not be exactly normalized. Our NLL calculations are thus only approximate (although computationally cheap). These limitations are common to current score-based and unnormalized energy- based approaches, including the related work of Yadin et al. (2024).
A.4
On frequency estimation with density models
Here, we describe counterintuitive property of density models in high dimensions which poses a challenge for estimating frequencies. We also refer the interested reader to Theis et al. (2016), which offers related observations.
Consider a mixture of two uniform distributions on two compact sets (classes) C₁ and C2 with respective frequencies fi and f2. The probability density is then constant on each class, with value d where V42 is the volume of C. The volume Vi typically scales exponentially with d, V ~ r where rę is the characteristic size of C,. The energy values on each class then
p: V.
U log p; dlog r;
log fi.
(20)
In high dimensions d >>>> 1, the energy values are dominated by the volume and only weakly depend on the frequencies of each class.
A concrete example of this phenomenon can be seen in the mixture of two Gaussian distributions considered in Figure 1, where the two classes correspond to the two spheres of radius vdo; for i E {1, 2}. The typical values of the energy U, = log (2mea?) log Sa are dominated by the entropy of the corresponding Gaussian.
The implications of this observation for -ich-iminniional density models are twofold. First, estimated densities are dominated by volumes of typical sets (entropies of the mixture components). more than frequencies of different categories (one-dimensional marginals). (Note that the latter is easily learned from data, while estimating the former is much more challenging task.) This observation explains empirical reports that energy models may assign higher probability to nut-of-disttibution samples (Nalisnick et al., 2019; Karezewski et al., 2025). Second, estimating energy up to a small relative error is not sufficient to capture observed frequencies. For our ImageNet64 model (d 4096), a relative error of ~ 0.1 dB/dimension (as estimated from the generalization experiment in Figure 2) is negligible compared to typical energy differences of ~ 10 dB//imension, but the required precision to estimate frequencies with relative accuracy of 10% is 1/d 2 x 10-1 dB/dimension.
18

---

B Proofs and derivations
B.1 Time score matching (proof of eq. (4))
We start from the expression of the energy (negative log probability) of noisy data:
U(y () log ply(();
(21)
p(ylt) / p)d
(22)
Differentiating the energy with respect to I yields
dep(ylt)
(,U
(23)
p(yl()
J (pp))((pa t) plyla, t)dæ
(24)
plylt)
Ex log (ogg((y 1) 4.
(25)
Note that the derivation exactly matches that of the Tilyasaw-TTweede identity by replacing 0 with Vy, and does not make any assumptions about the form of p(ylæ, t). Restricting to additive Gaussian
noise, log plyl xz t) = 2:. lly ==l 2 1og(22t), and substituting into eq. (25) gives the "time score-matching" identity of eq. (4):
d
11
2=1
(y, Ex
y
(26)
2t
26
B.2 Energy architecture
Suppose that the energy network is defined in terms of network so(y, Uo(y.t) (y. sy t)), where the score network is assumed conservative and homogeneous. Conservativity means that there exists scalar function á such that so () Vyd(y, L), which implies that the Jacobian V (y,)(y,t) Vooy, t) is symmetric. The homogeneity property requires that for all O, so(Ay, t) (sa t). Differentiating with respect to À and setting X = 1 yields
sys so(y.t).
(27)
We now calculate the gradient of the energy network:
v v,Uo(.tt)
(so(v. t) + Vyso(u,t)² = so(y, t),
(28)
using the conservative and homogeneity properties to derive that Vyso(y, f)y V vyso(y,
so s)) Note that even if se is not conservative, then it still holds that (v,
((yyso Sy t) which can be interpreted symmetrization of the Jacobian of sg.
We also remark that sy homogeneous, then Up is quadratically homogeneous (Ua t) X ()). Note that this does not correspond enforcing (asymmetric) Gaussian one-dimensional marginals (y, u). Rather, this enforces that the distribution of (y, u) conditioned on the orthogonal projection - (y, u)//lla/ ///l is (asymmetric) Gaussian.
B.3 Effective dimensionality (equivalence between eqs. (11) and (12))
We start by calculating the time derivative of the effective energy Ey (U,UU( t) xx). We have
0,E (v,,, 0 / plula, t)U t)dy
(29)
/ (3.p(ula, ()U plyla, ()a,U f))dy.
(30)
19

---

The first term is the derivative with respect to variance of a Gaussian distribution (z idd. It satisfies the diffusion equation:
O,p(y/2, 8) - - 2
(31)
Similarly, for the second term, we use the fact that Ug..t) f) - log p(ylt) where p(ylt) also satisfies the diffusion equation (as can be seen from the Fokker-Planck equation the variance exploding case (Song et al., 2021a)): 0,U (v.1)
Deplult)
(32)
p(ylt) A (
(33)
2p(ylt)
2p(vlt) V. (P(Ul))V
(34)
- ,((l()))
1
1 [v
(35)
2
This equation appeared in Lai et al. (2023); Bruna and Han (2024) and is a special case of Hamilton- Jacobi equation (Evans, 1993).
Combining eqs. (31) and (35) into eq. (30) and integrating by parts twice, we have
(U,,, (z] plule, t) A, t) (V
(36)
We recognize the expression of the mean squared error as given by combining Miyasawa-TTeedle with Stein's 's unbiased risk estimate (SURE). Indeed,
Ey ||= E2 [x |vll
= E,
|| - y + tV UV, 1]
(37)
Ey || - y/|² + 21 y, V U (v, ()) + + , 0)1
(38)
= td + fE 22 U (g. t) + V t)]l [
(39)
where we have used Stein's 's lemma in the last step. We finally combine eqs. (36) and (39) to obtain
_
By [| [x [w]
x
d - B, (v, 1) [z]
(40)
It is convenient to rewrite 2t0:- = 0g logt as d( logt) dr.
We define the common value in eq. (40) to be den (z,t). It is related (but not equal) to dimensionality measures estimated from the singular values of the Jacobian of a denoiser (Mohan et al., 2020; Horvat and Pfister, 2024). The limit of dett t)) when -> O has appeared in the literature under the name of local intrinsic dimensionality using (approximate) likelihood (LIDL) (Tempezyk et al., 2022; Stanezuk al., 2022; Kamkari et al., 2024). Its average over images is equal to the MMSE dimension of Wu and Verdú (2011).
C Experimental details
C.1 Energy network architecture
UNet architecture. Our UNet architecture sa is composed of 3 encoder blocks, middle block, and 3 decoder blocks. Each block is itself composed of 3 layers, each a sequence of bias-free convolution, normalization, and non- non-linearity, for total of 21 layers. The first convolutional layer of each encoder block but the first and the middle block has a stride of 2 (downsampling, in both vertical and horizontal directions), and the last convolutional layer of each decoder block (except for the last one and the middle block) is transposed with stride of 2 (upsampling). The output of each encoder block concatenated to the input of the corresponding decoder block (which comes from the output of the corresponding encoder block, via a "skip connection") ").. The number of channels is doubled in each block, starting from a base value of 64 at the coarsest scale. We replace ReLUs with GeLUs to ensure that V,U, is differentiable. Thus, Up is only approximately quadratically homogeneous. As result, and because training does not result in an exactly conservative 59, VyUo(y, 1) should be computed directly as opposed to using s8.
20

---

Normalization, We also replace batch normalizations, whose behavior during the backward pass is incompatible with a second paok-propagttion, with homogeneous version of instance normalization (Ulyanov et al., 2016). If the input x consists of C' channels (se)isesC, and its spatial mean is (x(x) E R,, each channel is normalized according to
2 p plx))1
Ze C²I (( (2c u(ze)).
(41)
Up to a small e >> 0 parameter for numerical stability, this ensures that after normalization all channels ze have equal norms while preserving the global norm |>]|| This normalization layer also homogeneous when E = 0, as the spatial mean is estimated from the input z. This normalization is followed by learned rescaling of each channel, Ze + Ye where y E RC is learnable.
Noise level conditioning. The noise variance is also an input of sg. As is standard in diffusion models (Nichol and Dhariwal, 2021), time embedding e(t) 2256 is computed with Fourier features (os(t),), sin(w)t) (we use 32 frequencies (wx)) ) that are linearly spaced in the log domain and ranging from 1/fmux to 1 mmim followed by shallow MLP. This time embedding e(t) is then used to condition the output of each normalization layer via gain control: Ie + ((1 + (We, (())) where We a learned layer- and channel-dependent vector R 206
C.2 Training hyper- hyper--prameee parameaaaeeess
To summarize Section 2.2, the dual score matching training objective is
I VVvi 5a + ("ocuatun) -
2 1(( 1 ||=11² d ))
{(0) = Ex,,,
(42)
where x ~ p(z) is the data distribution, 2 ~ NO, Id) is the noise, y = x + vtz is the noisy measurement, and logt (log Euin log maxx) In our experiments we use tmin = 10-9 and 'max =103, and the training image intensities are rescaled to have values in [0, 1].
We use the ImageNet64 dataset (Chrabaszez et al., 2017), with (only) horizontal flips as data augmentation. The models used in Tables 1 and and for the generalization experiment in Figure 2 are trained on color images, while the model used in Figures 3 to 5 is trained on grayscale images. Pixel values rescaled to [0, 1] by dividing by 255. All models are trained for 1M steps, with batch size of 128. We use the Adam optimizer with default parameters and an initial learning rate of 0.0005 (except for the generalization experiments which used a learning rate of 0.0002) that is halved every 100, 000 steps. All models are trained on a single NVIDIA H100 GPU, which takes about 5 days for ImageNet64.
C.3 Additional details
Gaussian scale mixture example (Figure 1). We generate n 100, 000 samples from a mixture of two Gaussian distributions, [W(00, ofd) + W0 oîId), with 01 = 1 and 02 4, in dimension d 1, .00. The true (normalized) energy is given by
1nl 2( ワ
log 10s(2=(0++)))
U 1) - log Σ
(43)
We parameterize the energy as a mixture of quadratics:
Uo(y,t) = -log Σ (000l
(44)
where the functions a, 6; are computed by 5-layer MLP with hidden dimension of 256 that takes log(t: + fmin as input. This network is trained either with single (space) score matching or with dual score matching, both across noise levels & E [tmin, (max]. Training otherwise similar to the ImageNet64 experiments (see Appendix C.2), for 20, 000 training steps with batch size of 512 and an initial learning rate of 0.0001, over noise levels from 4min = 10 to /max 10 We note that the energy learned after single score matching training is stochastic: as the relative energy between the two mixture components is not constrained by the data, its value is determined by the random initialization, so that rerunning the experiment will lead to a different value.
21

---

Histogram of log probabilities (Figure 3). The Gumbel fit is calculated by maximizing likelihood. The obtained parameters (in decibels/dimension)) are 9.57 for the location and 3.17 for the scale. Equivalently, p(z) follows Fréchet distribution, with a scale parameter equal to 9.05 and shape parameter equal to 1.37/d, d, where d = 64 x 64 4096 is the dimension of ImageNet64 grayscale images.
Computing dimensionality (Figure 5). Equations (11) and (12) provide two ways to estimate effective dimensionality from learned energy model. If the model exact, Up(y. t) Ug. t), then they coincide (more generally, they coincide when the model satisfies the diffusion equation), but in general they only approximately coincide. For numerical stability, we found it preferable to use the version defined from the mean squared eπror of the underlying denoiser (eq. (11)). This guarantees non-negative dimensionalities, and also has the advantage of being an upper bound on the true dimensionality (as any denoiser yields an upper bound on the minimum MSE).
D Additional results
D.1 ImageNet images according their probability
Figure 6: Highest probability images in ImageNet64 (test set).
Figure 7: Lowest probability images in ImageNet64 (test set).
Figure 8: Example images from ImageNet64 (test set) with inearly-spaced log probabilities, ordered from low to high.
22

---

D.2 Distribution of log probabilities on other datasets
To evaluate the pervasiveness of the models CelebA (grayscale and color) and CIFAR-10, with identical architecture, objective, and hyper-parameters as described in Section and Appendix C. Results are presented in Figures 11. In all cases, log-probability histograms share qualitative features with the Gumbel distribution, namely its skewness towards high-probability values, although the fits are less precise due to the presence of low- or probbabbbbylity outliers.
log probability values, we trained energy
10-1
Frequency
10"2
10-3
Data Gumbel fit
10
5f f 110 1 115 f 12b f 125
logp,(x) (dB/dim)
Figure Histogram of log probabilities of test images in the CelebA dataset (grayscale). Color-coded arrows indicate values for the example images on the right. The maximum-likelihon fit by a Gumbel distribution (red line) is perturbed by the low-probabilityy outliers, some of which correspond to corrupted images.
10
Frequency
10-2
MAL
10-3
Data
10-4
Gumbel fit
110 f1 h3 20
125
30
log (dB/dim)
Figure 10: Histogram of log probabilities of test images in the CelebA dataset (color). Color-coded arrows indicate values for the example images on the right. The maximum likelihood fit by a Gumbel distribution (red line) is perturbed by the low-probability outliers.
maximum-likelibodd
10-1
Data Gumbel fit
Frequency
-
1
10
15
20
25
30
Po(p) (dB/dim)
Figure 11: Histogram of log probabilities of test images in the CIFAR-10 dataset. Color-coded arrows indicate values for the example images on the right. The maximum-likelihon fit by Gumbel distribution (red line) is perturbed by high-probability outliers, which correspond to grayscale images.
23

---

D.3
Log probabilities of interpolated images
10
Log probabiity of Interpolated image (B,dda))
35
30
25
15
10
Bweet
Interpolated log probability (dB/dim)
Figure 12: Log probability of interpolated images. For each pair (Ii, 2:441) of consecutive example images (see Figure 3), we consider linear path between them, a E [0, + a)x; ar (shown in the right panel). In the left panel, we plot the log probability of the resulting images log pep(p11 a);; + a1++1) as a function of the interpolated log probability a log (x(a) + a) log 1. For medium- and high- gigh-probability images, this leads to sharp decrease in prob- ability even for close or the superposition of two natural images not natural image. However, for low-probabilityy images, this surprisingly leads to slight increase in probability: the interpolation between two densely textured images is another densely textured image but with slightly lower contrast, leading to higher probability values (as shown in Figure 4).
D.4 Optimizing log probability
Log prbability (c3lddm)
(+.
(-.0
(4 74)
20.6 82)
(9 199
21. 0)
(+ 0)
32.5 (-
37.0 0)
47.5 (t 0)
Maximizing (x())
Initial image
Minimizing log (L)))
Log prolubiilty (dBAlin) -1.6(-0.))
(.76)
3.5(-44.)
(-3.)
49(-124)
4.6-17..0
496 22.3)
544-27.1) 5.6(-144
3-0(-4449)
For each initial image (middle row), we perform 1,000 steps of gradient ascent (top row) or descent (middle row) on the log probability with fixed step size 10 10-2/d, clipping pixel intensity values to [0, 1] at each step. Due to numerical instabilities around high-probability images, we select the iterate with the highest or lowest probability. Resulting log probability values, and their difference with the starting value, are shown above or below the final images. Maximizing log probability leads generally smoother images (although with some visual artifacts), except for sufficiently high-probabilityy images, where gradient ascent fails to increase the probability. Minimizing log probability leads to noisier images, with regular geometric patterns for high-probability images.
Figure 13:
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Yes, we describe the scope and motivations for the work clearly in the abstract and introduction sections.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We discuss limitations and assumptions extensively throughout the paper, particularly in the discussion.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
25

---

Answer: [Yes]
Justification: Assumptions are stated clearly and proofs are provided either in the appendix or by citation.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
Experimental details are provided in Appendix C.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
26

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have only used publicly available datasets (ImageNet64). We release code to reproduce all experiments and pre-trained models at https://github. com/ aareaniifuut/DDaalsccraaaaa...
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.c cc/ public/guides/CCn
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https:
/ /nips. cwtppbliccguiceeeee wwtcc//guies//eee http:/uuuiddeeooSSuubbiiii SrndeeSubiissioPolll..
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: Experimental details are described in Appendix C.2.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Statistical error is negligible for the reported results. Figure 2 verifies that the trained energy model is not sensitive to its initialization or subset of the training data.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
27

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
These details are reported in Section 2.4 and Appendix C.2.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: newrips.....//
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work is not about applications but about foundations energy estimation, which makes the short-term societal impact negligible.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
28

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: Both the original and downsampled ImageNet papers are cited.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URI.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
29

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/com//ttt/o datasttss has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper doess not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
30

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
We did use LLMs this work.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https:/ /ttpa/:////e....///// ///Cooffrenccee/2225//// for what should or should be described.
31


Reviews:
[{"review_id": "H67yGJ2jU9", "reviewer": "Reviewer_Yarf", "summary": "This paper presents a framework for learning normalized energy-based models directly from image data, leveraging techniques from score-based diffusion models. Particularly, a dual score matching objective is proposed to match gradients of the energy network with respect to both space and noise level. The authors further propose a novel architecture for energy estimation by computing the inner product of an input image with the output of a score network. The resulting model is evaluated on ImageNet64 dataset and achieves comparable results to other methods on NLL.", "questions": "1. In the combined loss (equation (8)), why you normalize the two losses with dimensionality $d$? Especially there is no $d$ in the original DSM loss.\n2. Can the proposed method be applied to conditional generative models?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n1. The proposed dual score matching objective is interesting and looks effective for energy function learning.\n2. The proposed model achieves a good result on ImageNet64 in terms of NLL. In addition, extensive experiments and analyses are provided to illustrate the comprehensive ability of the proposed method.\n3. The overall paper is well-written and organized.\n\nWeaknesses:\n1. Compared to common diffusion objective, the proposed dual score matching requires to learn two objectives. The training process might be slower. Can you provide the convergence speed comparison between the DSM-only loss and the DSM-TSM joint loss.\n2. In the experiments, it would be great to also apply the method to higher resolution images such as ImageNet $256 \\times 256$ or CelebA. In addition, can you also report FID to evaluate the sample quality?\n3. In Table 2, it would be better to add at least one more energy-based models (EBMs). And adding one figure to compare the training stability between different EBMs would further increase the reliability of the proposed model.", "was_revised": true, "final_justification": "Thank you for the detailed and thoughtful rebuttal. The authors addressed my questions clearly and provided new results and clarifications. While I acknowledge that the method is technically sound and well-presented, my main concerns remain such as limited evaluation; so I maintain my original score.", "rebuttal": "Thank you for your review. We have run several experiments to answer your questions, which will help us to improve the paper.\n\n**Computational complexity of the training process:** Training an energy model with DSM+TSM losses converges in the same number of iterations as for DSM loss alone. The additional time complexity of the TSM loss is negligible: real world training time is 4d16h for DSM and 5d for DSM-TSM, for 1M steps on an A100 GPU.\n\n**Higher-resolution images:** We have now applied our model to CelebA images at 80x80 resolution, and find that it achieves a NLL of 1.94 bits/dim. We are not aware of reported NLLs on CelebA nor other higher-resolution datasets (ImageNet64 has been a de facto standard), but it is reasonable that the value is lower than for ImageNet64 due to reduced diversity in the images. Additionally, we verified that the histogram of log probabilities (as in Figure 3) is well-approximated by a Gumbel distribution (apart from a few outliers on the low-probability edge, due to corrupted images). Over the 20k test images, the range from highest to lowest probabilities is 17dB/dim. This is half that of ImageNet64, which is also reasonable given the lower diversity in image content. Finally, we mention that scaling dual score matching to higher-resolution images is bottlenecked by the memory cost of double backpropagation. This can be alleviated with standard techniques such as sliced score matching (Song et al., 2019), which is beyond the scope of this work.\n\n**Comparison with other EBMs:** We will add a comparison with the Routing Transformer (Roy et al., 2021) that achieves a NLL of 3.43 bits/dim (compared to our value of 3.36). But we were unable to find an energy-based model (i.e., a model which directly learns the energy) that reports NLL on ImageNet64 - all models we examined were auto-regressive, flow-based, or compute variational lower bounds, while older models consider lower-resolution datasets. We also were unsure what the reviewer meant by “training stability”? Our training procedure is quite stable, as it relies on two regression objectives. In addition, these objectives have the same minimum, so that there are no adversarial instabilities (as in e.g. GANs).\n\n**Normalization of losses with dimensionality:** The factors $1/d$ and $1/d^2$ in equation (8) set the relative scales of the two objectives. This is equivalent to minimizing $t \\ell_{\\mathrm{DSM}} + \\lambda t^2 \\ell_{\\mathrm{TSM}}$ with $\\lambda = 1/d$.\n\n**Extension to conditional models:** The proposed method is straightforward to extend to conditional generative models. Given conditioning information $c$, the network computes the conditional energy function $U(y,t,c) = -\\log p(y|t,c)$. The expression of the DSM and TSM losses are unchanged, as well as the energy architecture $U_\\theta(y,t,c) = \\frac12 \\langle s_\\theta(y,t,c), y\\rangle$. We will mention this in the paper, and we thank the reviewer for the suggestion.\n\n**FID:** The generative abilities of the score $\\nabla_y U_\\theta$ derived from our energy model $U_\\theta$ are equal to that of the base score network $s_\\theta$. FID thus measures the performance of the base diffusion generative model, which is orthogonal to the focus of this work (energy-based models).\n\n**Contributions:** We would like to note that the reviewer summary omits several aspects of the paper that we believe are important contributions, regarding analyzed properties of the learned density model. We show that images with highest vs lowest probability have starkly different visual features, and that the high-probability images are not the most typical images, as intuition might suggest (Figure 3). We also show that the support of _some (but not all)_ natural images can be _locally_ described as a low-dimensional manifold. We believe these observations to be novel and of value to the NeurIPS community (as mentioned by Reviewer hMws)."}, {"review_id": "wWGUnKgcCb", "reviewer": "Reviewer_hMws", "summary": "The authors propose a framework to learn a time-dependent energy function for a distribution given a set of images. This paper presents a dual score matching framework that enables direct learning of normalized energies using spatial score matching, based on the standard denoising target, and temporal score matching, a novel target, based on the gradient with respect to the noise level.", "questions": "- Will you publish code of your experiments?\n- Could you please expand on why single score matching fails in the toy Gaussian mixture example in Figure 1? A short 1-2 sentence  explanation would be sufficient.\n-  Can you comment on generalization to larger-scale images?", "limitations": "yes", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "The submission is written very clearly, especially the introduction, which does a great job of motivating the problem and explaining the core idea. In my opinion, this paper is a very good fit for NeurIPS. One of the main strengths is the deep analysis of different aspects of the learned image distribution, where the authors go beyond just reporting likelihood scores. These insights are really interesting for the community.\n\nSmall weaknesses:\n- The example in Figure 1 needs a bit more explanation. It's not immediately clear why single score matching fails there, maybe adding a short sentence in the main text or figure caption would help.\n- In Figure 3, it would be helpful to mention directly that Figures 6 and 7 (from the supplement) show the images corresponding to the marked arrows.", "was_revised": true, "final_justification": "The authors have answered all my questions. I appreciate that the code is going to be published and the explanation of the figures. I'm raising my rating from accept to strong accept.", "rebuttal": "Thank you for the kind comments! We especially appreciate your positive evaluation of our analysis of the learned energy-based model. We answer your questions below.\n\n**Clarification of Figure 1:** Single score matching fails on multimodal distributions because the objective does not constrain the absolute energy levels of separated modes. Indeed, one can add an arbitrary constant to the energy _that is different in each mode_ without changing the score in each mode. Of course, this requires “stitching” together the different constants, which adds a term to the score, but the transition can occur in-between the modes _where there is no data_ to constrain the scores. We hope this clarifies the issue: we will update the text to be more explicit.\n\n**Clarification of Figure 3:** The purple-to-yellow-colored arrows in Figure 3 correspond to the images on the right of Figure 3 (as mentioned in the caption), not the figures in Appendix. The brown and green arrows respectively correspond to a uniform noise image in [0, 1] and a constant image of intensity 0.5 (not shown).\n\n**Code:** We will publish code to reproduce all experiments, as well as trained model checkpoints.\n\n**Scaling to higher-resolution images:** The main bottleneck in increasing image resolution is the memory cost of double backpropagation. This is shared by all energy models that are trained with a score-based objective. It can be alleviated with standard techniques such as sliced score matching (Song et al., 2019), which here corresponds to projecting the DSM residual $\\nabla_y U_\\theta - (y-x)/t$ on a random direction $v$ and minimizing its square. This enables the use of forward-mode auto-differentiation to compute the scores, as only inner products $\\langle \\nabla_y U_\\theta, v\\rangle$ are required. This can be further combined with auto-differentiation software optimizations. For instance, the ICLR blog post “How to compute Hessian-vector products?” (Dagréou et al., 2024) reports up to 3x memory reduction with the use of optimized implementations (JIT compilation, etc) in the related setting of Hessian-vector products."}, {"review_id": "fP01RPLGsP", "reviewer": "Reviewer_QWMb", "summary": "The paper introduces a method for learning normalized energy-based models by extending score matching to account for both spatial and temporal gradients. Leveraging the structure of diffusion models, the authors propose a novel dual score matching objective that enables direct estimation of energy values rather than just their gradients. The method is demonstrated on a synthetic toy example and ImageNet64, showing consistent energy estimates and competitive negative log-likelihood performance. The paper also explores the geometry of the learned energy landscape.", "questions": "1. One of the main motivations for learning normalized energy is enabling tasks like OOD detection or uncertainty estimation. Have you tried applying your model in one of these settings? Even a small experiment could help illustrate the practical value beyond likelihood evaluation.\n2. The method is only tested on ImageNet64 (and a toy example). It would be helpful to know whether it generalizes to other datasets. Even a single additional benchmark could make the results feel more robust\t.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper tackles a key limitation of score-based models - their inability to estimate normalized densities -and proposes a direct way to recover normalized energies using both space and time gradients in diffusion.\n2. Learning a normalized energy opens the door to likelihood-based evaluation, OOD detection, and deeper understanding of image statistics — things that score-based models can’t directly do.\n3. On ImageNet64, the proposed method achieves negative log-likelihood (NLL) scores that are competitive with several strong baselines, including flow-based models\n4. The paper includes a detailed empirical study of the learned energy geometry, revealing properties like variation in local dimensionality and log-probability distributions across image content.\n\nWeaknesses:\n1. The paper lacks theoretical guarantees that minimizing the proposed objective leads to an accurate approximation of the true log-density, which makes some of the empirical findings (e.g., observed geometric properties) harder to fully trust.\n2. The paper evaluates only on ImageNet64 (and on a synthetic example). Including additional datasets would help demonstrate the generality and robustness of the proposed method.\n3. While the paper focuses on likelihood and geometric analysis, it would be valuable to demonstrate downstream applications (e.g., out-of-distribution or outlier detection) where normalized energy estimates are especially useful.", "was_revised": true, "final_justification": "While some concerns remain, such as the lack of comparison to another baseline on the additional datasets, the rebuttal provides clarifications and new experiments, which helped me better understand the strengths and limitations of the work. I therefore keep the borderline accept recommendation.", "rebuttal": "Thank you for your review, which will help us to improve the paper. We ran several experiments to address your questions, which we answer below.\n\n**Other datasets:** We applied our model to CelebA images at 80x80 resolution. Our model achieves a NLL of 1.94 bits/dim. We are not aware of reported NLLs on CelebA (ImageNet64 has been a de facto standard), but it is reasonable that the value is lower than for ImageNet64 due to reduced diversity in the images. Additionally, we verified that the histogram of log probabilities (as in Figure 3) is well-approximated by a Gumbel distribution (apart from a few outliers on the low-probability edge, due to corrupted images). Over the 20k test images, the range from highest to lowest probabilities is 17dB/dim. This is half that of ImageNet64, which is also reasonable given the lower diversity in image content. We also ran our method on CIFAR-10, obtaining an NLL of 3.16 bits/dim, and a similarly-skewed log p distribution (with a range of values spanning 23 dB/dim). We hope that these additional results alleviate the reviewer’s concerns. We believe that most photographic image datasets should lead to qualitatively similar results regarding the experiments of Section 3.\n\n**Downstream applications of energy-based models:** OOD detection based on probability values is notoriously difficult, due to the counterintuitive behavior of probability densities in high-dimensional spaces. As we show in Figure 3, the highest-probability images are not typical, and could be considered outliers (this is particularly apparent in the highest-probability image, which is blank). The following paper demonstrates that generative models can assign higher probability density values to images from other datasets (e.g., MNIST) than to those in their training set (e.g., FashionMNIST):  Nalisnick, E., Matsukawa, A., Teh, Y. W., Gorur, D. & Lakshminarayanan, B. Do Deep Generative Models Know What They Don’t Know? (2019).  \n\nNonetheless, probability values can be used to detect certain degradations such as additive noise. We verified that adding noise to the reference images in Figure 3 decreases their probability (this experiment is similar to the third panel of Figure 5, except that the value of $t$ given to the model is 0). In addition, one application of our model is that it can be used to estimate the variance of the noise $t$ added to an image $y$ when it is unknown. Indeed, the value of $t$ that maximizes $p(y|t)$ provides an accurate estimate of the true value of $t$. We will incorporate figures demonstrating these two experiments in the paper. Finally, we highlight that our analysis of the learned energy model revealed novel statistical and geometric properties of image densities learned from data, which we believe to be of general interest to the NeurIPS community.\n\n**On theoretical guarantees of denoising score matching:** We briefly touch on this subtle point in Appendix A.3. Obtaining theoretical guarantees requires establishing a Poincaré inequality on the joint distribution $p(y,t)$, which Figure 1 suggests could be true, but such inequalities are notoriously challenging to prove and have only been obtained in restricted settings. However, even if this inequality was shown to hold, in practice there is no guarantee that the trained network generalizes. Proving that this is the case is an even bigger open problem. So one always needs to resort to empirical verifications, which are complicated by the fact that the “true probability density” is not known for images. This is the goal of the _strong generalization_ test of Figure 2, which shows that energy models trained on _distinct_ images compute _nearly identical_ probabilities. Our energy-based model is the first for which such a stringent test has been performed. (We also conduct a more direct verification on synthetic data in Figure 1.)  We believe these tests to give as much credibility as can be realistically obtained to the other empirical findings of Figures 3 to 5."}, {"review_id": "mxfgbXTCFg", "reviewer": "Reviewer_JFVr", "summary": "The paper has two related but distinct contributions.\n\nFirstly, the paper proposes a framework called “Dual Score Matching” which is an extension of the original Score Matching (DSM). On top of traditional DSM, the authors used the matching of “time-score”s which is a quantity similar to score but the gradient is with respect to time $t$.\n\nSecondly, the authors did not model the score directly, instead they model “normalised energy” which is almost equivalent to log-probability. They came up with a new architecture that models energy as a conservative field.\n\nIn the experiment section, denoising performance and log-likelihood comparisons have been made. Additionally, a significant part of the experiment section focuses on analysing the image distribution using the learned model.", "questions": "Major questions:\n\n- Since you are modelling the energy using a neural network but still matching the “score”, you require a “double backpropagation” (once with respect to space and then parameters). It is quite evident from Eq. 43 of appendix C.2. This is clearly a computational drawback but not discussed in the paper (except one small line in the limitation section).\n- With regard to the “normalisation” in Eq. 9, do you need to invoke the network twice for every forward pass (i.e. one for a given $t$ and another for $t_{max}$) ? If so, it is also a computational drawback. Please clarify.\n- I feel like the two parts of the paper are sort of unrelated. Is it necessary to consider energy models in order to use dual score matching? Could you have done dual score matching with score models (one model for space and another for time)? There are no ablation or justification for these questions.\n\nMinor questions:\n\n- In several places in the paper it says \"gradient of the energy is the score”, which is technically untrue. The score is the *negative of the* gradient of the energy. Isn’t that true ?\n    \n    $p(x) = e^{-U(x)}/Z \\implies \\log p(x) = -U(x) - \\log Z \\implies \\nabla_x\\log p(x) = -\\nabla_x U(x)$\n    \n- A similar issue is with the term “normalised energy (log prob)” which is also technically not true. Log-probability is the “normalised *exponentiated* negative energy”. It seems that authors used certain terms in a very loose sense.\n- Equation 2 & 3 may have a typo — shouldn't it be $\\frac{x-y}{t}$ and not $\\frac{y - x}{t}$ ?", "limitations": "Mentioned limitation. But they are rather big drawbacks of the paper and should be discussed in the main paper.", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The authors have investigated a more accurate score matching by introducing a time score component. The authors claim that this modification lead to better denoising and log-likelihoods. While improving the fundamentals of score matching is commendable, the paper has several flaws in terms of novelty, presentation, focus and evaluations.\n\n- It is quite unclear to me why the time score matching is necessary. There is not much reasoning or justification given to support usage of the time score. Is it only to improve the estimated values at the network output? The time score does not seem to have any explicit usage in the paper (e.g. in log-likelihood calculation) — am I wrong in this assessment? The only piece of reasoning I found in the paper is (L42) “[it] .. ensure consistency of the energy estimates across noise levels”. What does this even mean? There is no explanation of this sentence. Why would a score network *without* the time matching part be NOT consistent across noise levels?\n- Authors claim the second objective (i.e. the time score matching) is novel (L42; “second is novel”) but clearly there is a 2022 paper by Kristy Choi et al. (it is cited indeed) that proposes the time score and the matching objective. There is no comparison or discussion about this paper.\n- I do not quite understand how section 3 of the paper is related to dual score matching, which is the central theme of the paper. Computing entropy, log density distribution, intensity range, dimensionality — these seem to be related more to the dataset and does not have much to do with the model or the algorithm presented in the paper. I'm sure any other model that computes log-likelihood and trained properly would behave the same way.\n- I'm quite doubtful about the usage of *conservative* energy models. In the Diffusion literature, it is quite an accepted fact that conservative energy models are NOT superior to unconstrained score models. The authors accepted that their finding is (L162) “*contrary to*” the results of Salimans & Ho [2021], where they attributed this to the specific choice of architecture. However, there is no ablation regarding the choice of architecture. Also, it is unfair to compare with Salimans & Ho [2021] as they investigated generation/sampling performance but this paper did not.\n- The authors also claimed that the (L178) “.. unique advantage of our method is that .. it is fast [for NLL computation] ..”, which seems to be an over-claim. The reason it is fast has not much to do with the specific parameterisation or the matching algorithm, but the specific choice of NLL calculation algorithm that they adopted from Karczewski et a. [33]. These important details only appear in appendix rather than mentioning clearly in the main paper.\n- There are a lot of approximations made in the theoretical calculations and justifications. On one hand you argue that your method is better due to “homogeneity” property but on the other hand you say that the network is *only approximately* (L508 appendix) homogeneous. On one hand, you argue about your energy functions being “normalised”, only to reveal later that it is only an approximation. Also, the explanation in appendix A.3 is quite hand-wavy, non-rigourous and does not conclude anything concrete. In general, it seems that the authors are quite uncertain about their rationals/claims.\n\nOverall, I am do not particularly see much novelty (or proper justifications for them) convoluted by the fact that the many rationals are quite uncertain and not concrete. Also, the paper lacks proper experimentation to validate the claims.", "was_revised": true, "final_justification": "The rebuttal clarified some core concerns. But the paper needs to be polished quite a bit. There are still downsides about the paper, regarding practicality and overall impact. I am increasing the score one point, but still slightly lean on the negative side.\n\nSee my rebuttal response for details: https://openreview.net/forum?id=wtYcS4kxpF&noteId=xYq2B0yBeF", "rebuttal": "Thank you for the detailed review, which will help us significantly improve the paper.\n\n**Clarifying contributions:** We believe that the primary contributions of the paper are twofold:\nA combination of architecture and objective that allow efficient training of an accurate normalized energy model. We model the energy directly (rather than the score) so that energy computation is done in a single forward pass. The purpose of the time score objective is to ensure that the energy values are accurate and normalized (see below). They go hand in hand: in our setting, time score matching is meaningless for score models (see below).\nAnalysis of the model (trained on ImageNet) to elucidate properties of natural image densities. We believe these observations to be novel and of value to the community (as mentioned by Reviewer hMws).\nThese two contributions are not disconnected. While most of the experiments in 2. could in principle be done with any generative model that can compute probabilities, they are much more easily done when the log probability is explicitly available with a single forward pass. Furthermore, the dimensionality measurements require access to the distribution of noisy images. Our energy model is the only one that provides access to the density of clean or noisy images in a single forward pass. We will incorporate these points into the introduction.\n\n**Novelty:** Beyond our combination of energy architecture and dual score matching objective which is novel (see below for time score matching), we are not aware that any of the results in Section 3 have appeared in the literature. Could the reviewer provide references that establish similar results?\n\n**Purpose of time score matching:** While the time score is not used at test time to evaluate the energy, it plays a crucial role during training to ensure that the learned _energy_ (rather than the “space” score) is accurate. This is demonstrated in Figure 1: the energy model trained with single score matching does not learn the ground truth energy, while dual score matching succeeds. The reason is that single score matching does not sufficiently constrain the values of the energy when the data is multimodal: one can add _different_ constants to the energy in disconnected modes without affecting the score. Matching time scores fixes these unwanted degrees of freedom. This is because the _joint distribution_ $p(y,t)$ is unimodal, contrary to $p(x)$, as visualized in the right panel of Figure 1. Matching its score, which contains both the “space” and time components, thus constrains the energy up to a _single_ additive constant, as opposed to one per mode of $p(x)$. We will clarify this central point in the paper, which was not sufficiently explained.\n\n**Novelty of time score matching:** While time score matching has been proposed before, our  particular objective is novel. In Choi et al., the time score matching objective uses a second derivative in time, which is more computationally expensive than our regression-based objective. The conditional expectation identity (4), and the resulting objective (5), are novel compared to this prior work. However, after submission of the paper we became aware of the concurrent work “Density Ratio Estimation with Conditional Probability Paths” (Hanlin et al., 2025) that derives the same identity and objective. We will incorporate citation and discussion of this into our paper.\n\n**NLL calculation:** The reviewer is mistaken. The fact that our NLL calculation is fast (a single forward pass of the model $U_\\theta$) is _specifically_ because of our choice parameterization and objective (this was the motivation for these choices). Our energy model can compute the probability of 50k images in 12s on an A100 GPU, whereas a diffusion model requires upwards of 3h20 as it requires at least 1000 forward (+ backward) passes per image (assuming 100 time steps and 10 noise samples for computing the divergence with Hutchinson trick or the MSE as explained in Appendix A.2). Appendix A.2 provides a literature review of other methods of estimating energies, **not** a description of our method. We do not rely on any algorithm from Karczweski et al. (who do not consider EBMs).\n\n**Normalization:** Normalization is computed only once post-training, by calculating the constant $c = \\mathbb E_y [U_\\theta(y,t) | t=t_{\\max}] + \\frac{d}{2}\\log(2\\pi e t_{\\max})$. Test-time normalized energies are then obtained with a single forward pass by calculating $U_\\theta(y,t) - c$.\n\n**Conservative diffusion models:** The goal of the paper is not to obtain better diffusion models by imposing them to be conservative. Rather, our primary goal is to learn the energy, a byproduct of which is a conservative score model. We merely point out that due to our choice of architecture, there are virtually no differences between the original score network and the gradient of the learned energy.\n\n**Computational drawback:** We do not understand the question. We mentioned the limitation of the added computational cost of the double backpropagation in the main text. To be more precise, the training time for a score model for 1M steps on ImageNet with an A100 GPU was 1d8h, while that of our energy model was 5d. We will mention this in the paper and apologize for the imprecise statement.\n\n**Approximations:** As in any empirical work, there are approximations involved. We will update the text to make these more explicit. It is common practice in energy-based models to resort to approximate normalization as exact normalization is typically impossible (continuous integrals can only be calculated by discretization unless a closed form is available). The purpose of Appendix A.3 is not to prove theoretical guarantees (we say this explicitly) but to provide a sketch of what would be needed to establish them for future research.\n\n**Terminology:** We apologize for the confusion. As an abuse of language, we called the gradient of the energy the score rather than negative score. Thus, there is no typo in eqs (2) and (3): $U$ is the _negative_ log probability, not the log probability (see equation (1)). However, we disagree with the second remark on terminology. Log probability is equal to the negative energy up to an additive constant $\\log p(x) = -U(x) -\\log Z$, not the normalized _exponentiated_ negative energy. We used “normalized energy” to refer to $-\\log p(x)$, a term that for instance is used in the review by Song and Kingma (2021)."}]

Author rebuttal:

Meta-review: This work proposes a new method for learning normalized energy-based models. One of the weaknesses of this work, raised by reviewers, is the much longer training time that might prevent its practical use. Authors discussed this limitation in the main paper and in the rebuttal and suggest it could be useful due to its faster inference. Overall, the AC agrees with the merit of this work for the various analyses and insights provided by the proposed approach in the paper. My recommendation is acceptance. Please incorporate the feedback from the reviewers on clarity and the terminology used in the final copy.
Decision: Accept (poster)