We are grateful for your extremely thorough and constructive review. Your detailed questions have helped us identify several areas for clarification. We address your comments in detail as follows.

W1: Pessimistic rates.

We agree that our convergence rates are conservative, especially for high-dimensional settings. As our work aims to offer worst-case guarantees under minimal assumptions, this is somewhat inevitable. Nonetheless, we believe our results fill an important gap in understanding what is theoretically achievable under weak regularity. We see breaking the CoD through further structural assumptions as an important future direction.

Q1: Sample-dependence of the score network

Thank you for this sharp and important question. You are correct that the score network constructed in our approximation results is sample-dependent. We clarify how this is handled in our analysis below.

1. The approximating network in Theorem 1. The neural network $\phi_{\text{score}}$ in Theorem 1 is indeed a random network. It is constructed to approximate the regularized empirical score function (a KDE-based estimator), and its parameters therefore depend on the specific sample used to construct the kernel estimator. The statement of Theorem 1 should be interpreted as follows: for any given sample $\bar{S} := \bar{x}^{(1)}, \dots, \bar{x}^{(n)} \sim P_0^{\otimes n}$, there exists a network with the specified architecture constructed from $S$ that achieves the stated approximation error. The expectation in the theorem's bound should be understood as an expectation taken over both the test point $x_t$ and the randomness of the sample $\bar{S}$, i.e., $E_{\bar{S}}[E_{x_t \sim P_t}[\cdot]]$. We will revise the statement of Theorem 1 to make this sample dependence explicit.

2. The learned network in Theorem 2. This sample dependence is a critical point that we carefully account for in our generalization analysis for the learned network $\widehat{\phi}$ in Theorem 2. Our proof strategy is as follows:

   • We introduced an independent "ghost sample" $\bar{S}$ to construct the target regularized empirical score function, $\widehat{\bar{\phi}}$ (Line 1689-1690 Page 79), as noted by the reviewer.

   • The "ideal" approximating network $\tilde{\phi}$ (Line 1725 Page 83) from Theorem 1 can then be considered to be the one constructed based on this ghost sample. This makes $\tilde{\phi}$ a random target that is independent of the actual training data.

   • The final network $\widehat{\phi}$ is obtained by performing empirical risk minimization on the actual training data.

   Our proof correctly bounds the excess risk by properly handling the randomness from both sample sets. Therefore, our generalization bound and proof of Theorem 2 remain valid.

We will revise the manuscript to improve clarity on this crucial point.

W2/Q2: Where $\tilde{\mathcal{O}}(n^{3/d})$ depth arises from, and whether it can be improved.

This is an excellent and constructive point. We agree that a network depth that grows polynomially with the sample size is a practical limitation. This can be improved within our framework.

The depth of $\tilde{\mathcal{O}}(n^{3/d})$ (for $k \geq d/2$) arises from a specific choice of parameters $N = 1, L = \lceil n^{1/(2k)} \rceil, k \geq d/2, \epsilon = n^{-1/k}, s=k$ in the proof of Lemma 12 (see Line 1447-1448 Page 67), which are made to satisfy a technical condition from Lemma 11 and 12 (i.e., $N^{-2}L^{-2} \leq t_0 \wedge n^{-1/k}$). However, this choice is not fundamental. Our framework allows for separate control over network width and depth. Following the reviewer's suggestion, we can revise our construction to achieve a more realistic architecture. By swapping the roles of width and depth, we can choose $N = \lceil n^{1/(2k)} \rceil, L = 1$ and construct a network with width $\leq \mathcal{O}(n^{\frac{3}{2k}}\log_2n)$ and depth $\leq \mathcal{O}(\log^2n)$. This new architecture is more practical and, crucially, our main theoretical guarantees remain valid:

• The condition $N^{-2}L^{-2} \leq t_0 \wedge n^{-1/k}$ required by Lemma 11 and 12 is still satisfied.

