Score-Based Diffusion Modeling for Nonparametric Empirical Bayes in Heteroscedastic Gaussian Mixtures
Authors: Gongyu Chen
Ying Cui
Venue: NeurIPS 2025 poster ()
Keywords: empirical Bayes
Gaussian location mixture model
g-modeling
heteroscedastic noise
diffusion processes
score matching
Primary area: probabilistic_methods
TL;DR: We propose score-based diffusion framework with explicit G-modeling, a nonparametric empirical Bayes method that achieves near-parametric score estimation guarantees and state-of-the-art denoising for multivariate heteroscedastic Gaussian mixtures.

Abstract: We propose a generalized score-based diffusion framework for learning multivariate Gaussian mixture models with homoscedastic or heteroscedastic noise. Our goal is to nonparametrically estimate the latent location distribution and denoise the observations. 
Departing from the conventional maximum likelihood approach, we reinterpret each observation as a temporal slice of a family of stochastic diffusion processes. This modeling choice enables a principled characterization of the additive noise structure and supports a multi-step denoising procedure grounded in reverse-time dynamics. We introduce a score-based objective that explicitly models the latent distribution and accommodates observation-specific noise covariances.
Theoretically, we establish that the score estimation error with $n$ independent observations achieves a near-parametric error rate of $\frac{\mathrm{polylog}(n)}{n}$, 
improving upon existing results in the diffusion literature. Empirically, our method outperforms the nonparametric maximum likelihood estimator in both density estimation and denoising fidelity,  especially in high-dimensional settings.
These findings suggest a promising direction for integrating nonparametric empirical Bayes with diffusion-based generative modeling for latent structure recovery.

Score-Based Diffusion Modeling for Nonparametric Empirical Bayes in Heteroscedastic Gaussian Mixtures
Gongyu Chen Department of IEOR UC Berkeley Berkeley, CA 94720 gongyu@berkeley. edu
Ying Cui Department of IEOR UC Berkeley Berkeley, CA 94720 Singcui@berkeley. edu
Abstract
We propose a generalized score-based diffusion framework for learning multivariate Gaussian mixture models with homoscedastic or heteroscedastic noise. Our goal is nonparametrically estimate the latent location distribution and denoise the observations. Departing from the conventional maximum likelihood approach, we reinterpret each observation as a temporal slice of family of stochastic diffusion processes. This modeling choice enables a principled characterization of the additive noise structure and supports a multi-step denoising procedure grounded reverse-time dynamics. We introduce score-based objective that explicitly models the latent distribution and accommodates observation-speciie Theoretically, we establish that the score estimation error with independent observations achieves near-par pear-parametric error rate of polylog(m)/n, improving upon existing results in the diffusion literature. Empirically, our method outperforms the nonparametric maximum likelihood estimator both density estimation and denoising fidelity, especially in nigh-dimensional settings. These findings suggest a promising direction for integrating nonparametric empirical Bayes with diffusion- based generative modeling for latent structure recovery.
noise covariances.
Introduction
In many large-scale scientific domains, such as the analysis of high-dimensional astronomical data, measurements are often corrupted by heteroscedastic noise, where error variances vary across dimensions and differ across observations [Kelly] 2012 Akritas and Bershady 1996 Anderson et al.) 2018). While this noise complicates downstream inference, is frequently well-chhracterizzd through the data collection process. A central goal in such settings is to recover the true latent signals and their underlying distribution from noisy observations by leveraging structural information inherent in the data. This problem can be naturally modeled using heteroscedastic Gaussian location mixtures, described by:
X; 0 +2
with Lid G* and Z ind (O, E.),
for i= L,....
where (Xi) C Rd are the observed data with distinct error distributions, (0) ] the unob- served latent signals drawn from a common unknown prior measure G" and (2:) are independent Gaussian noise with known positive-definite covariance matrices )
and
Model ( falls within the empirical Bayes (EB) framework [Robbins] 1951 1983), where the goal to estimate the prior distribution G and recover the latent vectors (0) ドー One way to achieve the goal is through the G-modeling approach [Efron] 2014], which is directly modeling the prior on the parameter scale and then plugged into the posterior calculation. A principled and widely used approach to estimate G*, without imposing any parametric assumption on it, is the nonparametric
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