• The generalization bound in Theorem 2 remains valid since the new network architecture's pseudo-dimension remains of the same order, which is bounded by $\tilde{\mathcal{O}}(n^{\frac{3}{k}})$.

We are grateful for this suggestion and will update the manuscript to reflect this significantly improved architecture.

Q2: Dependence on dimension $d$.

Thank you for the excellent question. This counter-intuitive relationship is a direct consequence of the interplay between the smoothness of $P_{t_0}$ and the network architecture in our theoretical analysis, which has a clear theoretical motivation:

• Higher dimensions force more smoothing. Our condition on the early stopping time (i.e., $t_0 \geq \tilde{\mathcal{O}}(n^{-2/d})$ in Theorem 1 and 2 as well as $t_0 = n^{\frac{-2}{d+2s}}$ in Corollary 1 and 2) forces $t_0$ to be larger for higher dimension $d$ to ensure theoretical guarantees.

• More smoothing makes the marginal distribution $P_{t_0}$ easier to learn. A larger $t_0$ means that $P_{t_0}$ is convolved with a Gaussian of larger variance. This makes $P_{t_0}$ inherently smoother and less complex, as fine-grained details are averaged out. As a result, its score can be approximated by a network whose size scales less severely with the sample size $n$.

In summary, while higher dimension usually implies increased statistical difficulty (as seen in the convergence rate $n^{-s/(d+2s)}$ suffering from the curse of dimensionality), the necessity of stronger smoothing to achieve uniform control in high dimensions makes the intermediate learning problem architecturally less demanding in terms of network size. This is a reflection of the bias–variance trade-off: we reduce variance (by smoothing), but incur bias, which is ultimately what limits the rate. We will clarify this point in the revised manuscript.

Q3: Use of $\rho_n$ in $\mathcal{G}$ Definition.

Yes, this should be written in terms of $\rho_n$. We will correct this typo.

Q4: Extension to sub-exponential or sub-Weibull distributions

This is an excellent suggestion for future work. Our current analysis relies on the sub-Gaussian assumption in two critical places:

• Lemma 7: To bound the score function's norm over low-density regions.

• Lemma 3: To bound the KL-divergence of the Gaussian smoothed empirical distribution, which is fundamental to our main error bounds.

Extending these results to broader tail classes like sub-exponential or sub-Weibull would require more refined control of tail integrals, but we agree it is plausible and an exciting avenue for future work. In particular, we believe that the extensions may benefit from the convergence rates of Gaussian smoothed empirical distributions under heavy-tail distribution assumptions. We will include the discussion in the revision.

Q5: Use of $f^{(1)}, f^{(2)}$.

These notations correspond to different components in the DNN decomposition (e.g., numerator vs. denominator in the score approximation), and are treated separately for modular construction. An alternative approach would be to construct a single sub-network for this expression and reuse its output via connections to different parts of the network. Since both methods result in a final network of a similar size and complexity, our notation was intended to reflect this separation.

Q6: Typo "estiamtors"

Thank you. We will correct the typo.

Q7: Clarification on “Standard” vs. “General” Gaussian-convolved Densities

"Standard Gaussian-convolved" refers to convolution with a standard Gaussian distribution $\mathcal{N}(0, \mathbf{I}_d)$, whereas "general Gaussian-convolved" refers to convolution with a general isotropic Gaussian distribution $\mathcal{N}(0, \sigma^2\mathbf{I}_d)$ for any $\sigma > 0$. We will improve this terminology for clarity.

We sincerely thank you for your thoughtful feedback and acknowledgement of the value of relaxing assumptions. We agree that many of the tools we employ, such as Bernstein’s inequality and Girsanov’s theorem, are well-established. Our primary contribution, however, lies in the novel synthesis and application of these tools to remove critical and restrictive assumptions from prior work, which required new, non-trivial technical developments. We address your comments in detail below.

W1 & W2: Use of standard tools and limited novelty of the contribution

While the individual components are standard, their integration within our framework is what enables our main results. Specifically, our key technical contributions are:

• A refined KL-divergence bound (Lemma 3): Prior work provided a KL-divergence rate for smoothed empirical distributions that scaled exponentially with the inverse of the noise standard deviation $1/\sigma$. As detailed on page 8, we refine their proof to establish a rate that scales polynomially in $1/\sigma$. This improvement is crucial for deriving meaningful bounds in the low-noise regime.

• New application of Bernstein's inequality: To remove the density lower-bound assumption from [22, 26, 27], we introduce a new proof strategy. Instead of analyzing the true score function, which may be unbounded, we analyze the excess risk for learning a surrogate regularized empirical score function. This surrogate is uniformly bounded by construction (Lemma 5), allowing us to verify Bernstein's condition and obtain a uniform high-probability bound without requiring a density lower bound.

The combination of these non-trivial technical steps that allows us to provide the a theoretical framework establishing nearly minimax optimal rates for neural network-based SGMs under the general sub-Gaussian assumption, without the restrictive requirements of Lipschitz continuity or a density lower bound. We believe this significantly expands the applicability of SGM theory and represents a valuable contribution to the literature.

We sincerely thank you for this valuable feedback. Your observation is correct and highlights a key feature of our framework: it accurately reflects the fundamental trade-offs in nonparametric estimation. The behavior you noted is expected and does not undermine our main results. We address your concern in detail as follows.

W1/Q1: Divergence of bounds with sample size $n$

While the reviewer’s observation is correct, this behavior is expected within our framework and does not undermine our main results. Specifically,

• It is expected with only a sub-Gaussian distribution assumption. Under only the sub-Gaussian assumption, we establish a convergence rate of $\tilde{\mathcal{O}}(t_0^{-d/2} n^{-1})$ for any $t_0 \geq 1/2\alpha^2n^{-2/d} \log n$. The rate indeed becomes vacuous (or diverges) if $t_0$ is taken too close to this lower limit, which is expected and consistent with the findings in [29]. In fact, [29] also shows that without assuming smoothness, the kernel density estimator cannot achieve meaningful convergence rates either. Our framework correctly captures this. It reflects the fundamental difficulty of estimating an arbitrarily "rough" distribution.

• The bound converges to zero if we select an early stopping time $t_0$ that decreases with $n$ slightly slower than $1/2\alpha^2n^{-2/d} \log n$. Noticing that we only require the early stopping time $t_0$ to be in a admissible range $t_0 \geq 1/2\alpha^2n^{-2/d} \log n$ and the bound converges to zero when the early stopping time $t_0$ is set slightly larger, for instance, $t_0 = \mathcal{O}(n^{-2/(d+\epsilon)})$ for any $\epsilon>0$. The convergence rate would be of $\tilde{\mathcal{O}}(n^{-\frac{\epsilon}{d+\epsilon}})$.

• $t_0$ can be tuned to achieve nearly-optimal rates with additional smoothness assumptions. In particular, as detailed in Corollaries 1 and 2, when we add standard nonparametric smoothness assumptions, we can tune the early stopping time $t_0$ to balance the bias $TV(P_0, P_{t_0})$ and statistical error $TV(P_{t_0}, \widehat{P}_{t_0})$. For instance, with a Sobolev density, we choose $t_0 = n^{-2/(d+2s)}$, which satisfies the admissible range $t_0 \geq 1/2\alpha^2n^{-2/d}\log n$, and we recover the nearly optimal convergence rate of $\tilde{\mathcal{O}}(n^{-2s/(d+2s)})$.

So in summary, our framework provides a general bound valid for all $t_0$ above the threshold, and when smoothness is available, the early stopping can be tuned to achieve the nearly minimax optimal rate.

We will clarify in both the abstract and main claims that minimax-optimality holds for Sobolev/Besov or otherwise regularized classes, while for the pure sub-Gaussian case we provide the convergence guarantee for $P_{t_0}$ when $t_0 > \mathcal{O}(\alpha^2n^{-2/d}\log n)$.