maximum likelihood estimator (NPMLE) Robbins 1950 Kiefer and Wolfowitz 1956 Tiang and
Zhang 2009 Gu and Koenker 2016 Polyanskiy and Wu 2020]. The NPMLE seeks a distribution G E P
Rd) that maximizes the marginal likelihood of the observed data {X
GMLE E arg max GEP(3") n log f,,,,
where P(R4) is the family of all probability measures on Rd and fG (x) := 4s (z 0)dG(O)
is the marginal density of X, under G with yr, (v) (2et(2 exp (- 3y x y). Once
G is obtained, each latent signal 0g can be estimated using the empirical Bayes posterior mean:
0, = Bgl0 X, E
X; Σ log fā (X;),
where (*) is due to the Tweedie's formula Efron 2011 linking posterior means to V log fā, (X.). the score function of the marginal density. This nonparametric approach flexibly adapts to latent structures without strong modeling assumptions on G Theoretical guarantees for the NPMLE in this model well established in the homoscedastic setting Σ, E the univariate setting d 1: it achieves rear-parametric sample complexity rate Onn (log n) (1)d)), i.e., parametric up to logarithmic multiplicative factors, for both density estimation and Bayes regret risk [Jiang and Zhang 2009 Saha and Guntuboyinal 2020]. Recently, the results have further been extended to multivariate
and heteroscedastic setting in Soloff et al. 20255].
Despite these guarantees, the NPMLE and its associated prredie-baseed denoising procedure suffer from inherent statistical and computational limitations that render their practical effectiveness in high-dimensional settings: (a) Statistically, Tweedie's posterior mean estimator often has an
X
over-shrinkage bias. As an example, consider a joint Gaussian model for (X, ("):
0"
21a la
1)- It follows that marginally l" ~ NN, I2), and X00 ~ (N08", la la
However, the posterior mean of 0* is E| X/2 (0, Ia/2), thus underestimating the variance of the latent signal. (b) Computationally, although problem (2) is convex optimization in G, it infiniteddimensional and difficult to solve in practice. Since it known that there exists solution GALE with support on at most n atoms, one common way to address the issue is to restrict G' to be discrete measure of the form Ei= wey with =i=1 W; = 1 for some S Optimizing over {(w, m1)) C R4 Rd leads to finite dimensional but nonconvex optimization, often tackled via the inspecttiio-mmaxmmzzatio slow convergence and may converge to local maxima due to the nonconvexity. To mitigate these issues, prior work has suggested fixing the atoms { / using either grid points over a compact region [Koenker and Mizera 2014 Feng and Dicker 2018 or the observed data itself (Lashkari and Golland 2007]. Yet the grid-based based approach scales poorly with dimension, and the data-driven strategy is suboptimal under heteroscedastic noise [Soloff et al. 2025). Figure further illustrates the over-shrinkage and computational challenges that arise from these limitations.
(EM) algorithm [Laird] 1978). However, EM suffers from
Jobby Duta
Figure 000 twn-dimensional observations X; (in blue) sampled from X: A/ (0.3) uniformly distributed unit circle (in black), The posterior
red), obtained via the NPMLE and [ over-shrunk relative the unit circle. Right (udapted from al. 2025 Figure 7]): a real- world 19-dimensional astronomy data (in black). The estimated posterior means (in blue) from NPMLE poorly reveal the underlying chemical structure Section for the denoising results from methods).

---

1.1 Our Methods and Contributions
To address the statistical and computational limitations of the classical NPMLE and Tweedie-based posterior denoising, we propose a novel perspective that reinterprets the data generation mechanism in model (I) through the lens of stochastic diffusion processes. Our method built on two core ideas: replacing the NPMLE with score matching to estimate the prior within the G-modeling approach, and applying multi step posterior denoising to recover latent signals. Below, we present an overview of our approach in the simplified homoscedastic setting; the full treatment of the heteroscedastic case developed detail in Section 2
Modeling the data measurement process. Consider the homoscedastic version of model (T. where Σ, Σ for all 1,... n. In this special case, the observed data { in are i.idd samples from the convolution of the latent distribution G* with the Gaussian kernel (0, Σ):
(4) Rather than treating each observed data X; as the result of static additive noise model, we propose interpret {= as snapshots from a forward stochastic diffusion process. Specifically, we model the measurement process as zero-drift stochastic differential equation (SDE) that continuously perturbs samples from the latent prior G
X; ~ G" * (N(,E),
for i = 1,..., 71
dx = g(t Σ dd with (() {0;1:=1 ~ G",
(5)
where B, Rdd denotes standard Brownian motion and g R2+ R4 is the diffusion coefficient
controlling the noise scaling. The function g(t) can be adapted to reflect assumptions on the measure-
ment process; for simplicity, one natural choice is g(t) V=nee = (tt) = (Jo g(8) =
Under this model, the marginal distribution of (+) evolves according to:
222 e Σ. == at)
x G N (0,E.)).
where I0)
In particular, at time = 1, we have E(1 = E. Therefore, we can indeed view each X; as sampled from the diffusion process (5) at 1, i... (X1 x(1..
Denoising the observations. The forward-time process in (5) is associated with a reverse-time diffusion process, which can be approximated by a deterministic flow that shares the same marginal distributions over time [Anderson, 1982 Maoutsa et al.) 2020 [Song et al., 2021a]. This flow is
governed by the following ordinary differential equation (ODE):
dx - -t)) V. log for yo (X ")) dt.
(6)
Given estimator G of the prior G*, this reverse time ODE naturally defines generalized multi- step Tweedie-type empirical Bayes denoising procedure. Specifically, replacing G* with & and discretizing the ODE backward in time from f = 1 t = 0 with step size At yields the iterate:
0 g(0) + 2
At g())² E Va log fa (0()) with initial condition 0:10 = X;-
The final denoised estimate is taken as 0, = 0,00 i n. Notably, this method generalizes Tweedie's 's formula in (3) which corresponds to a single-step approximation with At 1. Our multi-step denoising procedure mirrors optimization dynamics, particularly Langevin-type flows over probability spaces, where smaller step sizes yield more controlled updates. This approach allows gradual refinement of the denoised signal, helping mitigate the over-shrinkage bias observed in one-step Tweedie denoising methods.
Learning the latent prior G*. Our stochastic framework also naturally induces denoising score matching (DSM) objective [Vincent 2011] for learning the prior G". We propose to further perturb the noisy observations X; using the forward SDE (5) to generate conditional samples from the kernel
Pie( (X,) = fg _ - -X,), up to some stopping time T. We then match the score of the model's marginal density to the conditional score induced by this perturbation:
1 1 Emi(ix) [a(t)² | log fa (z) - Vlog dt, (8)
G = arg min
GEP( 1=1
3

---

where | := x: I x' denotes weighted fz-norm, extending the likelihood-weighting idea of Song et al. [2021b]. Intuitively, under our stochastic view, the observations at 1 already contain calibrated noise relative to G", and additional perturbation simulates conditional transitions. Matching model scores to these local conditional scores allows global recovery of the latent prior. Unlike standard diffusion models with the DSM objective that directly parameterize the score function with neural networks [Song et al. 2021a], we model the distribution G' explicitly in ( Such G-modeling is crucial in deriving new theoretical guarantees for the score estimation risk. Computationally. although the objective D is nonconvex, diffusion-type squared loss objectives have been empirically observed to exhibit smoother optimization landscapes and avoid poor local minima [Xu et al., 2024] This empirical tractability is further supported by theory on the benign optimization behavior of overparrameterized neural networks trained over meen-squaree losses that can typically satisfy the condition or behave linearly in the neural tangent kernel regime [Liu et al.
2022].
In summary, we propose a novel integration of score-based diffusion models with the classical EB framework via G-modeling. The rest of the paper is organized as follows: In Section 2 we present a stochastic process framework to model heteroscedastic Gaussian mixtures, where each observation is treated marginal realization from distinct diffusion process with its own score function. Section B eestablishes the theoretical advantages of our approach, showing that G-modeling combined with diffusion yields near- peaa-parametric rates for score estimation. We further discuss how to numerically solve the G-modeling based score matching objective in Section [ In Section [5] we present experimental results on both synthetic and real astronomical data, demonstrating substantial improvements in denoising and prior recovery over NPMLE and one-step Tweedie- based methods, supporting the improved computational scalability of our framework in high-dimensional settings.
1.2 Related Work
Score-based models for mixtures. Recent theory has provided provable guarantces for diffusion and score-based models in learning Gaussian mixtures, yet these results neither handle the heteroscedastic measurement setting we study nor aim at recovering the latent prior itself. Works such as Gatmiry et al.) [2224] Chen et al. [2024], [Shh et al. [2023) propose diffusion- -based algorithms for mixture models under either finite discrete structures or continuous latent mixtures with homoscedastic noise. These approaches typically assume access to an infinite supply of i.i.d. samples and yield slow convergence
guarantees of order ((logn) -1/p) in score or total variation distance, for some exponent p > 1. Importantly, these works directly parametrize the score function, using either neural networks or piecewise polynomial surrogates, without attempting nonparametric recovery of the latent mixing measure G or incorporate explicit G-modeling into the objective. Because of this direct score parametrization, these methods cannot be readily applied to our core task of latent recovery and are incompatible with our de- de-bi multi-step denoising procedure. A notable exception is Ghosh et al.. [202]]. which presents the first analysis of combining score matching with G- moddelngg achieving a near parametric rate Olog n/n) for score estimation risk. However, their results are limited to the univariate, homoscedastic setting (d = 1). In contrast, our work tackles the more challenging regime of multivariate, heteroscedastic Gaussian location mixtures with arbitrary dimension d 2 1. By coupling G-modeling with a stochastic process perspective, we develop an explicit estimator of the latent prior G* and establish near-pane sear-parametri rates for score estimation (up to small smoothing). along with imite-sample guarantees for recovering the latent signals.
Image denoising in computer vision. While our problem and the proposed method may appear superficially related score-based diffusion approaches for image denoising, the settings and statistical challenges are fundamentally different. Most notably, our task lacks access to clean ground truth data, precluding the use of pre-trained score models, and all inference must be performed from noisy observations alone. In contrast, supervised denoising approaches such as BM3D [Dabov al. 2006) requires access to large number of paired noisy-clean images, while self-supervised methods like Noise2Score [Kmm and Ye 2021) rely on one-step Tweedie approximations that ignore the unknown prior G and leverage only partial model information Σ, These methods are either unrealistic or statistically inefficient in scientific domains where latent structure is present but clean data are unavailable, such as the astronomy measurement data. Even moderate dimensions (e.g., d = 20), empirical Bayes is a highly nontrivial statistical task due to the data sparsity: each observation X, corresponds to a unique latent parameter 0ž and an individual covariance E,, yielding just one single noisy sample to infer each 0".
4

---

2 SDGM: Score-based Diffusion with G-Modeling
In this section, we present our general treatment of the score-based diffusion approach with G- modeling under the heteroscedastic setting, The key idea is still to reinterpret each noisy observation as time-1 snapshot from latent stochastic diffusion process. While our homoscedastie formulation serves as motivating special case, the heteroscedastic setting introduces a crucial distinction: each observation X; corresponds to a different SDE, despite all sharing common latent initialization. Specifically,, we consider the following family of zero-drift linear SDEs:
dx -ttt) dBe with x; G", for - 1..,..,
where {9:(t) C R are the time- denendeenddntt diffusion coefficients. Without loss of generality, we assume ai( = S0 gi (s) satisfies a; (;( 1 for all i. Then the time-dependent covariances are given by 21(0) a;; (t) Ei, and the marginal law of X becomes
x!! ~ G * *N E for t> 1=.,....... n.
(10)
Thus, each observation X; in the original mixture model ( can be viewed as a marginal sample from this diffusion process at i..., X; _ x! G (0, E,) for all i 1,..., As in the homoscedastic case, each forward SDE in p is coupled with the reverse-time ODE:
xx = 90 E, V log x11) 2 (X,) dt,
for I....,n. n.
Given any plug-in estimator 7 for G substituting E(0 x01) in defines multi-step score-based denoising procedure that remains consistent with the EB framework. Furthermore, as A/ - 0, the dynamics in 77 induces deterministic transport map from the noisy empirical measure toward
latent measure G = G This perspective aligns with the design principles behind recent optimal transport approaches for empirical Bayes denoising (Soloff et 2025 Zhang et al. 2024 mitigating the over-shrinkage limitation of one-step Tweedie-type procedures.
Importantly, in our setting, the SDEs in serve as belief model for the true data bata-generating process, and the noisy observations in T are treated as marginal realizations from latent diffusion paths. This fundamentally differs from the generative score-based diffusion models like (Song et al. 202Ia], where SDEs are artificially constructed for the purpose of noise injection and score network training. In an idealized setting where the true latent measure G is known accessible via unlimited samples, one could apply the framework of Song et al.) [021b to learn a time-dependent score network so(- t), mapping G* pure noise through an invertible path. However, setting departs significantly from this setup. First, we observe only a finite collection of noisy samples at forward time / 1, each generated by unique, heteroscedastic SDE rather than homogeneous, artificially perturbed noise. Second, our goal is to explicitly estimate the latent distribution G*, as opposed to modeling it implicitly via score functions. Finally, because no data available at intermediate times t e (0, 1), the standard inne-indexed score network sa(-, t) cannot be trained in our setting.
To address these challenges, propose new objective that aggregates information across all diffusion paths to estimate G Our method, which we term score-based diffusion with G-modeling (SDGM), obtains the estimator G by:
G E arg min (G) GEP(RA)
/' (.(t)
fi V log (e V log g (x)
(11) dt, E.
where, to 2 1 is the minimum diffusion time, and to simulate the homoscedastic behavior, we consider the rescaled Gaussian transition kernel: i, ~ Ke(- := (((aa(t)),) (1
(a() defined fort 2 to 2 1, and define the mixture densities (4) and de to be:
n
(empirical) ( == n Σ K.( (i).
(G-induced)
de
n Σ E (E!)) 1/2
where T#G denotes the pushforward measure of G under map such that for every set C Rd, (t) (A) G(T-IA), A)). The defined mixtures satisfy that Equ. 4g. Intuitively, 46 aggregates
(0)
5

---

the global behavior of the latent measure G across all SDEs, while dh serves as its empirical counterpart based on noisy observations. By matching the score functions of these two mixtures, we guide G toward an accurate approximation of G In the homoscedastic special case with Σ; = E and gi - 9, the mixtures simplify rescaled version of the objective (IT) that is equivalent the score matching formulation in (8) up to constants.
Theoretical Guarantees
In this section, we provide guarantees on the statistical performance for both homoscedastic and heteroscedastic noise models. We primarily study three statistical properties to quantity how well the SDGM estimator G of (IT) performs as a plug-in estimator for the true prior G*.
The successful convergence of diffusion- baased generative models is underlaid by low score esti- mation error (Yakovlev and Puchkin, 2025 [Bortoli] 2022]. Define the following time-cumulative average Fisher divergence for all path measures μ and
B (ulle) = Σ 1=1
B 2
50
V log (1) V log e0 122d.
Then, the specialized metric 3,c, (Ja- p200) ((ā measures the time-cumulative risk on how gge along the forward-time SDE path beyond the observable time. It is consistent with the generalization risk concerned in Yakovlev and Puchkin [2025]
V log go approximates V log fc. gee
For density estimation, G also provides a natural estimate for the mixture marginal law IG- aggregated from the paths of (9). We measure the quality of density estimation by the standard squared Hellinger accuracy, 25 (p, := [(Vp / between densities and q. In particular, we are interested in the following squared Hellinger distances:
(t)
(homoscedastic) 5y² (fa 200 fa (heteroscedastic)
(1) (t)), (96- .8
which captures the global of the plug-in estimator G for estimating the joint law of {X (x(") 1121-
The third metric of interest we establish is the squared Wasserstein distance (see e.g.. Villani (2008)) ). which is defined for all μ., V E P(R") that W (14, v) == inf B(z,y)>>1llaa y v with [1,, v)
denotes the set of all joint measures over μ4, v. The distance W G*) naturally characterizes the direct consistency of G for the true G
We make the following assumptions on the latent structure and the data generation process (9) In the assumptions, the support conditions are widely considered in literature [Jiang and Zhang 2009 Saha and Guntuboyina 2020 Ghosh al., 2025). They encompass a wide range of realistic scenarios, including discrete mixtures, compactly supported manifold -uspooed- latent structures, and representative experiment settings in Section [ Furthermore, well-conditioned assumptions on noise covariances are ubiquitous in prior works [Soloff et al. 2025) and are essential to ensure smoothly- defined score everywhere and to enable information-thhorrtiaall meaningful latent recovery.
Assumption 1. The true latent measure in model ① is compactly supported, i.e., G* E P(--M, M]d) for some finite radius M > 0. There exists < a S σ < oo such that the co- variances satisfy gla = E; 3 Fla all i == .n. The diffusion coefficient satisfies a;(
(such that (2)) = Σx) and other mild regularity conditions specified in the supplemental.
Our first theorem characterizes the performance of our SDGM framework for the canonical ho- moscedastic model 44) which is specialized from our guarantee for the heteroscedastic model.
Theorem (Homoscedastic noise). Under Assumption D aad further assume that Σ; = E, 9:
for all i = =,... Let G minimizes objective (8) constrained over the family P[[-M. (|]) and with parameters to and T satisfying large, the following bounds hold with probability at least - n -2
log n S a(to) S ((T)) < 2. Then provided n sufficiently

---

Score and density estimation:
2 SM.d - (log n 2d+3
(12)
(eoe,,, SMddaa) (log n) 22++3 )
(13)
Deconvolution:
(G 1 logn
(14)
Here, we write f(n) S g(n) if there exists universal constant C independent of n such that ['m) I Cg(n), Vn. We write S0,.0p O0.... to emphasize that the constant C' Cp also depends on the fixed specification parameters 0,'s. In the following, we discuss implications of this result. Under a latent manifold setting with homoscedastic isotropic noise and neural parameterization of the score function, recent work Yakovlev and Puchkin 2025 proves the rate ET N7. (fa- yo (lfr, =
On - pn)) (ollyooomm where do the intrinsic dimension of the manifold and / denotes a smoothness parameter. To our knowledge, this is the sharpest rate known for diffusion-bar -baed score estimation under such assumptions. By leveraging explicit G- G-modeling to impose the structure on the estimated scores, result achieves improved O( leading dependence n (up to logarithmic factors), under more general assumptions. Moreover, our Hellinger accuracy bound (13), modulo a small distributional shift from / = 1 lo, matches the rate achieved in NPMI.K .E literature Saha and Guntuboyina, 2020]_ If the target is to estimate fa- xx at the observational noise level, an unavoidable smoothing cost due to diffusion arises. In this regime, the task incurs the inevitable minimax rate On (2((44))) that cannot compete with NPMLE [Wibisono et al. 2024). Finally, in the Gaussian deconvolution setting with fixed noise design, the minimax rate deteriorates to the logarithmic regime O log n) due inherent ill posedness induced by Gaussian smoothing and it is the best rate possibly achieved by both NPMLE and our framework [Soloff et al. 2025 Dedecker and Michel 2013). We further discuss limitations of the current theory in Section
sample size
Our next theorem generalizes the results to the heteroscedastic settings. The core contents of Theorem 2 parallels that of Theorem [ but additional care is needed to account for observation- specific noise covariances and SDEs. To make the readers easy to follow, we present the main result informally below and defer detailed technical derivations and proofs of Theorem and 2 to the supplemental material.
Theorem 2 (Heteroscedastic noise; Informal). Under Assumption 7 provided solution to constrained over the family P([-M, M)) with (lo, T) satisfying well-behaved smoothing condition, then with high probability.
(log g) 20+3
(log n) 2d-3
(to) (96-
110
(2a)
x
(ta) E5y
(50)
EW?
1
[E4)
112
4/2
log
4 Solution Methods to Solve 11))
The optimization problem in (II) is an infinite dimensional nonconvex problem that is challenging to solve numerically. We discuss practical strategies for approximately computing the SDGM solution.
Strategy one: discretize the support of G. A classical approach to deal with the infinite dimensional variable G, as commonly used NPMLE, is to restrict the class of priors to discrete measures with finite number of supports: G = Wyywy, with Σ;= Wy = 1, wy 2 0. This reduces the optimization to dinee--imessional, albeit nonconvex, problem over {(wj, k,))
Strategy two: continuous modeling with normalizing flows. An alternative approach is to model G continuously using normalizing flows [Rezende and Mohamed 2015 Dinh et al. 2017 Durkan
jt al. 2019]. normalizing flow represents a flexible distribution as the pushforward of a simple base distribution (c.g., standard Gaussian) through invertible neural network. We parameterize G
7

---

using a Neural Spline Flow [Durkan et al. 2019 implemented via the nflows library (Durkan et al. 2020]. allowing efficient sampling and density evaluation. In our framework, we draw Monte Carlo samples from the flow to estimate fa and its score during training. Sampling from a continuous prior theoretically justified by empirical Bayes approximation results: given a continuous G, there exists discrete approximation G' that closely matches fG and its gradient [Saha and Guntuboyina, 2020 Lemma D.3]. Furthermore, continuous modeling enables richer recovery of the underlying structure of G, which is crucial for scientific interpretation and denoising applications [Shen and Wui 2022
Batched and conditional estimation of qG- The full mixture (G plays central role in our theoretical guarantees, but evaluating its score Vz log 9G is computationally intensive for large n. A natural relaxation is to approximate qG at each iteration using only a mini-batch of obser- vations. In particular, we focus on the case where the mini-batch size B = 1: at each training step, we align the rescaled conditional score log f cx (x) with the empirical conditional score
v, log. (ar; X;, (t)(()) 1) This corresponds to a conditional score matching update, mirroring the classical intuition of denoising score matching (Vincent, 2011], while offering computational sim- plicity and stability in stochastic optimization. In fact, subsampled estimation for mixture inference has also been used in early works [Celeux et al. 2001 Thler et al. 2003
5 Numerical Experiments
In this section, we evaluate the empirical performance of our proposed method synthetic het- eroscedastic datasets and real-world astronomy data. We also investigate whether the empirical scaling behavior aligns with our theoretical guarantees, particularly refer the readers to the supplemental for more implementation details and results.
vary the problem dimension. We
Baselines and metrics. We evaluate two variants of our proposed method: SDGM-C, which uses a continuous prior modeled via normalizing flows, and SDGM-D,, which uses discrete prior. We compare our approaches against the ofte--of-thearr NPMLE solvers for the problem implemented in literature. Specifically, we compare with (1) NPMLE-N, the discretized NPMLE with supports taken fixed at noisy data (thus the objective (2) is optimized over only weights and becomes convex),
as advocated in Lashkari and Golland| [2007] |; we adopt the Newton- based augmented Lagrangian method implemented from [Zhang al. [024].
(2) PEM, the partial EM method which is generalized EM method proposed in Zhang et al. 2022 that used tighter sequential relaxations.
In addition, we consider hybrid heuristic called SDGM-P to further improve the likelihood value defined in (2). This method first trains the SDGM-C model to learn a continuous prior, from which we sample support points. These support points are then used for a post-processing step via the convex NPMLE-N solver, effectively refining the prior estimate to maximize the likelihood.
We evaluate the performance of different methods over following metrics: (1) Fisher divergence
(FD) defined by (f. f.) := " |ol far.a--... fa-,,, (X,)
log fā,x, (X.) (2) Negative
log- ligliikl (NLL) defined by
2i=1 log fā.z, (X+), which serves as a proximal for the KL
divergence between fā.5, and fa P (3) Squared Wasserstein-2 distance (W3) between estimated G and G*, approximated by the entropy-regularee
optimal transport distance [Cuturi] 2013 Genevay
al. 2018 (4) Mean squared error (MSE) between the denoised and clean data EE 10, 0;
Synthetic datasets. We evaluate our method on two synthetic datasets with oiw-dimensional manifold structure, both of which are widely used in the NPMLE literature (Saha and Guntuboyina 2020
Zhang et al. 2024]. Specifically, for d 2 2, we generate clean latent vectors 0: that exhibit nontrivial structure only in the first two dimensions, with the remaining dimensions padded with zeros:
Example (cirele): The first two coordinates of each 0 are drawn uniformly from a circle of radius 6 centered at the origin.
Example 2 (uppercase letter "A"):: The first two coordinates of each 0 are drawn uniformly from one of five line segments that form the shape of the capital letter "A." These segments connect the points (-4,6), (-2, 0). (0, 6). (2, and (4, 6), with each segment selected with equal probability 1/5.
We then corrupt each 0, with independently generated heteroscedastic Gaussian noise Z, (O, E.). where each Σ, is a diagonal matrix. The diagonal entries of E.; are sampled independently and uniformly from the interval (1.2)).
8

---

10", The indicales that lower values are betier. We report Table 1: Models irained with dimension d 16 and sample Denoising of the baseline methods NPML F-N and PEM performed vin generilization risks evalunted on a separate test dataset size 10*, the one-stop Tweedic's formula, while our methods the solver. Rosults stochastic methods averaged over independent runs.
circle
uppercase letter "A"
metric
FD
NLL
W
MSE
FD 1
NLL
W,
MSE
NPMLEEN
1443
3.57
25.10
3.04
PEM
4,71
24.89
4,52
1.26
3.78
24.40
4.09
1.35
SDGM-C
1.27#0.04 04
21.82+0.01
2.07
1.13+0.01
1.23++0.11 11
21.33++.001
1.68+0.31
1.32=0.08
SDGM-D
3.7220
24.88
3.9610
1.04:0.03
6.2611.05
24.41
3.94
1.73
SDGM-P
0.68g0.06
24.78.0.00
1.07::.00
4.20
0.49.0.02
24.30=0.00
1.25g0.07
1.45g0.02
The results are summarized in Table and Figures and where G of our methods are plotted using the marginal density estimated through geodesic- KDE Viincent and Bengio [2002]. As shown in T SDGM consistently outperforms the NPMLE baselines across metrics, with the additional post-processing step (SDGM-P) further improving the quality of prior estimation in most cases. The
denoising effects of different methods are evident Figures 2 and 3 Specifically, the estimated
signals 0 from NPMLE and one-step Tweedie denoising appear overly blurred and shrunk, failing to recover the underlying manifold structure. In contrast, SDGM, equipped with multi-step ODE-based denoising, is able to clearly capture and preserve the latent geometry of the data.
Real-world astronomy data. We further demonstrate SDGM-C on challenging world scientific problem: estimating the latent distribution of stellar chemical abundances and denoising observations from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) Data Release 14 Abbofaahh et al. 2018 Majewski et al,, 2017 With pre-processing in Ratclifte et al. 2020). this dataset, comprising (e.g., [Fe/H], [Mg/Fe], [Si/Fe]), features heteroscedastic,
27, 238 stars with d 19 chemical abundance measurements
correlttd Gaussian noise for each star. The true underlying abundances Or are expected reside on
manifold
governed by astrophysical nucleosynthetic processes (Ting et al., 2012 Weinberg et al. 2019)
Accurately revealing the underlying physical structures through the denoised abundances crucial for reliable studies of Milky Way formation and evolution [Freeman and Hlandhoa 2002
As depicted earlier in Figure [ applying the NPMLE across the full 19- 19-dimensional space computationally prohibitive and produces unreliable estimates, as evidenced by Soloff et al. 20255),
Zhang al. 2024]. Notably, applying SDGM-C the APOGEE data (Figure 4 reveals sharpened astrochemically significant structures. Taken together with our theoretical guarantees, these empirical results showcase that SDGM couples superior denoising fidelity with graceful nigh-dimensional scalability, making it a promising tool for robust inference in complex scientific datasets.
Figure 2: Projected plots of the circle experiment onto the first two dimensions. NPMLE-N (columns & 3) performance significantly degrades when the dimension increases from (top row) to (bottorn sow). proposed SDGM-C (columms & stays efficient when grows.
6 Discussions, Limitations, and Future Work
This work develops a stochastic-proeess reinterpretation of Gaussian location mixtures that bridges classical empirical Bayes inference and modern diffusion-based modeling. Through this lens, we introduce diffusion-type G-modeling objective and a multi-step generalized Tweedie denoising procedure that jointly estimate the latent measure and recover the underlying signals from het-

---

Figure 3: Projected plats of the letter A experiment d 16 onto the first two dimensions. NPMLE and PEM (top row) struggled learning only noisy latent prior, and causing denoised data to collapse within the region between the two logs the A, illustrating over-shrinkage. The batlom row shows our SDGM variants. Notably, even is estimated by SDGM, one step Tweedie demoising (middle of the bottom row) still over-shrinks By contrust. proposed denoising scheme better preserves the A structure,
points contined
of the A.
[SSie Plane
MMlleHMggHHs Plane
Figure 4: Left: Naisy APOGEE observations prujected onto twa 2.D planes. Middle: The demoised data reveals fusilform like or hanana like
structure in the central region and a thin munifold structure in the upper region. Right: The marginal density of fitted suggests strong abundance correlation among chemicals. These results significanlly improve the performance NPMLE and denoising in Figure
croscedastic mixtures. Unlike standard score-based diffusion approaches that directly parameterize the score function, method explicitly models the latent prior distribution G, which is crucial for consistent recovery of the latent law. In contrast to traditional NPMLE solvers that suffer from slow convergence and over-shrinkage bias, framework achieves comparable improved denoising and prior recovery performance, with better empirical scalability in high-dimensional regimes.
Limitations. Our theoretical analysis builds assumptions comparable to those in prior NPMLE studies [Saha and Guntuboyina, 2020, [Soloff et al., 2025) and attains a near-parametric rate of n (log These results rely on standard regularity conditions, including compact support of the latent measure and well-conditioned covariance matrices, which are common in NPMLE heavy-tailed ill-conditioned settings. Although the rate matches the minimax lower bound in Hellinger distance for Gaussian mixture estimation [Kim) and Guntuboyina 2022]. the dependence on the ambient dimension remains unfavorable and limits its practicality in very high dimensions. This reflects the intrinsic difficulty of estimating a d-dimensional latent measure and recovering each 6; from a single noisy sample X,. Since our experiments show that the boree-bsed G- G-modeling approach is more robust and computationally scalable than existing NPMLE solvers, and diffusion models have been observed to capture low- dimensional data manifolds [Pidstrigach] 2022
theory but restrict general applicability
promising future direction is to derive rates d of the manifold supporting G*, following current solution methods rely on heuristic sampling and approximations that do not guarantee convergence to global minima of the difinttedimensiinall nonconvex objective. The formulation in Section 2 was primarily chosen ensure theoretical tractability under the heteroscedastic setting and is not yet optimized for large-scale computation. Future work should aim bridge this gap between theory and implementation by refining both the modeling formulation and algorithmic design. supported by further analysis of optimization and approximation errors within the empirical Bayes framework.
depending only on the intrinsic dimension dp < assumptions similar to Yakovlev and Puchkin (2255] Algorithmically.
Broader implications. The line of work suggests that integrating diffusion- based generative modeling with nonparametric empirical Bayes principles can yield unified approach statistical deconvolu- tion, manifold recovery, and inverse inference under heteroscedastic noise. Future investigations can further extend this framework toward manifold-aware score estimation and other related directions that enhance its theoretical scope and practical applicability.
10

---

Acknowledgment
The authors are partially supported by the National Science Foundation Division of Mathematical Sciences [Grant DMS-2416251 and the National Institutes of Health [Grant IROICA287413-01]
References
Brandon C Kelly. Measurement error models in astronomy. In Statistical challenges in modern astronomy V. pages 147-162. Springer, 2012.
Michael G Akritas and Matthew A Bershady. Linear regression for astronomical data with measure- ment errors and intrinsic scatter. Astrophysical Journal, 470:706, 1996.
Lauren Anderson, David W Hogg, Boris Leistedt, Adrian M Sheee-Weelan, and Jo Bovy. Improving gaia parallax precision with a data-driven model of stars. The Astronomical Journal, 156(4):145, 2018.
Herbert Robbins. Asymptotically subminimax solutions of compound statistical decision problems. In Proceedings of the second Berkeley symposium mathematical statistics and probability, volume 1, pages 131-148. University of California Press Berkeley, 1951.
Herbert Robbins. Some thoughts on empirical Bayes estimation. The Annals of Statistics, pages 713-723, 1983.
Bradley Efron. Two modeling strategies for empirical Bayes estimation. Statistical science: review journal of the Institute of Mathematical Statistics, 29(2):285, 2014.
Herbert Robbins. A generalization of the method of maximum likelihood: Estimating mixing distribution. The Annals of Mathematical Statistics, 21 21(2)3114-315, 1950.
Jack Kiefer and Jacob Wolfowitz. Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters. The Annals of Mathematical Statistics, pages 887-906, 1956.
Wenhua Jiang and Cun-Hui Zhang. General maximum likelihood empirical bayes estimation of normal means. The Annals of Statistics, pages 1647-1684, 2009.
Jiaying Gu and Roger Koenker. On problem of robbins. International Statistical Review, 84(2): 224-244, 2016.
Yury Polyanskiy and Yihong Wu. relf-regul property of nonparametric maximum likelihood estimator in mixture models. arXiv preprint arXiv:2008.08244 082444 2020.
the Bradley Efron. Tweedie's Formula and Selection Bias. Journal the American Statistical Association,
1664496):160-11644 2011.
Sujayam Saha and Adityanand Guntuboyina. On the nonparametric maximum likelihood estimator for gaussian location mixture densities with application to gaussian denoising. The Annals of Staristics, (8(2):338-763 2020.
Jake A Soloff, Adityanand Guntuboyina, and Bodhisattva Sen. Multivariate, heteroscedastic empirical bayes via nonparametric maximum likelihood. Journal of the Royal Statistical Society Series B: Statistical Methodology, 87(1):1-32, 2025.
Nan Laird. Nonparametric maximum likelihood estimation of a mixing distribution. Journal of the American Statistical Association, 73(564):805-811, 1978.
Roger Koenker and Ivan Mizera. Convex optimization, shape constraints, compound decisions, and empirical bayes rules. Journal of the American Statistical Association, 109(506):674-65,, 2014.
Long Feng and Lee H Dicker. Approximate nonparametric maximum likelihood for mixture models: convex optimization approach to fitting arbitrary multivariate mixing distributions. Computa- tional Statistics & Data Analysis, 122:80-91, 2018.

---

Danial Lashkari and Polina Golland. Convex clustering with exemplar based models. Advances in neural information processing systems, 20, 2007.
Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their
Applications, 12(3):313-326, 1982.
Dimitra Maoutsa, Sebastian Reich, and Manfred Opper. Interacting particle solutions fokker-
planck equations through gradient-log- lotfiest-logddessitty estimation. Entropy, 22(8):802, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben
Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021a.
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computa- lion, 1367):1611-1677, 2011.
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.
Tianshuo Xu, Peng Mi, Ruilin Wang, and Yingcong Chen. Towards faster training of diffusion models: An inspiration of a consistency phenomenon. arXiv preprint arXiv:2404.0P9 072444447946 2024.
Chaoyue Liu, Libin Zhu, and Mikhail Belkin. Loss landscapes and optimization in peeeptaramteriized non-linear systems and neural networks. Applied and Computational Harmonic Analysis, 59: 85-116, 2022.
Khashayar Gatmiry, Jonathan Kelner, and Holden Lee. Learning mixtures of gaussians using diffusion models. arXiv preprint urXiv:2404.18800 2024.
Sitan Chen, Vasilis Kontonis, and Kulin Shah. Learning general gaussian mixtures with efficient score matching. arXiv preprint 24X::244441883 2024.
Kulin Shah, Sitan Chen, and Adam Klivans. Learning mixtures of gaussians using the ddpm objective. Advances in Neural Information Processing Systems, 36: 19636-99649, 2023.
Sulagna Ghosh, Nikolaos Ignatiadis, Frederic Koehler, and Amber Lee. Stein's unbiased risk estimate and Hyvärinen's ' score matching. arXiv preprint arXiv:2502..0112. 2025.
Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian. Image denoising with block-matching and 3d filtering. In Image processing: algorithmss and systems, neural networks, and machine learning, volume 6064, pages 354-365. SPIE, 2006.
Kwanyoung Kim and Jong Chul Ye. Noise2score: Tweedie's approach to self-supervised image denoising without clean images. Advances in Neural Information Processing Systems, 34:864-874, 2021.
Yangjing Zhang, Ying Cui, Bodhisattva Sen, and Kim-Chuan Toh. On efficient and scalable computa- tion of the nonparametric maximum likelihood estimator in mixture models. Journal of Machine Learning Research, 25(8):1-46, 2024.
Konstantin Yakovlev and Nikita Puchkin. Generalization error bound for denoising score matching under relaxed manifold assumption. arXiv preprint a 13662, 2025.
Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypoth- esis. Transactions on Machine Learning Research, 2022. ISSN 2835-8856. URL https Expert Certification.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2008
Andre Wibisono, Yihong Wu, and Kaylee Yingxi Yang. Optimal score estimation via empirical bayes smoothing. The Thirty Seventh Annual Conference on Learning Theory, pages 4958 4958-4991. PMLR, 2024.
12

---

Jérôme Dedecker and Bertrand Michel. Minimax rates of convergence for wasserstein deconvolution with supersmooth errors in any dimension. Journal of Multivariate Analysis, 122:278-22 291,, 2013.
Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows.
International Conference on Machine Learning, 2015.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. In International Conference Learning Representations, 2017.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. In Advances in Neural Information Processing Systems, 2019.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios, nflows: normalizing flows in PyTorch, November 2020. URL
42962877
Yandi Shen and Yihong Wu. Empirical bayes estimation: When does g-modeling beat f f-modeling in theory (and in practice)? arXiv preprint arXiv:2211. 2022.
Gilles Celeux, Stéphane Chrétien, Florence Forbes, and Abdallah Mkhadri. A component-wise EM algorithm for mixtures. Journal of Computational and Graphical Statistics, 10(4))697-712, 2001.
Alexander Ihler, Erik Sudderth, William Freeman, and Alan Willsky. Efficient multiscale sampling from products of gaussian mixtures. Advances in Neural Information Processing Systems, 16, 2003.
Yangjing Zhang, Ying Cui, Bodhisattva Sen, and Kim-Chuan Toh. On efficient and scalable computa- tion of the nonparametric maximum likelihood estimator in mixture models. arXiv:2208.075144 2022.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with sinkhorn diver- gences. In International Conference on Artificial Intelligence and Statistics, pages 1608-1617. PMLR, 2018.
Pascal Vincent and Yoshua Bengio. Manifold parzen windows. Advances in neural information processing systems, 15, 2002.
Abolfathi, D. S. Aguado, G. Aguilar, C. Allende Prieto, A. Almeida, T. Ananna, Anders, F. Anderson, B. H. Andrews, B. Anguiano, and et al. The Fourteenth Data Release of the Sloan Digital Sky Survey: First Spectroscopic Data from the Extended Baryon Oscillation Spectroscopic Survey and from the Second Phase of the Apache Point Observatory Galactic Evolution Experiment. The Astrophysical Journal Supplement Series, 235(2):42, 2018.
doi:
S. R. Majewski, R. P. Schiavon, P. M. Frinchaboy. C. Allende Prieto, and et al. The Apache Point Observatory Galactic Evolution Experiment (APOGEE). The Astronomical Journal, 154(3)994, 2017. doi: 10.3847/1538-333 /aa-5488888///8844
Bridget Rateliffe, Melissa K Ness, Kathryn V Johnston, and Bodhisattva Sen. Tracing the assembly of the milky way's is disk through abundance clustering. The Astrophysical Journal, 900, 2020.
Yuan-Sen Ting, K. C. Freeman, C. Kobayashi, G. M. De Silva, and J. Gland-Hawthone Chemical tagging with APOGEE: The role of dimensionality and metallicity. Monthly Notices of the Royal Astronomical Society, 227:882-893, 2012. doi: 10. 10.1111/1 1/j. http///1.....1. 0666-2966.2012.22002..
David H. Weinberg, Jon A. Holtzman, Sten Hasselquist, Mclissa K. Ness, and et al. The Chemical Evolution of the Milky Way from APOGEE Abundances of a-FFlements and Iron. The Astrophysical Journal, 876:65, 2019. doi: 01444471538-4357/abbl-
Ken Freeman and Joss BlwnhoH The new galaxy: Signatures of its formation. Annual Review of Astronomy and Astrophysics, 40(1):487-537, 2002.
13

---

Arlene KH Kim and Adityanand Guntuboyina. Minimax bounds for estimating multivariate gaussian location mixtures. Electronic Journal of Statistics, 16(1):1461-1484, 2022.
Jakiw Pidstrigach. Score-based generative models detect manifolds. Advances in Neural Information Processing Systems, 35:35852-35865, 2022.
Lanfeng Pan, Yehua Li, Kevin He, Yanming Li, and Yi Li. Generalized linear mixed models with gaussian mixture random effects: Inference and application. Journal of multivariate analysis, 175: 104555, 2020.
Saharon Rosset and Ryan J Tibshirani. From fixed-x random-x regression: Bias-variance decom- positions, covariance penalties, and prediction eπro estimation. Journal of the American Statistical Association, 2020.
Peter Bartlett, Olivier Bousquet, and Shahar Mendelson. Local Rademacher complexities. The Annals of Statistics, 13(4))1497 1537, 2005. doi: 10.121400090566650000000222222 URL https //doi.org//l /0/ 02111//00000000000000000000
Yunwen Lei, Lixin Ding, and Yingzhou Bi. Local rademacher complexity bounds based on covering numbers. Neurocomputing, 218:320-330, 2016.
Charles M Stein. Estimation of the mean of a multivariate normal distribution. The annals of Statistics, pages 133-11551, 1981.
Bernt Oksendal. Stochastic differential equations: an introduction with applications. Springer Science & Business Media, 2013.
Ricky T. Q. Chen. torchdiffeq. 2018. URL htt.:tttee//ghhhbbo.oorrr/r//eer
14

---

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]
Our abstract and introduction set out an overview of our methods, the the- oretical guarantees, and empirical improvements that are detailed in Section and respectively. We clearly introduce the scope of study and background context in the first few paragraphs of the introduction.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It is fine 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 discussed several non-positive implications with our theoretical results following Theorem [
We also discussed computational challenges of our proposal in Section 4, with efforts to alleviate them.
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 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 (e.g., independence assumptions, noiseless settings, model well-specification, 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 a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect 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-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency 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 fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren'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 not penalize honesty concerning limitations.
3. Theory assumptions and proofs
15

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: We detailed the most significant assumptions in the main text and presented proof sketch. We aim to present complete technical derivations in our supplemental.
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 in the statement of any theorems.
The proof' can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix 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 it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Our major contribution a statistical modeling perspective together with an optimization objective that can be implemented with various approaches. We discuss our approaches to approximately solve our proposed model. We detailed datasets, metries, base- lines and model variants in Section 5, and pointer to the supplemental for implementation details, code, and more results.
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, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or 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 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 all 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 the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution 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 instructions for how to construct the dataset).
16

---

(d)
The authors should provide scripts to reproduce all experimental results for the new reproducible, they 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 be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: 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?
Answer: [Yes]
Justification: We detailed the complete synthetic data preparation process in the main text. For the APOGEE data we used, is publicly available at hthttp://daat sdss org/ sas/
dr14/appogeevva/appogee rc/cat/ We provided the code together with usages in the supplemental
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.com cc/
Sobllic/gides/Cee httpe///oooosbbiiisioooooooooooo
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/pips..cc/publc//gune http://uiiieeeseeeeee /
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,
proposed method and baselines. If only subset of experiments 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 supplemental material (appended the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Dataset generation procedures, noise levels, model variants and evaluation protocols are given in Section 5. We will put the architecture and hyperparameters for the continuous normalizing flow model in the supplemental.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary 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]
17

---

Justification:
Table I reports mean and standard deviation over 5 independent runs for each metric, if the method is run under stochastic optimization methods.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
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, 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 or the standard error of the mean.
is OK to report I-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 to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables plots, The authors should explain in 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 to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We plan detail these details together with other implementation details in the supplemental.
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 of 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 failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: hett://nnssissse wwweurpps.ccppuulicc//ooooo httpppblii//tthcccuuueeeeel.eeee
Answer: [Yes]
Justification:
We acknowledge that we have reviewed the Code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS 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 (c.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
18

---

Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
We discussed the impact of the problem in scientific domains. We demonstrated that our framework help to mitigate challenges in traditional methods, which has the potential to benefit future scientific inference tasks.
Justification:
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 a 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 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 to attacks, mechanisms for monitoring misuse, mechanisms 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 or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work does not release high-risk model or dataset.
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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We cited the source of APOGEE data and relevant papers. The data is publicly available hithtps:// /www. sdss4. org/
19

---

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 or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
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.
If assets released, the license, copyright information, and terms of in the package should be provided. For popular datasets, paperswi thcode. com/ datasets 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 it has changed) should be provided.
If this information is 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: Our paper does not release any new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model 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 is 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. Crowdsourcing 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 as details about compensation (if any)?
Answer: [NA]
Justification:
The work involves no human- human-subject data collection.
Guidelines:
The answer NA means that the paper does 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 the main paper.
According to the NeurIPS Code Ethics, workers involved in data collection, curation, 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 an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
20

---

Justification:
No numan-subject research is conducted.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor 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.
recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: 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, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality the research, declaration is not required.
Answer: [NA]
Justification:
The core original method are developed without any LLMs.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy https///eeeeessccc/oo//eeeeee//e for what should or should not be described.
21

---

Supplemental Material for "Score-Based Diffusion Modeling for Nonparametric Empirical Bayes in Heteroscedastic Gaussian Mixtures"
A
Code
Our code is available at/https://gtthbb com/ wwwo///hhoogg///gnnnoooooo neegggm-neurippp git
B Proofs of Results in Section 3
B.1 Formal statements of Theorem 2
Assumption 2. For the SDEs defined in assuae
(latent measure) E P([-M, M| ) for some radius M E (0,
2. (covariances) gla 3 E; for all ,....,n.
3. (diffusion coefficients) For all = ,,..., n. 9i R4 - R is a continuous function. Let a,(t) S/ ( ds. We require the normalization condition (1(1) 1, and additionally the well-behavion condition: fix a T > 1, then for every E [1, T), there exists (g, 9) 0 such that g S (t)) S g.
Theorem 3. Provided conditions in Assumption 2 les (la, T) satisfy to S T and
log S $/, (((()tt & { (1)d
Let Gn be an optimal solution to Objective constrained over the measure class P([- M, M] 7
then provided n sufficiently large, with probability at least
for score estimation:
E3 to. m (96
119
=
E Bengge V log loggo!
- (x) V log 4
(t)
(x) 1at
(Ca.M(o.)))(a.) (log ( 2d433
for density estimation at = to:
E 11)
(ta)) SC SCa.M
(log n) 2d+3,
for the deconvolution risk,
EW2
1 M (E:)
-1/2
>> (E.) -2/2 #/22
log
B.2 Control of Fisher risk
Our arguments are built upon the "random-X" regime with fixed design model (i.e., sample size latent prior G* and covariances {E,)r are fixed). Conditioned on the fixed design, we can resample the dataset X,, [X.)%, from stochastic processes {G" (0, E,) This view
is similar to, e.g., |an et al. [2200] where the observed data as treated random realizations from fixed design model; see also
Rosset and Tibshirani 2020 for discussions on random versus fixed-X
statistical viewpoints.
22

---

We recall the following definitions of the rescaled mixture laws and its empirical estimator:
1 - Ki(-18)
(8)
4G
fa
Ex. n 19h
9G
where Ga,, = ) 1/2 #G denotes the pushforward measure under linear map induced by the covariance structure Σ at time t, and the rescaled transition kernel that perturbs the i-th data pair (X,, E,) is:
Kcl- N ((a, (1 (1- (1-aa(t)-+)1))
(*) We abbreviate the score function for 9G and ( respectively as:
V logae
(1) V log ql0
By definition, the cumulated Fisher divergence E ES T (4a- Ma between G and an arbitrary deterministic measure G is naturally related to the population score matching risk defined by:
J.(G) :- ー Σ I" oe (t)
[1188(0) - 3/ (x)
dt.
10
It is clear that the Fisher risk can be upper bounded by J. via rescaling of diffusion and noise coefficients:
(E,) J.GG).
sup sEln).eElto.T!
Notice however that J. (G) needs not be the expectation of the finite-n empirical risk defined in eq.(1) (T]. i.e., Ex, (G) * (G) in general. This is due to that (1) the score estimator is biased ratio estimator (though the density QG- and gradient is unbiasedly estimated): 1xx, 8, * &G-, and (2) both the integrand functional 8,, and the sampling law K, are X,edeependent thus they cannot be trivially decoupled.
Thus, we further define a decoupled risk with respect to population-level score sg. and data-dependent sampling law K-(-)X)):
B. = 1 Σ n1
|18 - (a) dt.
(t)
1m1
Then it is clear that Ex, [.) J. (G) for every deterministic G, because E(K)(+()) = for
For the random estimator G., the optimality condition of Jn (G..) implies a basic inequality that allows us link J. and Jn through J., detail below.
Lemma (Basic inequality). Every optimal solution G, to satisfies:
J.(G.,) S 4 Ja(").
(15)
Proof of Lemma 7 (i, ()) to ease the integral representation J..
We start with the following definition that presents joint view over the tuple
Definition (Joint law over index, time and sample). Fix the data {(x,, Ei) and set
n
Σ
(t) dt.
Define the probability measure Pn on [l....,n} nn x [to. T x Rd by
(t)
Pú k.t dt, dx)
n Z 9 (11-0((0))))
dt dx.
23

---

Equivalently:
Pr(i k}
g Ho sf(() du'
i k) has density
(x i k, 1) Ki((k).
With the above definitions, we first note that we can rewrite the objective 11 by:
[n(G) = B()) (x) (0) (x)
2 dl
x
=Z E (1t
|1 1/ (z) - (z)
By the optimality condition, provided G E P(R4) is minimizer of 11 and G" is feasible, we have
Tn (G) In")"). Using the quadratic expansion:
||a c|L |6 clB ||e b|b - 2(b - c) [a - b) s 0,
where we let G b s and = (1) we obtain:
E ((tta) |1 l I| s 2 (
[00 -100-1107 // E -
(16)
By Cauchy- Cauchy-Schwarz Schwarz inequality applied to the inner product term (T6) and rearranging terms,
E (L..)))PP
|| (x)
2 (x))
[
4 (ML.)
[1 ( (x) (x)
E.
which recovers the claim by rewriting the integrals back Jn and J.
By the basic inequality eq. 15
can now relate the empirical risk to the population risk follows:
= -
J..)...
(17)
The first term is the uniform deviation between the true and empirical J. over the feasible set of G, and the second term most four times the score empirical regret (G*) conditioned on the observed data X,, Below we control the two terms in expectation respectively.
B.2.1 Control of (G))
We first have:
T
n T1
(z)) 2 -000)())))
[n()
gi (t) 1
dt
[1-200)
(00 (x) | dt.
max (ot Amax (EA)E g0 iEin|
To control this term, we integrate and generalize tools from Wibisono et al. (2024 Saha and Guntuboyina (202)] prove our next lemma.
Lemma 2 (Control of score empirical regret). For each t [to, where (to, T) are specified in Theorem
((0) (x) 12] < Ca
1 nlog
Ex,
E
(0)
$ (x)
(18)
We now build up the necessary toolbox prove the above lemma. Following Jiang and Zhang 2009) Saha and Guntuboyina [2020), it is beneficial to consider the following regularized score rule. Given
an arbitrary density p. define its associated regularized score function with parameter > 0
Vp
PVc

---

where p V e denotes max c}. For fixed t, when it is clear in context, we now drop the index (t) the score. Expand the empirical regret via regularization:
Egg.
= 3Ey. [||sc- sc:- 112
|s --.--s,,11 s,,11
11 s
As useful in the following proof, we first define general class of heteroscedastic mixtures, in which both the base measures G and the Gaussian smoothing kernels can be heteroscodastic across components. We then show that their associated regularized score function is uniformly bounded. The definition will be useful to characterize the behavior of the empirical score 8n, in which each component is different Dirac measure "(a,E.)-1/2x smoothed personalized transition kernel (0, K,) with K; = (1 - aj "a.
Definition 2 (General heteroscedastic mixtures). Let {Ki) be positive definite covariance matrices such that kla K₁ 3 klg. Let G/n) {G:) be collection of probability measures R& Define the heteroscedastic mixture p. E G. * (0, l.K,) for bandwidth parameter
h > O.
Lemma (Uniform bounding of regularized scores), Fix a mixture density defined Definition For every E Rd
(2Thk)
log
Pm Ve
hk
2
for 0 < <E< (2thk) e/26-1/2 1/2
Proof. For each component P: G. * NO, hK,) in the mixture, consider the transformation such that if X Pi. let X; (hK,) 1/2 X; and let p. = G. (0 Ia) be the transformed measure governing X, Then we have the following equations:
(x(x) = |4 1 ((xK))-1222)
(p.(2)= = |/K | ((hK))
For Ps being convolution of measure with standard Gaussian kernel, Saha and Guntuboyina 2020
Lemma F.1] directly applies:
12T)
log
Pi
1122 [(hK,)
(2t))
-d
(2thk) -d
log
log
Pi
Pr
lk
m2
lik
P
We observe that the function
//oloooo -x is concave E (0, VB),
thus by Jensen's inequality:
1
(2mhk)
|VPlll Supowann
log
lk
P2
and
I PeGmAAll
Vp.l
(2Thk)
log
PGpn,,
PG
The claim then follows from same argument as in proof of (F:2) in Lemma F.1 of Saha and Guntuboyina [2020) to relate the regularized rule with the unregularized bound derived above under the two cases h and P e.
We now turn our attention to tools for bounding the discrepancy between regularized scores:
By [1l =,,113]
via their corresponding Hellinger distance. We rely specializing Theorem E. of Saha and Guntuboyina (020) and Lemma of Wibisono al. 2024 the following form that is closely
related to setting with heteroscedastic mixture.
25

---

Lemma 4 (Bayes discrepancy of regularized scores). For two arbitrary n-collections of probabiliry measures (pn) and [[n] on R Let the mixtures P and PH be defined in Definition then
[2]
[2]
Cd max/ { {1os. (2Thh)) dd/2 -d/2
$
max
log ((GG PMWa)aa] P 53
(porn) ( PHn))
hk
E 2//2-112. 1/2.
if hk & and < <ES (2thk)
Proof of Lemma 4
Examining the proof of Saha and Guntuboyina (2020 Theorem E.1], we have:
|| - s siHm))bb
1
1/2
$ Ti + T2.
with
d 13 S
7 T - 25'
and
A
4
1=1
where
[P)a
A /
(0h 10xt (po - 1h)
r
sup ZER4 De H(a)
(o)) Ve
PG V PH Ve
With integration by parts, it is shown there that A's satisfy a recursion that for h 2 1,
2 S 2T A
With generalized Lemma F:2 of Saha and Guntuboyina 2020 and some work around the recursion, we can deduce that for every a Z V(2k 1)/(k)):
A 2(2mhk)
4/2
2 (a² (22/hk) (hk 13/2 hkaa
2* 53 (PG PH
PHm h)) + V
(26-1122
and further with & = ko+12 2 we have
hk 22 and taking oN Th max (222i + I,, / log (prm. Po,n)]
2
e hkaa < (P) PHHp))h)
and
A.1 S C max {2kg Y, e(1 + aa S(mm, Po.h).
Plugging in the uniform upper bound on Y by Lemma 3 and some basic algebra around the choice of ku and merging of terms yield the desired bound. Because the last algebraic steps exactly follow from [aha and Guntuboyina 2020 Proof of Theorem E. we refer our readers to their work for more details.
The above result implies that we can further reduce the discrepancy between regularized scores to bounding the Hellinger distance between (G- and its empirical estimate 9n- We recall the following result from Wibisono et al. (2024] that controls the Hellinger convergence rate of smoothed empirical distributions with the vanilla radial Gaussian kernel.
Lemma 5 ([Wibisono et al. 2024 Lemmas 2 & 31). Let p be an a-subgaussian measure on IRd
whose score s(x) V log p(x) L-Lipschitz. Let ph p sph where sp is the density of N(0, hla). Then for every E Rd
p(x) S exp exp(d (dLh/)).
(19)
(v(
Let Pn be size-n empirical measure and Pa, = Pn * N(0, Ila). min a?}.
Then provided h s
1
Ca" log d/2
4d
(p pn)] <
+
71. h
n
where C' > 0 is a universal constant.
26

---

We extend the above lemma for empirical heteroscedastic kernel densities.
Corollary 1 (Hellinger convergence of heteroscedastically smoothed empirical measures). Let {Yi)== be independent random variables, each drawn from an a-subgaussian, L--lggsmooth mea-
sure Ps on IR, Consider the heteroscedasttcc mixture Pn with density: Pn, (2)) == +/4 + Ia
Yi) and denote its population mixture Pa,h Elpnn,]
Ei= ( hK;). Assume each
kla K; kT and hk min(a?, , 1/21. Then,
Cg n
logg h
d/2
4d
E (p
+
(20)
n
for a constant C that depends on a, k.
Proof of Corollary 7 By proof of Lemma in Wibisono al. 2024 let B Rdd
$ (((p,, Pn, (x)) de + Bc Pa,h
(Pm,hr (P,...P...)))
(,,)(z)
Fix each I, then for each i-th component of the empirical mixture shk, (x - Yi), by their trick:
Epak, (x Y,)² -tahkK, Yi).
Thus
1
Varlpn,
Varlsax, (x Yil
(2 - VJ) - Ex v2)3) (4Thk) 4/2
n p=.)//(1))
t=1
Combining with above,
< (4thk)-4/22 E 50 (po,h,
/ P..n/2la) dx + Pa.k-
B (x.
B
Notice that P.,a is subgaussian with parameter region with center Ep,, [Y] and radius //lo n, standard concentration inequalities for multivariate subgaussian variables yield S Pt,h S and for the first term we can bound:
va² Q + hk. By letting B be the box
[ Pa,h/2(() dx < (B) sup Pa,n/2(z) (Ao² log n) 4/2 (x)
P=,A/((z)
sup x
x (r,)(x))
P-,a
where the last term can be bounded by applying immediate extension from (19 and Lemma of Wibisono et al. (224). For each component the mixture, consider the rotation Pi,
K; 172 #p: (O, hla). Then notice V VlogPi. =K,
K./² Pi.h- For all h > 0 and R. let
K; x such that:
log
P&d. (I:)
log
Pi.h K 1/2
PL2h
Pi.zh
< llg ( log Ph (ī vhZil
-E-A((0.La)
(-Z, V log (T uZ)iddu
<E f 1211 log uZ) log (T)llldu S kLdh, which implies pa. h. ( (x)) S exp(k.
Combining the results, we obtain:
402 + hk) log n/ 4/2
4d
Ca / 4a (/22/h1 /h + To
d/2
4d
E (p <
log
n
k
n
27

---

Lastly, we study the discrepancy due to regularizing 9n, and gG-, which can be controlled via adapting results from Saha and Guntuboyina
[020]. Specifically, we recall the following result:
Lemma 6 (Saha and Guntuboyina 2220 Lemma 4.3)). For probability measure G E P(RR)) and
let E (0, (2x)- --12j. Let M(e)
log log( (2xx (2).
Then for every compact set S Rd.
A(G, = (1-Jah
Ja. fG.1.Ve Ve
IV 2 S CC
fG.ua
where N(6,S) S) is the regular Euclidean ô-covering of set S.
Intuitively, the score and its regularized counterpart differ only in the tails.
Corollary 2. For an arbitrary n-mixture of probability measures Giul 1 Di= G. E P( and n Let E E (0, (22
let the heteroscedastic mixture PG defined as in Defimition 2 and define L(e)
log( (2mhk))n222) For compact sets S SI,...,,,
Sp S C
4 h (.((5))
E [1 - o m
C.E.E n &C
Σ
Proof of Corollary 2 The strategy prove this slightly generalized version is similar proof of Lemma 3 the Proof of Lemma F.7 in Saha and Guntuboyina (2020) and proofs in Section D.2 of Soloff et al. 2025). Specifically, for every fixed x. we can let G, denote the distribution of
(hK,)-1720, where 0: ~ G; and write = (hK,)-1/22, This allow us to rescale the terms in the
score function:
(z))) =
|hK 1/2 f8.. (m:).
i=1
[pG
|KK, 1/2,,11// fovaaaaa fān (zi).
n
By convexity,
VP (1)
1 n
||KK,1 ||2KK11//((K,))--11101>))) -12 f8 14 (I:))1
PG (x))
[hKi 1/2 1/2 (xi)
< (d/2++112+11) 1
IV VVa.de (2:2))
S (5:)
and it is trivial to verify that for every E [n]]
MT) 2 (li)) 4/2 -d/2 (x.)
PG PG
$ (1 1-
(( (2i)
2
PG (x)
Sa (2x) V m(hTydd/26
Thus,
S6 16, 1 111-
PG ((())
2 VPGGa (2)]12
Epg n [|lsamn/
dx
pG ( (xx) Ve
((r))
(1:)
2
|1V fö. (c:)l1?
$ C. do nh Σ
dx.
Se (2x) V n(hit)d/2e)
fā,,1
= -CaKT nh 7=1
A ((hk))aa)
28

---

where C,,,ET = K 4/2 /h /1./2211 Now, Lemma directly applies to each summand. Provided 1 hu log (2mhR)"n²??) then for compact sets S i
(0, (2nhh) 11/20-1/22 11 and L(c) V-
1. ,,. n,
E Pe mn G0
f. II C. nh Σ (44/ dG.,)))
where we applied the covering number bound and the change-of-variables identity for pushforward measures under linear maps to re-express the terms defined on and S
We are now finally ready to integrate the above pieces together to prove Lemma] Proof of Lemma
Provided
n
On
x)(-1))
Σ
6
(a. 1/2x, *N(0 (1 -a,
n1
i=1
9 n
(
n
fo
1
Σ
for
*N(0,, a -1-a0(1-a;)11
72
1=1
where each Gū == and qa can be viewed as specialization of the mixtures in Definition 2
1//2 and 2: = (z)E,, (Le gf (s)ds)E,, It highlights that qu We can easily verify that is at most / 4 + I-subgaussian and 1-log-smooth. By construction, 1- a-a," e fo t=1 11.
By the earlier decomposition,
TEx,, F.g.. [|sc- - 8n [|] S 3R.
||5G- - .6. 61] + +3Ex, E.g. x. [1 s/11 + 1|8 [1
For the first term, the test function is deterministic with respect to the data, thus we can reduce the expectation to its population quantity. On the other hand, for terms involving Sn, it is random function of the data samples, thus we evaluate first the inner expectation, conditioned on a fixed data.
The first term can be controlled via direct application of Corollary 2
with (2T) (2m2-d/2e e//2-112/ /n
and the sets S. = ) 1/20 [-M, so (S,) = Then we have:
Eg. ||sa- - -6 12] SCalsup N (4, 1 s
(1
MVd
1
2(min; az()))
SCa.Mo
The last term can be controlled similarly. Standard Gaussian concentration yields that: provided 2 2d, with probability least -1 max || I Tn M + /logg())). Now let the compact set in Corollary 2 be S. == ) z x [-rn, r..]d. Conditioned the high- probability event, each (1/2x, is compactly supported on S, and thus satisfies G,(ST) = 0:
otherwise GS,)) < is uniformly bounded. Let -2/2 - 1/2 /n and we obtain:
(logn)d
Ex, Bgg. x, [1Im. - [1]
sup N(4/L(), iE(m)
Lastly for the Bayes discrepancy between regularized scores, by Lemma Corollary [ and the same strategy reduce the expected Hellinger distance Wibisono et al. 2004 Lemma 4),
Ex, E xx., [116-- 11
<C,,EE., [max ( ()ogn)* log (gg.. qo-)13 (qn. 9c-
(log flogn)d/2+3
<CA,
29

---

Combining the terms, we obtain the desired result that for all t E [8v, [ T),
1 8 (x) s Ca.. (logm))+3
Ex,
E snggh (t)
(c)
B.2.2 Control of
We control the term J. (Gn) T. (G.) in 17 via empirical process theory together with the regularization trick introduced above. In particular, we erucially rely on local Rademacher complexity analysis following Bartlett al. 22005] Lei et al. [2016]. Localization also appears essential obtain
sharp learning rate empirical Bayes guarantees [Saha and Guntuboyina 2020 Soloff et al. 2025
Stein 1981] Under the uniform bounding assumption, it is without loss of generality to let 9i for all E [n] to ease the notations.
Denote the pointwise loss &i with:
(a G) G): = =1 (z) - (a)
2
and let the regularized point-wise discrepancy et where we write the shorthand notation s6. == V (t) (d6 V e):
:- - 56
We also let:
ZZG), G); B2--No,1- + Z:G). such that for all fixed G, denote X, ==
ha; G) -2
Eih,, Gi] = Bax-saz da | (X:
5.(G)= (00))o0(()) 9(1)2 Σ (Rx G)dt, J.(G) = (G)
Consistent with Bartlett et al. [2005), we write for a function f + IR
Paf = Σ f(X), Pf=Bf(X) Raf 1 Σ (X.).
Here 01.-... Un independent Rademacher variables. For function class F, its Rademacher complexity is ER, = E. [sup]eFR R.
To facilitate metric entropy toolkit derived in Saha and Guntuboyina 2020) Soloff al. 2225 consider the following pseudometric that restricts the localized score discrepancy to compact set St
sup m (G, H) == =E[n] zes
(1(),(()))) ,;),
It is natural to consider the following compact set: A := {o E Rd le (x) 2 e, VG e P[[-M, in which we guarantee 8G $6z for all measures G. Denote
(g, ox) = min minn (I(()), max max (2( [Ejn)
:E[n]
It is straightforward to show that a sufficient radial condition that implies & E A is:
||>||| <r. < (2l0g((2m)) Md
Let B (x E Rdd ||0||| S r.} then B S A. By standard Gaussian concentration, provided
S = (2Tr) 4/2 exp ( (M Vd/g, Vd 2 log
30

---

then
P(X E A) 2 P(X E B) 2 1--n- n VX Ja Vi [n].
(21)
With the above construction, we define the diserepancy function he, [n] R + that restricts the regularized discrepancy to the high-probability set A::
hi x)) -B 9) I (2 (; ZE AH]
Then for arbitrary fixed G. separate:
J.(G) JJ(G)- [" ((0))( Paha,idtt
+ I s(t) 1 Σ e 1[1 - Roo 2 ) 2 E.7] dt
to
1 f o(t)|Phcc - P.ha. + {" g(t)²
Elub, Hd
to
where P and Pn are defined with respect to the product measures Σ>= 8; & Scr and EE= ( defined over the augmented space [n] x Rd. We drop the last term since a/'s are
nonnegative.
Control of [P- Pa.T.. We apply the following results that establishes generalization analysis with local Rademacher complexities:
Lemma (Theorem 3.3 of Bartlett al. 2005)). Let F be class of functions with range in [a, b] and assume that there are some functional T + R. and some constant B such that for every T Varlf] T[ < B Pf. Let yp be sub-root function that satisfies the following properties:
x: 10, x)) - [0, oo) non-decreusing and - =pr)/ /r is non-increasing for >
pp is continuous on [0, oo) and the equation yp(r = r has unique positive solution r". For all r >, r = r 2 p(r)).
If for any r 2 r", vp satisfies,
((r) 2 TER S())
then for any K > 1, and any t > O, with probability at least - exp(-tt)
K
K " + (e2(b - ())(b-a)+00) a) c))c)BK)) E
PJS K-1 Pa. + P
where Co. C1, 02 are numerical constants.
is convenient then apply Lei bounds in terms of metric entropy.
Lei et al. 2016) results that presents the local Rademacher complexity
Lemma 8 (Corollary of Lei al. 2016). Let F. be a function class with sup Il S If F can be finitely covered with respect to the metric || (l12(P.) (f that is there exists positive numbers 5, such that log N Nu,, (u,F, F. F11.111)) == sup, supp, log N (u, F, (1-1l22(P..) S
8 (r/((r//) for any E 7), then for any E (0,
8r lorlog (2 ER F: Pp <r) Cbppr [
log (2-m2/2) n
(22)
n
In particular, the right hand side of (22 defines valid sub-root function whose fixed point satisfies
r* < cô log(n)p/ /n (see Section B.3. of Lei et al. (2016).. Thus, it suffices for to show the function class of interests satisfies the required covering number bounds.
For every fixed I, the (relaxed) function class of interests is He = (ha. : E P(RR)). We have the following:
31

---

For all h E Na, there exists (e( such that 0 S h S v0' I 4supe |1&6 - bc(c). By Lemma
let
bc(c) = 2 So/lo[(11) 2].
The variance condition is satisfies by letting T(h) == Ph Then Var[M] S Thh by definition and T(h) S PP by the uniform bounding above.
The covering number of H. satisfies N (u, H4, 11-112) S ( (zitee P(RR), m For G, H )( and (i, E Pr. Th..- hin < 2h, (<) + Z:G)1/2 1 + 1/2 1(x + 2 E A}}
By the reverse triangle inequality that |||a cll b clll s ||a bll and conditioned on y = x + ZE A:
(u; G) ( H)
(11) s6., (, sa- [1)12,
(1)) (4)(())11, S m (G, H).
Thus for every and Pa.
ha, h S 2b,(c)m, log (u, He, || 112) S log N u
(G. H)
P((R)), m
2br(c)
By the same strategy in Soloff et al. 2025], we have:
(C, So, (2ee)-bbb()
sup
190 (a) 9H (x)
(t)
zEA
sup + 1/2 6/26-1111
| (t)
90 (x)
V (t) (x)112
zEA
So. 12 (26) (006)-''be b.(e)
max sup JC (x) J Sa (x)
(E(n) ZEA
+ o 1/2 1/2 c- max sup iE(n) ZEA
[ (x) - (T)l..
Lemma (Lemmas 5 and 6 of Soloff et al. 2025)). ). Define the semi-norm over S c Rd:
fl.c. fH.ooloc.s == max sup (c) fHe.)a())),
*E[n] #es
lllc. fHH,olv,s max sax V fo....aa (x).. (E(n] aeS aES
and let Sa == fxr E Rdd inf || ||x-
-ol--ll vll S a). For F = (If.,K.ii (G. E P(RR)) kTa 3
{K 3 kla), there exists positive constants Ca and Ca such that for all compact S c Rd. M > 0, and sufficiently small nj such that
d+1
Cd,E log 17
log N1,, F, lloo.s) S Sa)
d+1
Cd.
log (1,, (|v,s) & CaN((a
log
where a -2k log
/ kA1 (2mk)d/21
Apply Lemma to our metric, by letting
u" J. 1/2 (b)(e22 1)
32

---

provided u² sufficiently small, we have log Nu", P(RR), m²) I log N21/2, [l log (1v,,))
Cd g S Ca)
7
Caa 2oo4 ca(bi(e)/2+1) 1) yd+1
where the last step is due to Lemma F.6 of [Saha and Guntuboyina condition for x E A is given by
2020 Since a necessary radial
d(x,
M M
1d)2
inf
|| v112 S (2gg((2))) (022))))d2221)) 4)26-1)= =: (c),
vo
ME(
M
we have that:
C [T((e) + M M / +
3
a Vol(A800/2) (00l(A30/2]))) a CC.
2
SdMaa (Fa(c) VI) log (5(2T)-
d/2
1
This implies that:
u²* log 20)(e)) P(RR),,)) m SaMaa. (Tr()) log C"(b.(e)/2 + 1)b,(e) 3d/2+1
en"
By choosing
o(n-.), we have be(en) x vlog n and (i((e)) x logn, provided > Nu.M,.... sufficiently large such that terms involving n dominate constant terms. Thus, one can apply Lemma 7 (with K = 2) and Lemma (with xx (br(e) + ())x d x (Fr(e) V 1)), P = 3d/2+ + 1) to show that for arbitrary k 2 1 and G E P([-M, M, M) with probability at least 1 n
Pha.- - Patiaa S Patica+ + 2(Tr(E) V 1)ad log(n) log(n)3d/++1 3d/2+1
+ c2be(e) log((n)
be(c)²
n
n
- - g(t) ((t))((pi))) (Pha.. Paha_. ldt s (....) + O (logn)2d
)
since Tio, (1, S h1,e G) and the time integral of Paht, (x;; G) recovers G.( by definition. Control of tail. We remain to bound for every and i: B (er (r E A
=B((ys, (z; G) (x; G)) E A + ElUi, G)){x A9).
The second term is trivially bounded via uniform bounding (e( of every regularized score functional gie and eq. (21) that IPAA)) S n² provided en chosen above. Thus,
E Blvf, Gz E A & (e))n=
By standard algebraic operations via triangle inequality, let 51 maxee(a,)) ||c cf|| and 52 (a,.e(a,b) ||c-||, la bll - ||a" - b5l + l|a - a² || + 1bb - b'l = ||a² - b5ll + 20
|||a b112 - -Ilaa 1 <|I|a a- bll ||a" - 6 ||| (Ila -1l + ||a - I))
<28 (2) b51 + 261) - 4(0 + 28,
Letting s b
s s a- 5 be 6 . and applying Cauchy-Schward inequality, the first
term in above becomes:
E[(V', ((A,) =E[ [[]
sc (12000 - |1 $ |2 |(44
max 40, [[[[|sH - - + 2||sH siill || [
1 ((-MEP
max 46, Bi - s (c -
HEP([-M..
bx(c) A 1/2
< 45,
A
(ii
n
33

---

where we again used the uniform bounding. The discrepancy A := E| +11l1 can be bounded by Corollary 2 with change of measure. By the Schwarz inequality.
A A / + (Bnn || - ffx 1
where the fourth moment of the score discrepancy can be analogously controlled as the second moment with only the change of the leading constant and the dependency from T.d(e) to (c). The chi- squared distance x" (pllg) f e is jointly convex and we have (f upper bounded by constant. Consequently, given G .. compactly supported, with the choice of e S n -, we have that:
S Ca,Ma
(log n (d+2)/4
Recall by(c) = Vlogn, n, so combining (i) + (ii) yields:
[ g(t)²= Σ Blu (A 2 (logn)(4+2)4
Combining the high-probabilityy part and the tail, we can conclude this section by asserting:
(G - -..(Gn) Sd.M, J. (G.) + polylog(n)
By further upper bounding (G.) by In (G*) from the basic inequality and applying Lemma 2 we recover the first claim in Theorem 3 on score estimation guarantee that with probability at least 1 n k for arbitrary 1,
log(n) BJ.(GG) SEJ,(G") + ()24))241122
1 (logn)0(d)
n
n
n
B.3 Control of Hellinger divergence
Our next result controls the KI. divergence between te and 96- by the Fisher risk controlled above. Without loss of generality, assume gt()) g(t) (t) att). Recall that we define:
1 n 71
1 #G N(0,La
Let
(0)
(Σ,) 12
9g
n
V log 46
(1) (x)
V log 46 (x/ (o()).
Va(t)
Lemma 10. For every fixed E P(RR),
KL(g (1)) [1g(
S 3 too
(g (t)
|996
Proof. Denote the path space by Ω C(to, T): Rd). Under the latent measure G and the same conditions on Σ, and g; as in Assumption I define the path law QG to be the unique solution for the following reverse-time SDE on [to, T):
(x -tt))² log de (x di + g(f) dB,,
(t))
TG
(23)
Similarly to Song et al. 2021b Proof of Theorem 1], consider the measurable Markov kernel K: Ω B 0, with K (w, x) dw(())= then we have:
(ta) 4a
QcK = / K(w,))
34

---

that shows that the marginal mixture law qG at to is the push-forward of the aggregated path measure Qc under the kernel K. By the data processing inequality,
KL
||0 (tt
(L K k)) & KL(Qo-1QQa)
In addition, the KL divergence is invariant with the linear rescaling:
KL (to) (
(L (
Since the SDEs associated with Qe and Qa- differ only in the drift coefficient, we can apply the Girsanov theorem [Oksendal 2013 to yield:
kl Qa) T
[ g())²
log aG- V log
9G
(1)
2
dt
g(1)
E
V log log 9G- - V log 96
(e)
(2)
dt
2
a(t)
2
|| log (2) - v log
(0) 96
B Bawsor
dt
2n
2
S 8 ( 11
(1)
Plugging in the first claim that bounds F1.. 7| (9G-1 ll9 inequality that 57² < KT. yields the desired second claim. Lastly, the third claim on the Wasserstein divergence:
(t) in Lemma 10 and use the basic divergence
EW2
(1)- (E.))
1
log
follows from direct application of Theorem 10 of Soloff et al. 2005 with the upper bound of
55² (9G 96 from above.
(to) (ta)
B.4 Reduction with Homoscedastic Noise
Proof of Theorem
When Σ; = E are all equal, We have:
da = f(50)) fa.m = 100 (1500)-1/22)
Vlog Saxg (x) = (5/1))-1/22 log 16 ( (111)))1122)
We have the following equivalence:
8 ( S Amin (
KT = 114h)); W H) VAmax
This shows that Theorem 1 direct corollary of Theorem with only adjustments to the leading constants.
C Implementation Details and Additional Results
Experiments compute resources. We our models and performed all experiments on an internal HPC cluster. Jobs were scheduled GPU-enabled nodes, each using one NVIDIA GeForce GTX 2080 Ti GPU and accompanying Intel Xeon Skylake 6130 @ 2.1 GHz CPUs with 96 GB of RAM per node. Individual training runs require maximal of 2 GPU hours and 4 CPU hours, while
35

---

the denoising runs require maximal of 1 GPU hour and 2 CPU hours. The entire experiment takes approximately 120 hours. These estimates do not include preliminary or failed experiments.
Hyperparameters and architectures. For our implementations, we consider two variants SDGM-C and SDGM-D, for continuous and discrete latent model, respectively. In the continuous case, the flow uses four repeats of the following block: ActNorm - random permutation + MaskedPiece- AaaaaaionnnQuuaaaiiAAAtttttooo
Transform. where each autoregressive transform conditions on all previous dimensions via a 2-layer MADE network (128 hidden features), uses eight linear-tail bins (tail bound 10), enforces minimum bin width/height and derivative of x 10 leaky ReLU activations, residual blocks, and no dropout or batch norm. The base distribution is standard normal distribution. the discrete case, we use a trainable mixture of 2,000 components (initialized from random normals or random draws from the observations) and corresponding randomly initialized weights.
For the baselines, used the official implementation of Dual-ALM provided by the authors of ZZhng al. (2022) at
For NPMLE-N, we use all available training data samples as the fixed supports. Thus, the latent mixture is given by: Σラ=1 w,6x, with wy's trainable. For PEM, we fix a latent mixture of 2,000 components, i.e., G Ej=1 Wr6p, with both {(w;, /3)) 2 optimizable. We randomly sample points from the dataset the initial supports and uniformly initialize the weights.
com/Yangl hwtp//aannnnngggnnnaaa
SDE setting and training objective. Throughout our experiments, we consider the following family of variance-eae exxploiix SDEs:
aY = g(t)_///dB, dB, with y-(1) G"
In particular, we use a shared diffusion coefficient g()) = V²2_ Under this SDE, the marginal distribution at any time 2 is y(()) G N(0, =.()) where the effective time-dependent covariance is D = -1 Σj Our observed data X; correspond to Y,()).
The training objective for our SDGM models is based on DSM to match the conditional score of fā with (1x,. For each data point (X,, Σ,) in training batch of size NB:
1. A time 4 > sampled (specifically, i, drawn uniformly from + e, Tmax) where we use = 10-3 and Tmax 1.5 in our experiments).
2. A perturbed data point, X/. is generated by corrupting X, with additional noise: X = X; Z, where Z; ~ (0, I(ah ) This X is sample from the SDE path time t, conditioned on X, being the state at =
To compute the score of our model's current estimate of the data distribution at time t. we first represent our learned latent prior, G.
For SDGM-C, ō represented by neural spline flow. We obtain a discrete approximation by drawing = 2, 000 samples {1y} from this flow and assigning uniform weights w 1/m.
For SDGM-D, Ō Σイ wydng, where the support points {uy}y= and their corresponding logits (which determine weights wy via softmax) are directly trainable parameters, with m = 200.
The score of the marginal distribution conditional on i-th SDE at time 4 implied by ă and noise covariance 2x(()) denoted V log (x This score is computed as Σjn P(,,), ((())) X(), where fax(o()) En wyNN I(")) and XX, I(t.)) are the posterior probabilities of the mixture components given the perturbed data X/. The objective function minimizes the expected squared difference between this model score and -(>)(e) -
the true score of the perturbation kernel, Vx; log Pie, (X X.)
Z,. For batch,
36

---

the loss is:
N2 Σ 9 (4) NB
[(G)
log (xx( log (x)
All our experiments are performed with the importance sampling strategy as Song et (022bb that further reduce the variance in estimating the time integral by directly sample t; according p(t) x gt) The optimization is performed using the Adam optimizer. For the experiments presented, we used learning rate of 5 10-5 (unless specified otherwise for particular experiments), 31 0.9, Adam c c 10-8 no weight decay (weight decay 0.0), and gradient norm clipping at 1.0.
Denoising. For denoising, we run the "dopri5" ODE solver implemented in torchdiffeg [Chen] 2018) with absolute and relative tolerance fixed at 10-1 from 0.999 to 0.001. We the same procedure as in the training process to obtain a discrete approximation of the learned latent prior for estimating the score along the reverse-time path.
More numerical results. We performed an additional synthetic experiment in which the nontrivial structure in the first two dimensions is discretely supported:
Example (discrete): The first two coordinates of each 0g are drawn uniformly from the discrete ensembles of three points {(0, 0), (0, 6), (6, 0)).
The noises are then drawn with the same procedure as described earlier to produce the observations.
dimensions dimensions d (8, 16) and traniing sample 2 x 10", SDGM-C stays effective in estimating the discrete locations as d grows. Recent NPMLE solver in Zhang et al. 2022 demonstrates signilicanlly degraded performance such
Figure S: Projected plots of the discrete experiment
We also performed additional circle experiments with d further increases.
Figure 6; Projected plots of the circle experiment onto the first two dimensions in higher dimensions (21, and training sample size
n = 10". SDGM-C stays relatively stable in recovering underlying circular structure d grows to dimension for which numerical NPMLE salvers struggle.
37


Reviews:
[{"review_id": "oaLJB1X94y", "reviewer": "Reviewer_9mdT", "summary": "This paper presents a novel framework that integrates score-based diffusion models with G-modeling for empirical Bayes estimation, particularly targeting heteroscedastic Gaussian mixture models. The authors reinterpret each noisy observation as a terminal state of a diffusion process, enabling the estimation of a nonparametric prior distribution $G^*$ without assuming parametric forms.", "questions": "I have some questions regarding Theorem 1:\n\n1. The requirement to target distribution is only 'compact support'. If we have better property, like $G^*$ has strongly-convex score, how will the bound in Theorem improve?\n\n2. The noise is supposed to be lower bounded by $\\underline{\\sigma}$. Do you require $\\underline{\\sigma}$ to be strictly positive? Or how the result depends on the condition number of $\\Sigma_i$? Also, can you verify these assumptions by some examples in section 5?", "limitations": "Yes.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "The paper delivers a technically sound and mathematically sophisticated contribution to both diffusion modeling and empirical Bayes literature.The use of score matching in G-modeling is novel and theoretically grounded, with convergence results rigorously derived. As far as I know, this is the first work to extend score-based methods to a heteroscedastic Gaussian mixture setting with provable statistical guarantees.\n\nThe only weakness I see is that the main theoretical result (Theorem 1 in Section 3) is not fully interpretable in its current form: the constant dependencies (on $\\sigma$) in the bound are not made explicit, making it difficult to assess practical implications or generality. Also, the assumptions are not verified or illustrated by the experiments. See details below in the Question part.\n\nI would raise my score if the authors are able to (1) clarify the structure and implications of Theorem 1 with precise dependencies and (2) confirm whether the assumptions are satisfied in the presented synthetic and real-data examples.", "was_revised": true, "final_justification": "All of my concerns have been addressed in the rebuttal and so I raised my score.", "rebuttal": "Dear Reviewer **9mdT**,\n\nWe are grateful for your careful recognition of our framework. Indeed, the novel integration of the diffusion reinterpretation with the empirical Bayes G-modeling naturally permits both efficient score-based learning of latent measure and the multi-step Tweedie-type denoising procedure as coherent ingredients. A classical score-based approach with direct parametrization of the score function cannot immediately achieve the two challenges in empirical Bayes we focused on. Our results demonstrated that the proposed framework tackled to a satisfactory extent the over-shrinkage bias and scalability challenges of traditional likelihood-based methods.\n\nWe also want to thank reviewers for bringing constructive feedback based upon which we would like to further improve our manuscript. For each question, we intend to both clarify our understanding, and discuss plans to update our manuscript if we have the opportunity to submit a final version. Due to space limits, we may need to refer you to our responses to other reviewers’ related questions. We apologize for such inconvenience and we are happy to follow up with any persisting issue in the next review phase. \n\nAn overview of our revision plan based on all reviewers’ feedback is as follows: \n* We will better clarify the assumptions, parameters, and implications of theorems in Section 4.\n* We will present a full discussion of limitations and future directions in a new Section 6.\n* We will add numerical results that better explains the method’s scalability and the interplay between $n$ and $d$.\n* We will add explanatory sentences as promised in our rebuttals to reviewers, and be sure to clear all typos noted.\n\nFor papers mentioned below that are cited in [name, year] format, the citations are available in our manuscript’s reference; if instead a paper is cited by number such as [1], we will provide the complementary reference at the bottom. We respond to your questions and concerns as follows.\n\n---\n\n\n***W1. Are assumptions realistic***\n\nThe compact support assumption both covers a vast variety of realistic scenarios and is commonly adopted in prior works. Mixture model (1) is used when researchers believe there appears to be a certain structure to unveil from noisy data. In particular, true signals of our synthetic datasets are constructed to reside in low-dimensional compact manifolds (e.g., a 2D circle or letter “A”) embedded in a high-dimensional space.  The astronomical data is preprocessed by astronomers and statisticians in [Ratcliffe et al., 2020; Sec. 2.2]. Specifically, they retain only those stars that have abundance features confined within the box $[-1, 1]^{19}$, and treat the stars violating the rule as outliers due to anomalous measurements. This crucial step implies the scientific belief that the features of a \"proper\" star concerned in the downstream study should live in a compactly supported space within the box. We discuss more about the downstream relevance in our responses to @fmp1 (Q1). \nIn short, we can certify that all examples in Section 5 indeed satisfy the compact support assumption. These examples are representative of the common scenarios that the model (1) applies to and researchers are interested in. We will make these statements explicit in our final version.\nWhile a mild generalization may be obtained for measures with decaying tails (see e.g., [Soloff et al., 2025]), the rates are vastly similar and we believe compactness is sufficient for conveying the core message of our work and is also satisfied for experiments we performed. \n\nConditions on noise covariances are ubiquitous, important, and satisfied in our experiments. In various scenarios concerned by the model (1), including scientific measurements, covariance information is available or can be estimated from device-specific calibrated noise models, as we justified this in Line 21. For numerical examples in Section 5, the APOGEE dataset comes with known standard deviation estimates bounded within [0.01, 0.62] and has mean 0.05. Our synthetic examples are also constructed with random bounded variances, as shown in Line 316.  Thus, our examples all satisfy the assumption on the noise covariances. \n\nWe now discuss why it is important to assume well-conditioned covariances. Intuitively, we require an upper bound because if $\\overline{\\sigma} \\to \\infty$, then it is information-theoretically impossible to extract useful knowledge from the data. On the other hand, a lower bound on $\\underline{\\sigma} > 0$ is needed to ensure the scores are well-defined everywhere across the space. Since we assume no parametric form of $G^\\*$, it can potentially take a discrete form. Without any perturbation it is nondifferentiable at its component locations. To guard against these malicious possibilities, a positive yet finite noise scale is expected to smooth out degenerate locations and ensure well-behaved scores along an arbitrary SDE path for any $t > 0$. \n\nFurthermore, both the compact support and noise assumptions are commonly used in empirical Bayes literature [Saha and Guntuboyina 2020, Soloff et al. 2025] that proved convergence rates under these conditions. We did not invent any new conditions that are special or unrealistic. We presented these assumptions to enable comparison with existing established results. \n\nRegarding explicitness of constants, we would like to first overview the implicit part that involves $M, d, (\\overline{\\sigma}, \\underline{\\sigma}), g$. As we can imagine, the diffusion coefficient $g$ serves to rescale the covariance and thus can be absorbed in the $\\sigma$ terms. Then we have roughly $C_{M, d, (\\overline{\\sigma}, \\underline{\\sigma})} \\approx C\\cdot d(M^2\\overline{\\sigma}/\\underline{\\sigma})^{d/2}$, which involves scaling due to condition number of covariances $(\\overline{\\sigma}/\\underline{\\sigma})$, and $M^d$ corresponds to covering the high-density region. Our supplemental contains relevant derivations in detail.\nWe also want to emphasize that the near-parametric sample complexity part $\\text{polylog}(n)/n$ is where most literature compares against and focuses on improving, while the leading constant dependencies are usually problem-specific and appear regardless of which learning method is adopted. Our risk presentation style is consistent with both empirical Bayes literature such as [Soloff et al. 2025] and various papers in this venue to highlight and convey the core takeaways. Sample complexity provides the core practical guidance on efficient utilization of collected data. We further discussed the sample complexity in our responses to @DsAZ (W2).\n\n\n\n***Q1. Stronger latent assumptions***\n\nIndeed, if there exists better properties to leverage, it is possible to derive improved bounds. We interpret that your mentioned idea of “$G^\\*$ has strongly convex score”  means that its density $G^\\*(x) \\propto \\exp(-g(x))$ for some strongly convex function $g$. Notice that the score is a vector and it is ambiguous to define a “strongly convex” score. The score-based approach already leverages such convexity properties of Gaussian smoothing that allows efficient sampling and denoising due to good isoperimetric constants (see e.g., [1]). Unfortunately, we cannot think of any obvious way to leverage this additional structural information if the true latent measure is also strongly log concave.\n\nWe do present some conjecture if we can improve the rate with the knowledge of intrinsic low-dimensional manifold dimension in our responses to @DsAZ (W2), which is supported by our experiments and the common belief that diffusion models can detect manifolds.\nOnce again, we want to make assumptions consistent with existing empirical Bayes literature to deliver comparable results. We look forward to developing more advanced theories in future work. We will be sure to envision such directions in the new Limitations section.\n\n***Q2. Covariance assumptions***\n\nWe mentioned above that yes, we do require $\\underline{\\sigma}$ strictly positive for well-defined scores, and the dependence on the condition number. We have also certified the assumptions for our examples above.\nWe will be sure to clarify $0 \\prec \\underline{\\sigma} I_d \\preceq \\Sigma_i$ in our Assumption 1, and discuss the dependency of these problem specifications better after Theorem 1.\n\nIn summary, we appreciate the reviewer to require further clarification on assumptions and realistic implications of our theorem. Our assumptions are pretty mild, commonly used in empirical Bayes, and are satisfied for the examples in Section 5. In addition, statistical literature including prior work in empirical Bayes chooses to keep dependency on these constants less explicit, and pays more attention to the sample complexity as we increase the sample size $n$.  We plan to carefully integrate more explanations and emphasize suitability of these assumptions to suitable locations in our final version, based on the discussion above. We also expect to extend our work further if we can leverage stronger properties. \nPlease let us know if you have any further concerns and we are happy to address them in the next phase.\n\n---\n[1] Koehler, Frederic, Alexander Heckett, and Andrej Risteski. \"Statistical efficiency of score matching: The view from isoperimetry.\" arXiv preprint arXiv:2210.00726 (2022)."}, {"review_id": "dlEH7s6lD1", "reviewer": "Reviewer_fmp1", "summary": "This paper proposed a non-parametric empirical Bayes approach, which is based on a generalized score-based diffusion framework, for learning multivariate Gaussian mixture models with homoscedastic or heteroscedastic noise. The main idea is to view each observation as temporal slices of a family of stochastic diffusion process, which enables a multi-level score-based objective for learning the latent location distribution. Theoretically, a near-parametric error rate was established for the score estimation error. They also show that the proposed method outperforms the non-parametric maximum likelihood estimator in both density estimation and denoising fidelity.", "questions": "-\tIn the real-world astronomy experiment, what constitutes an astrochemically meaningful structure in the denoised output? Is the goal to recover scientifically interpretable features, or is sharpness alone the main criterion for success? Some clarification on the downstream relevance would help assess the practical value of the method.\n-\tThe constants hidden in Theorems 1 and 2 are stated to depend on the dimension ddd. Could you specify how this dependence scales (e.g., linearly, polynomially, or exponentially)? Also, how does the method behave in higher dimensions? The experiments are limited to $d=19$ -- do you anticipate scalability to much higher dimensions?\n-\tHow is the upper bound $T$ for the diffusion time chosen in practice? Does its value affect the convergence rate or the bias-variance tradeoff in the score estimation? Any practical heuristics or theoretical guidance would be appreciated.", "limitations": "The main limitations are the computation of the marginal density $f$. Additionally, the paper lacks discussion and empirical validation of performance in higher-dimensional settings, limiting insight into the method’s scalability.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\n\n- The use of score matching, rather than the classical NPMLE, to estimate the prior within the empirical Bayes framework via G-modeling is novel and insightful.\n- The introduction of a multi-step denoising procedure based on reverse-time diffusion dynamics offers a principled generalization of the traditional one-step Tweedie estimator, mitigating over-shrinkage bias.\n- The theoretical analysis is rigorous and well-developed, providing near-parametric guarantees for score estimation and extending to the more realistic heteroscedastic noise setting.\n\nWeakness:\n\n- While the overall methodological development is well-presented, the paper lacks a clear explanation of how the marginal density $f_{G,\\Sigma_i}$  computed in practice. Given its central role in both the score matching objective and the denoising procedure, a more explicit description would be helpful for reproducibility and understanding.\n- The rationale for preferring score matching over nonparametric maximum likelihood estimation (NPMLE) is not fully justified. As mentioned around Line 233, the convergence rate achieved by score matching merely matches that of NPMLE. Moreover, it seems plausible that NPMLE could also be applied to perturbed observations across multiple noise levels, as in the multi-step denoising setting. Since score matching still relies on computing the marginal score, it is not immediately clear what practical advantages it offers in this context. A deeper comparison—both theoretical and empirical—would strengthen the case for using score-based methods.\n\nTypo:\n- Line 69: scoring matching  -> score matching", "was_revised": true, "final_justification": "The author has adequately address my concerns, and I maintain my positive evaluation of this paper", "rebuttal": "Dear reviewer **fmp1**,\n\nWe sincerely thank you for your kind words about our novelty and technical strengths. Our idea of integrating diffusion reinterpretation with empirical Bayes G-modeling jointly and naturally permits both efficient score-based latent learning and multi-step Tweedie-type denoising. This cannot be achieved with a classical score-based approach with direct parametrization of the score function. Our results demonstrated that the proposed framework did tackle the over-shrinkage bias and scalability challenges of traditional likelihood-based methods.\n\nWe also want to thank you for bringing constructive suggestions based upon which we would like to further improve our manuscript. For each question, we intend to both clarify our understanding, and discuss plans to update our manuscript if we have the opportunity to submit a final version. Due to space limits, we may need to refer you to our responses to other reviewers’ related questions. We apologize for such inconvenience and we are happy to follow up with any persisting issue in the next review phase. \n\nAn overview of our revision plan based on all reviewers’ feedback is as follows:\n* We will better clarify the assumptions, parameters, and implications of theorems in Section 4.\n* We will present a full discussion of limitations and future directions in a new Section 6.\n* We will add numerical results that better explains the method’s scalability and the interplay between $n$ and $d$.\n* We will add explanatory sentences as promised in our rebuttals to reviewers, and be sure to clear all typos noted.\n\nFor papers mentioned below that are cited in [name, year] format, the citations are available in our manuscript’s reference. We respond to your questions and concerns as follows.\n\n---\n\n***W1. Computation of marginal density***\n\nThe computation of the marginal density and its log-gradient is straightforward in practice. If $G$ is represented by a discrete measure with finite number of supports $\\sum_{j=1}^m w_j \\mu_j$ (strategy one in Section 4), then $f_{G, \\Sigma_i}(x) = \\sum_{j=1}^N w_j \\varphi_{\\Sigma_i}(x - \\mu_j)$ where $\\phi_{\\Sigma}(\\cdot)$ is the standard zero-mean multivariate Gaussian density function with covariance $\\Sigma$ that we gave in Line 38, and $f_{G,\\Sigma_i}$ is basically sum of exponentially transformed matrix products. Similarly we can obtain $\\nabla_x f_{G, \\Sigma_i}(x) =  -\\sum_{j=1}^N w_j \\varphi_{\\Sigma_i}(x - \\mu_j) \\Sigma_i^{-1} (x - \\mu_j)$ and the score $\\nabla \\log f_{G,\\Sigma_i} = \\nabla_x f_{G, \\Sigma_i} /  f_{G, \\Sigma_i}$ via chain rule, which is a weighted average of vectors $\\Sigma_i^{-1}\\{x - \\mu_j\\}$. \nAs we explained in Section 4, if $\\widehat{G}$ is represented by a continuous normalizing flow, Monte Carlo sampling can serve as a theoretically backed tool to approximate the density, and the computation can be done with the same formulas above. \nNote that the computations involve mostly matrix operations and can be efficiently performed on GPU with proper tensorization. Our code in supplementary presents how these steps are actually performed and allows reproducibility.\nWe plan to explicitly display the above formulas in Sec. 4 to assist understanding.\n\n***W2. score matching vs. NPMLE***\n\nOverall, we believe that either score matching or NPMLE contributes as a component of a full empirical Bayes framework. We believe that it is our overall framework from a stochastic diffusion perspective that stands out, with the score matching as a coherent and natural ingredient. We are also deeply interested in the possibility to improve solution methods of NPMLE by leveraging established knowledge from literature and our current manuscript. \n\nFor our defense on theoretical convergence rate, please refer to our responses to @DsAZ (W2) where we jointly answered some concerns over how we interpret the rate as a function of $d$. In short, our theoretical analysis aims at achieving the rate of convergence under comparable assumptions as existing work. The same order of dependence on $d$ reveals the fundamental difficulty of the task. Empirically, we observe improved scalability compared to SOTA solution methods for NPMLE (as visible in Figure 2). We believe in possibilities to derive sharper theoretical results if we leverage more structural information.\n\nRegarding the multi-level perturbation for NPMLE, while it seems plausible, it is not clear how NPMLE can benefit from such perturbation. It does not seem to solve either important issue faced by likelihood-based approaches we detailed in Introduction. \nOn the other hand, perturbation is very natural for score learning, and is justified theoretically and proven useful in diffusion literature [Vincent 2011, Song et al. 2021a]. Intuitively, score matching aims at learning a vector field over space thus perturbation is beneficial for exploration spanning the whole space.  \n\nFurther, whether the multi-step denoising is successful heavily relies on accurate estimation of the latent measure. This is because the reverse-time ODE effectively defines a “transport” mapping from observation to the latent measure. However, as our numerical results demonstrate (Fig. 2), existing NPMLE solvers are highly unscalable and the estimated $\\widehat{G}$ is very low-quality even in $d=4$. On the other hand, our numerical results suggest that score-based approach indeed achieves practical advantages in both numeric metrics and visible denoising effects. \n\nWe plan to dedicate a paragraph in the Limitations and Future Work section that summarizes the above discussion and compares the two approaches in detail. We hope our work not only presented a new method for empirical Bayes, but also shed some light on the advances in likelihood-based approaches.\n\n\n***Q1. Interpretation of astronomy experiment***\n\nOne particularly relevant downstream application appears in the cited work [Ratcliffe et al., 2020]. The goal there is to cluster the stars by their chemical similarity using the same 19-dimensional APOGEE chemical abundance dataset, and then examine the resulting groups in age and galactic position to trace disk assembly over time. Before clustering, they quantify pairwise correlations and apply PCA/Isomap to identify influential features and prevent overfitting in high dimensions, yet one difficulty they mentioned explicitly (Sec. 4.2) in the paper is that large uncertainties (i.e., noise we aim to reduce) in the measurements reveal extremely weak correlations among the elements. On the other hand, it is assumed that “our data lie on an embedded nonlinear manifold within the 19-dimensional space” (Sec. 3.1.3). \n\nThese points suggest that a denoising step that reduces measurement uncertainties and recovers the manifold geometry can improve the estimation of true abundance correlations and help downstream clustering mitigate over‑fitting due to noise and strong feature correlations. Consequently, we believe that while sharpness may not necessarily be a criterion for success, the visibly thin, coherent manifolds recovered from the denoised data in our Figure 4 suggests a success, following the above astronomers’ assumptions. Such qualitative criterion is used widely in cited NPMLE works that worked on this data [Soloff et al., 2025, Zhang et al., 2024]. Notice that a clear manifold structure can be successfully recovered in prior work only if they analyze a 2-dimensional projection (e.g., Figure 6 of [Zhang et al., 2024]), while our work successfully applies over the original 19 dimensional space. \n\nFor our final version, we plan to add the following sentence after Line 333 to better demonstrate the practical value of our method:\n“In particular, successful denoising and recovery of low‑dimensional manifold structure in chemical abundance space are crucial for downstream tasks such as clustering stars by mitigating over‑fitting driven by measurement noise and strong inter‑element correlations [Ratcliffe et al., 2020].”\n\n\n***Q2. Scalability in high dimensional settings***\n\nWe once again refer you to our responses to @DsAZ (W2) and @9mdT (W1) for discussion about theoretical scalability. For computational scalability, our method is proved to be much better compared to SOTA NPMLE solvers. Recent paper [Zhang et al. 2024; Fig. 14] (our baselines) reported their struggles even in computing a nontrivial solution for $d\\geq 10$. Our results (e.g. Fig. 2) confirms visible recovery and denoising in $d=16$. The additional experiment in our response to @DsAZ (W2) scales to $d=32$. Thus, we believe in a positive answer to your question.\n\n***Q3. Choice of upper bound $T$***\n\nRecall that $T$ controls how further we wish to smooth the empirical density in the training objective. We did require an upper bound on $T$ to derive the desired guarantees in Theorem 1 & 2, which was shown in the Appendix. Your intuition is right: the core role $T$ plays is indeed related to the convergence rate of the smoothed empirical density from its population counterpart.  It in turn affects the score estimation error. Technically, we pick an upper bound so as to ensure balance among error terms. Intuitively, it follows the same philosophy as in kernel density estimation where stronger smoothing ($T$ larger) leads to both empirical and population distributions collapsed to flatter densities thus the variance reduces, yet we lose more information over the high-density regions leading to high bias. In our final version, we will specify the upper bound more explicitly in the main text. \n\nIn our experiments, guided by the above theoretical intuition, we heuristically found that $T$ need not be set super large for the training to work successfully so we fix $T$ at 1.5 for our experiments. Such hand design is ubiquitous in diffusion literature, e.g., [Song et al., 2021a].\n\nWe would love to further discuss any follow-up questions during the incoming phase."}, {"review_id": "PnqD0ZiNwW", "reviewer": "Reviewer_DsAZ", "summary": "The paper propose Score-based Diffusion with G-Modeling for non-parametric empirical Bayes in multivariate heteroscedastic Gaussian‐location mixtures. The key idea is to reinterpret each noisy observation as the time-1 snapshot of a forward SDE whose latent initial state is drawn from an unknown prior G. This perspective yields a reverse-time ODE (Eq. 7/9) that supplies a multi-step denoising rule generalising Tweedie’s one-shot formula. For the homoscedastic case the authors prove that the plug-in estimator achieves near-parametric Fisher-risk and Hellinger error O(polylog(n)/n) and Wasserstein error O(1/logn) (Thm 1). They extend the guarantees to heteroscedastic noise (Thm 2). Two practical solvers are described: a discrete-support EM-style optimiser and a continuous normalising-flow parameterisation. Experiments on synthetic and real-world data show improvement over NPMLE baselines.", "questions": "See weaknesses above for points about the paper which I am concerned about. I would appreciate clarifications on them to increase my score.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths\n1. They paper proposes a reverse-time denoising procedure and score-matching objective that yield near-parametric Fisher-risk and Hellinger error rates, improving on prior diffusion-model analyses as well as NPMLE estimators. \n\n2. Practical denoising advantage – Multi-step reverse denoiser corrects the classic Tweedie over-shrinkage and demonstrates tangible MSE gains in experiments.\n\n3. Empirical results – On synthetic and real data, SDGM achieves lower Fisher divergence and negative log-likelihood than NPMLE and PEM baselines, supplementing their theory results.\n\nWeaknesses\n1. Requires compact support, well conditioned heteroscedastic covariances, but it is not clear why the assumption holds in practice.\n\n2. Error bounds carry $log^d(n)$ factors; the actual exponent is unspecified, so guarantees may become impractical for large $d$.\n\n3.Each gradient step touches all n observations (or all obsevations in the batch); the paper provides no wall-clock or memory benchmarks, leaving practical scalability unclear. In context of Diffusion models, there was a recent [paper ](https://arxiv.org/pdf/2310.12395) which tried to come up with a non-parametric form of the score function, but they faced certain challenges when compared with neural networks regarding scalability and overfitting. I am curious to see if the rich representation power of neural networks, that has been applied very successfully with diffusion models for score representation, could be applied here as well.\n\n4. Experiments do not compare against recent diffusion-based mixture learners mentioned in the manuscript. It would help to state why these methods don't apply to the authors' experiments directly.", "was_revised": true, "final_justification": "I trust the authors to include the explanations about the dependence on the ambient dimension and their rationale for not comparing with other diffusion-based mixture learners as baselines since their problem setting is different. I would therefore like to maintain my original score.", "rebuttal": "Dear reviewer **DsAZ**,\n\nWe sincerely thank you for recognizing the technical strengths of our proposal. Indeed, the novel integration of the diffusion reinterpretation with the empirical Bayes G-modeling naturally permits both efficient score-based learning of latent measure and the multi-step Tweedie-type denoising procedure as coherent ingredients. A classical score-based approach with direct parametrization of the score function cannot immediately resolve the two challenges in empirical Bayes we focused on. Our results demonstrated that the proposed framework tackled to a satisfactory extent the over-shrinkage bias and scalability challenges of traditional likelihood-based methods.\n\nWe also want to thank reviewers for bringing constructive feedback based upon which we would like to further improve our manuscript. For each question, we intend to both clarify our understanding, and discuss plans to update our manuscript if we have the opportunity to submit a final version. Due to space limits, we may need to refer you to our responses to other reviewers’ related questions. We apologize for such inconvenience and we are happy to follow up with any persisting issue in the next review phase. \n\nAn overview of our revision plan based on all reviewers’ feedback is as follows: \n* We will better clarify the assumptions, parameters, and implications of theorems in Section 4.\n* We will present a full discussion of limitations and future directions in a new Section 6.\n* We will add numerical results that better explains the method’s scalability and the interplay between $n$ and $d$.\n* We will add explanatory sentences as promised in our rebuttals to reviewers, and be sure to clear all typos noted.\n\nFor papers mentioned below that are cited in [name, year] format, the citations are available in our manuscript’s reference; if instead a paper is cited by number such as [1], we will provide the complementary reference at the bottom. We respond to your questions and concerns as follows.\n\n---\n\n***W1. Whether assumptions are realistic***\n \nDue to space constraints, we kindly refer you to our responses to @9mdT (W1). If you find it unsatisfactory for your question, we are happy to follow up in the next phase.\n\n***W2. Dependency on the ambient dimension $d$***\n\nWe presented $O(d)$ in the exponent of $\\log n$ with the aim to reduce burden to understand the core takeaways of the theorem. Our technical supplemental provides the exact exponent as $(5d+2)/2$. This matches previous work with a similar exponent order (Corollary 3.2 of [Saha and Guntuboyina, 2020] gives 2d or 3d under various comparable assumptions). \n\nWe agree with the reviewer’s concerns that the bound $(\\log n)^{O(d)}/n$ becomes impractical for very large $d$. This on one hand reveals the fundamental curse of dimensionality difficulty of the task: We estimate a $d$-dimensional measure and denoise each $\\theta_i$ from a single observation $X_i$. In fact, Kim and Guntuboyina [1] show the minimax rate in Hellinger distance for estimating Gaussian location mixtures is bounded from below by $(\\log n)^d/n$, confirming the difficulty of the task.\n\nOn the other hand, our theorems attempted to demonstrate a fair comparison under comparable assumptions as these work. Such rate is not available in literature before our work for score-based empirical Bayes approach in high dimensional setting.\nAs we mentioned in Related Work (Line 131), whether a non-NPMLE method can achieve a comparable rate with leading $O(1/n)$  factor is only known recently in Ghosh et al. [2025], and we are the first work that extends their $d=1$ setting to arbitrary dimension $d\\geq 1$. In addition, establishing our results is nontrivial, requiring carefully adapting, reinventing, and integrating results from various related fields. Consequently, we believe that our results are significant, which presents a baseline guarantee and, together with our numerical findings, further arouses interest if we could attain better rates if further structural information is leveraged. \n\nIndeed, our experiment results (Figure 2) suggest that a score-based approach is more robust in computational scalability than solution methods of NPMLE in higher dimensions. Based on our experiments and recent belief and observations in literature that diffusion models can identify data manifolds (e.g., in [2]), an educated conjecture is that, if we indeed know the intrinsic dimension $d_0$ of the low-dimensional manifold, instead of the ambient dimension $d$, then we may be able to derive sharper results that avoid the dependence on $d$ in the exponent of the logarithmic factor.\nFollowing reviewers’ suggestions, we conducted additional experiments to investigate two questions: (1) if our method exhibits consistent scalability in higher dimension, and (2) if the empirical interplay between $n$ and $d$, with a particular focus on the latter, performs at least better than the theoretical convergence rate.  We constructed the same synthetic circle dataset appeared in our Sec. 5 with a fixed $n=10^7$ and varied $d = 8, 16, 32$. We applied our SDGM-C, running 5 trials of 6 full data passes for each setting. We report two metrics from our paper: the Fisher divergence (FD) and the negative log-likelihood (NLL), with their normalized values to facilitate comparison across $d$. \n\n| metric \\ d | 8 | 16 | 32 |\n|----------|----------|----------|----------|\n|  NLL   |   13.525   |  24.824    | 47.456  |\n| (normalized)   |  1 |    1.84   | 3.51|\n| |\n| FD   |   0.018   |   0.024    |  0.036 |\n| (normalized)    | 1  |  1.36    |  2.05 |\n\nThe results show both metrics scales (mostly) linear with $d$ in contrast to the unfavorable $(\\log n)^d$ suggested by the baseline theory.  It may motivate the hypothesis of an improved rate such as $C d{(\\log n)^{d_0}}/{n}$. \n\nFor theoretical justification on improving risk rates, we look forward to better leveraging the structure and exploring various different proof tools as a future work.\nFor the current manuscript, based on the above discussions, we plan to (1) revise the theorem implication paragraph (Lines 222 to 242) and add the higher dimensional experiments with varying $(n, d)$ in Section 5 to better demonstrate and interpret the scalability issue, and (2) fully explain the theoretical limitations and envision the future directions in the added section. \n\n***W3. Resource benchmark and neural network***\n\nWe appreciate that the reviewer shares the insightful paper with us. It looks like their goal is to circumvent training, while our method, though nonparametric, did perform parametrization of $G$ and requires training. In essence, the term “nonparametric” as in various literature means that we do not assume any parametric form of the true target (i.e., $G^\\*$ is not assumed to belong to any specific parametrized measure class), but we can use various parametrized models to serve the learning task.\nIndeed, our implementation already makes use of deep neural networks and benefits from GPU speedup. As we explained in Section 4, the normalizing flows are a class of neural networks that we used to parametrize $G^*$ continuously and our SDGM-C tends to perform better than the discrete counterpart. \n\nFor a quick glimpse into the resource benchmarks, we report  average wall-clock time per step and the peak GPU memory usage over the training tasks with batch size 1024 and 500k Pytorch model parameters performed on 2x Nvidia V100 GPUs:\n| d  | Time (second per step) | Peak GPU usage (GB) |\n|----------|----------|----------|\n| 16    |    0.53      |    2.0     |\n| 32    |     1.17     |     5.5     |\n\nFor a very rough comparison, the statistics reported in the official repo for [Song et al., 2021a] shows their Pytorch model runs 0.56 second per step and uses 20.6 GB memory in total on 4x V100 GPUs. \nWe believe there exists room for optimizing our implementations and we expect to explore these practical optimization considerations as a future work. \n\n***W4. Lacking other diffusion-based mixture learners as baselines***\n\nWe did not compare with diffusion-based mixture learners [Shah et al. 2023, Kim and Ye 2022] primarily because of two reasons.\n* Resultant denoisers from these methods can only be attained through the one-step Tweedie and the unfavorable over-shrinkage bias is unavoidable. \n These methods are not compatible with the multi-step denoising procedure, due to they all directly parametrize the score function without G-modeling. As we mentioned in Line 180, without G-modeling,  scores are not defined between $t = 0$ to 1 and are not learnable because no data exists in the time interval. It is implausible for the score network learned in forward time ($t > 1$) to interpolate for the backward time ($t < 1$), thus the multi-step denoising steps would not work. \n* All these learners do not learn the latent distribution. By learning the scores directly, they only permit efficient generation of samples from the noisy mixture, yet no $\\widehat{G}$ can be recovered. In such cases,  half of our metrics (NLL and W2) are uncomputable and cannot be reported for a fair comparison.\n \nGiven the above reasons, we believed that while relevant to our work, they are too brute-force for serving the core task, which has evenly weighted focuses on both the recovery of latent measure and on denoising the observations. These approaches didn’t excel in either task. Consequently, we didn’t include them as a valid baseline. If the reviewer finds some strong reasons to include them, we are happy to further discuss them.\n\n---\n[1] Arlene K. H. Kim, Adityanand Guntuboyina \"Minimax bounds for estimating multivariate Gaussian location mixtures,\" Electronic Journal of Statistics, Electron. J. Statist. 16(1), 1461-1484, (2022)\n[2] Pidstrigach, Jakiw. \"Score-based generative models detect manifolds.\" Advances in Neural Information Processing Systems 35 (2022): 35852-35865."}, {"review_id": "u3mdEMBb87", "reviewer": "Reviewer_KpwC", "summary": "The paper proposes a method to integrate nonparametric empirical Bayes with score-based diffusion generative  models. The goal is to estimate an unknown prior distribution from noisy observations. In doing so, the paper addresses both statistical and computational limitations of the classical approach, which combines the nonparametric maximum likelihood estimator (NPMLE) and Tweedie-based posterior denoising. The proposed method generalizes and improves upon existing approaches, which either do not account for heteroscedastic measurement noise or lack theoretical guarantees in dimensions higher than $1$. \n\nThe key idea is to reinterpret each noisy observation as a time-1 snapshot from a latent stochastic diffusion process. This perspective naturally induces a denoising score matching objective for explicitly learning the prior distribution $G^*$. The resulting optimization problem is infinite-dimensional and nonconvex. The posterior mean is then estimated by combining this approach with a multi-step Tweedie-type empirical Bayes denoising procedure. \n\nAssuming access to the minimizer of the non-convex optimization problem, the paper provides theoretical guarantees on the statistical performance of the method. Finally, numerical experiments on synthetic and real-world datasets are presented to validate the theory, with particular attention paird to robustness as the dimension of the problem varies.", "questions": "- I appreciate a lot how the paper is organized. However, I wonder why the authors do not emphasize more strongly the result regarding the interplay between $n$ and $d$ in the presentation, since dimension plays a role also in the experiments. \n\n- I am confused by the example about the over-shrinkage bias. Isn't the reduced variance simply a consequence of observing $X$, which provides information and hence reduces uncertainty in the posterior estimation? \n\n- How difficult would it be to quantitatively incorporate the optimization error of the nonconvex estimator into the statistical analysis of the method?", "limitations": "See Weaknesses.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strenght:** \n\nThe paper is well-written and provides both a theoretical treatment and numerical experiments that illustrate the theory.  One interesting feature of the analysis is that the error bounds highlight the interplay between the number of observations and the dimensionality of the problem – something that previous literature didn’t address in this context. Finally, some assumptions are needed to make the analysis clean; these are carefully discussed with reference to prior literature.\n\n\n**Weakness:** \n\nTo my understanding, the paper does not provide strong algorithmic guarantees. To me, this is not a major issue, since the authors are transparent about their assumptions; however, these limitations should be discussed explicitly in a “Limitations” section, which is currently missing. For example, assuming solvability of the nonconvex problem is, I believe (though I am happy to change my mind based on the authors’ response), quite strong, as the optimization can easily get stuck in a local minimum. It would also be interesting to incorporate in the error analysis of Theorem 1-2 the fact that $\\widehat{G}$ may not be a good  approximation of the minimizer of the nonconvex objective. \n\nI also noticed a typo on line 182: a repeated “that”.", "was_revised": true, "final_justification": "As of writing:\n\n- The paper is well-written and provides both a theoretical treatment and numerical experiments that illustrate the theory. \n- The authors promised to extend the discussion on the interplay between $n$ and $d$, which I think is a relevant aspect of the paper.\n\nI am still unsure about the assumption concerning the solvability of the infinite-dimensional nonconvex problem. However, after reading the author's rebuttal and the other reviews, I realize that this concern might not represent a major obstacle.", "rebuttal": "Dear reviewer **KpwC**,\n\nWe truly appreciate your recognition of our organization and contributions. Indeed, it is the diffusion reinterpretation integrated with the empirical Bayes G-modeling that jointly and naturally permits both efficient score-based learning of latent measure and the multi-step Tweedie-type denoising procedure. This is not achievable with a classical score-based approach with direct parametrization of the score function. Our results did showcase that we tackled to a satisfactory extent the over-shrinkage bias and scalability challenges of traditional likelihood-based methods.\n\nWe also want to thank reviewers for bringing excellent suggestions based upon which we would like to further improve our manuscript. For each question, we intend to both clarify our understanding, and discuss plans to update our manuscript if we have the opportunity to submit a final version. Due to space limits, we may need to refer you to our responses to other reviewers’ related questions. We apologize for such inconvenience and we are happy to follow up with any persisting issue in the next review phase. \n\nAn overview of our revision plan based on all reviewers’ feedback is as follows:\n* We will better clarify the assumptions, parameters, and implications of theorems in Section 4.\n* We will present a full discussion of limitations and future directions in a new Section 6.\n* We will add numerical results that better explains the method’s scalability and the interplay between $n$ and $d$.\n* We will add explanatory sentences as promised in our rebuttals to reviewers, and be sure to clear all typos noted.\n\nFor papers mentioned below that are cited in [name, year] format, the citations are available in our manuscript’s reference; if instead a paper is cited by number such as [1], we will provide the complementary reference at the bottom. We now begin our responses to your questions.\n\n---\n\n***W1 & Q3. Algorithmic guarantees***\n\nWe thank the reviewer for this thoughtful comment. The current manuscript aims to propose a new framework for nonparametric empirical Bayes and to analyze its statistical guarantees. We agree that a more thorough discussion of the algorithmic aspects of our approach is valuable and acknowledge that our current manuscript does not fully address this important point. We will explicitly discuss these issues in the limitation section.\n\nIn particular, we agree that assuming solvability of the nonconvex optimization problem is a strong assumption. In the existing literature, global optimality guarantees for training neural networks are often established under overparameterization, for example via the neural tangent kernel regime [1]. In such settings, the objective behaves almost linearly and can satisfy favorable conditions such as the Polyak–Lojasiewicz inequality. Incorporating this line of analysis that linking optimization theory with approximation error of deep neural networks and statistical estimation bounds would provide a more comprehensive understanding of the proposed method. We will include these points in our revised manuscript.\n\n***Q1. Interplay between n and d***\n\nWe kindly refer you to our responses to @DsAZ (W2) where we bear in mind your question in answering there and it should cover most aspects concerning such interplay. In short, we plan to incorporate an additional experiment that explicitly play around these dependencies.  \nWe are happy to further discuss if you find it insufficient to answer your question here.\n\n***Q2. Understanding the over-shrinkage bias***\n\nYour intuition from the information-theoretic perspective is certainly correct. The example demonstrates that the posterior mean (or one-step Tweedie) estimator $\\mathbb{E}[\\theta |X]$ is not sufficient for explaining the full variation in the distribution of $\\theta$. Indeed, using the total law of variance, we have:\n$ Var(\\theta) = Var[\\mathbb{E}[\\theta | X]] + \\mathbb{E}[Var(\\theta | X)]$. This shows that the posterior distribution of random variables $\\mathbb{E}[\\theta |X]$ can only partially explain the actual variation in $\\theta$, and thus leads to over-shrunk estimators for $\\theta_i$. \n\nIn contrary, if we apply the multi-step Tweedie, then on one hand the ODE defines a deterministic flow path that “transport” the observation distribution to the (approximated) latent distribution, thus we can expect the resulting denoised estimator to distribute according to the latent distribution without any reduced variance. On the other hand, related to the above analysis, we may again understand from the information-theoretic perspective. Assume we know the true latent distribution and follow Equation (6) in our paper to denoise, then the denoising sequence of random variables $\\\\{X_i\\\\} = X^{(t=1)}, X^{(1-\\Delta t)}, \\dots, X^{(\\Delta t)}, X^{(0)} = \\\\{\\theta_i\\\\}$ defines a Markov chain connecting latent $\\theta$ and observables $X_i$ and we wish to explain the variability in $\\theta_i$ through the sequence. Then by the data processing inequality, we have $Var[\\mathbb{E}[\\theta | X^{(1)}]] \\leq Var[\\mathbb{E}[\\theta | X^{(1-\\Delta t)}]] \\leq \\cdots \\leq Var[\\mathbb{E}[\\theta | X^{(\\Delta t)}]]$, i..e, we gain enhanced explanatory power as we move along the ODE path. This clearly demonstrates why we prefer the multi-step procedure.\n\nProvided more space in the camera-ready version, we will attempt to incorporate the above analysis into our final manuscript.\n\nWe would be happy to further discuss any questions or concerns that may arise during the discussion phase.\n\n---\n[1] Jacot, Arthur, Franck Gabriel, and Clément Hongler. \"Neural tangent kernel: Convergence and generalization in neural networks.\" Advances in neural information processing systems 31 (2018)."}]

Author rebuttal:

Meta-review: During the rebuttal period the authors addressed several concerns raised by reviewers, who all agree that the paper has important merits and deserves to be published. Hence I recommend an acceptance.
Decision: Accept (poster)