We sincerely thank you for your insightful comments and constructive feedback, which have helped us improve the quality of our paper. We address your comments in detail below.

W1: More discussions about the technical challenges for tackling the curve of dimensionality; Suggestion to discuss [1*, 2*].

We thank the reviewer for this constructive suggestion. Indeed, our current theoretical guarantees suffer from the curse of dimensionality (CoD), as we explicitly acknowledge in the conclusion. Incorporating structural assumptions, such as manifold assumption explored in [1*, 2*], or low-rank structures) is a crucial next step to mitigate the CoD. This would require non-trivial extensions to our analysis, particularly in adapting our empirical Bayes smoothing techniques and neural network approximation theory to efficiently exploit low-dimensional geometry. In particular, our analysis relies on an isotropic Gaussian kernel for the empirical score estimator. To effectively leverage a manifold structure, this would need to be replaced with an estimator that respects the underlying geometry, such as one using anisotropic or manifold-intrinsic kernels (e.g., heat kernels), to prevent smoothing data off the manifold. A manifold-aware analysis, following, e.g., the path of [1*, 2*], would require establishing new bounds for score estimation and approximation that depend on the intrinsic dimension of the data, rather than the ambient dimension. While this is beyond the scope of the current work, we view it as a promising direction. We thank the reviewer for this direction and will add a substantive discussion in the revision, referring to the cited works for context.

W2: Discussing the tightness of Theorems 1 & 2.

Thank you for pointing this out. While we do not provide a full minimax lower bound for the neural network-based score estimation in Theorems 1 and 2, we establish the near-optimality of our framework in two key aspects:

• Our rates in Theorem 2 match the minimax-optimal score estimation rate for regularized kernel-based estimators (see Lemma 1), up to logarithmic factors. In particular, our results align with the optimal rate derived in [29] for a similar truncated estimator, indicating that this foundational part of our analysis is indeed tight.

• When we supplement our sub-Gaussian assumption with standard nonparametric smoothness assumptions (i.e., Sobolev or Besov classes), our framework demonstrates that neural network-based SGMs can achieve nearly minimax optimal convergence rates for distribution estimation in total variation distance. This is detailed in Corollaries 1 and 2. These results show that our framework is fundamentally sound and capable of achieving optimal rates when appropriate conditions are considered. We will clarify this tightness explicitly in the revised manuscript.

W3: Adding experiments to illustrate the validity of the established bounds.

We appreciate the suggestion. As our work is purely theoretical and aims to answer fundamental questions about SGMs under minimal assumptions, we prioritized analytical clarity. We agree that illustrative experiments would be valuable for validating bounds empirically and will consider adding such results in future versions or follow-up work.

Q1: For the score estimation via empirical Bayes smoothing, what is the technical challenge to remove the Lipschitz score assumption compared with [30]?

The primary challenge lies in controlling the score matching loss without causing the error bounds to depend sub-optimally on the noise variance $\sigma$. In particular, Lemma 2 of [30] established a convergence rate of the smoothed empirical distribution in Hellinger distance $H^2$ under the Lipschitz score assumption. As discussed in [30] following their Lemma 2, one can instead bound the Hellinger distance by KL divergence by the fact $H^2 \leq KL$ without the need of the Lipschitz score assumption. However, the SOTA convergence rate of the smoothed empirical measure in KL-divergence established in existing works, see [42], has a constant scale exponential in $1/\sigma$, i.e., $\mathcal{O}(\exp(1/\sigma)\frac{\log^dn}{n})$. We improve this KL-divergence convergence rate in Lemma 3 for the constant from exponential dependence to polynomial dependence on $1/\sigma$, i.e., $\mathcal{O}(\sigma^{-d}\frac{\log^{d/2}n}{n})$. This crucial improvement allows our bounds to hold even for small noise levels without needing the Lipschitz assumption to control the error. We believe this result on the convergence of smoothed empirical measures may be of independent interest.

[42] Rate of convergence of the smoothed empirical Wasserstein distance. 2022