1 Introduction↩︎

Feedforward neural networks (FNNs) with rectified linear unit (ReLU) activations have received much attention in nonparametric regression due to their ability to efficiently approximate functions with latent low-dimensional structures [1]–[3]. This structural flexibility has facilitated their applications to a broad range of statistical problems, including robust regression [4], [5], survival analysis [6], and factor and interaction models [7], [8], among others.

While ReLU FNN estimators achieve minimax-optimal convergence rates in the \(L^2(P)\) norm [4], their uniform convergence properties remain less understood. Yet, uniform convergence is a fundamental prerequisite for many downstream statistical theoretical analyses and decision-making applications, such as ensuring reliable individualized treatment recommendations [9], constructing valid confidence bands [10], and enabling transfer learning [11]. Current theoretical guarantees for the uniform convergence of ReLU FNNs are notably limited: they are often restricted to shallow ReLU network estimators with univariate covariates [11] or rely on the availability of a uniformly consistent pilot estimator [12]. This critical theoretical gap raises significant concerns regarding the reliability of ReLU FNN estimators in contexts where uniform guarantees are essential, which restricts their practical applicability.

To investigate this open problem, we first characterize a fundamental bottleneck in the uniform convergence of ReLU FNNs. Specifically, we establish the first theoretical lower bound in the literature demonstrating that least-squares ReLU FNN estimators inherently suffer from the curse of dimensionality in their uniform convergence rates, even when they simultaneously achieve the minimax-optimal convergence rate in the \(L^2(P)\) norm. This underscores that the optimal \(L^2(P)\) convergence of ReLU FNNs does not translate to the uniform reliability required for downstream statistical theoretical analyses.

Beyond these theoretical limitations, the focus on ReLU FNNs increasingly diverges from modern deep learning practice. The piecewise linear ReLU activation has largely been replaced by \(C^\infty\)–smooth alternatives, such as the sigmoid linear unit (SiLU) [13], [14], and the Gaussian error linear unit (GELU) [15], driven by their superior empirical performance. Furthermore, residual architectures (ResNets) [16], [17] have fundamentally reshaped neural network design. By resolving the vanishing gradient and degradation problems inherent in training very deep networks, residual connections have become the indispensable backbone of modern deep learning. Although the statistical theory underlying these components remains limited, the integration of residual architectures with \(C^\infty\)-smooth activations forms the structural basis of current powerful foundation models, from Vision Transformers [18] to state-of-the-art large language models like LLaMA [19] and DeepSeek [20].

To overcome this bottleneck, motivated by recent architectural developments in deep learning, we focus on deep neural networks with smooth activations (smooth DNNs), encompassing both feedforward and residual architectures, and develop a rigorous theoretical framework for their uniform convergence analysis. Through a comprehensive theoretical analysis together with supporting simulation studies, we show that smooth DNNs can substantially mitigate the curse of dimensionality in uniform convergence and provide strictly stronger uniform guarantees than ReLU FNNs. These results position smooth DNNs as a principled and theoretically grounded alternative for statistical learning tasks in which uniform convergence is essential.

The primary methodological and theoretical contributions of this work are threefold. Together, they bridge a critical gap between the empirical success of smooth DNNs and the current lack of a general statistical theory for their uniform convergence:

• Characterizing the curse of dimensionality for ReLU FNNs in the uniform convergence. We establish the first theoretical lower bound demonstrating that ReLU FNNs inherently suffer from the curse of dimensionality in the uniform norm. While achieving minimax-optimal \(L^2(P)\) convergence [4], Theorem 1 reveals that the uniform convergence rate of the least-squares ReLU FNN estimator can be bounded below by \(n^{-\frac{1}{d+1}}\) and above by \(n^{-\frac{1}{d+2}}\), where \(n\) is the sample size and \(d\) is the feature dimension, regardless of the Hölder smoothness of the target function. Furthermore, our analysis based on interpolation inequalities in Section 2.2 suggests that this limitation stems fundamentally from the limited smoothness of the ReLU activation, which prevents ReLU FNNs from exploiting higher-order regularity.

• A foundational theoretical framework for smooth DNNs. The statistical theory for deep neural networks, particularly ResNets, with smooth activations, remains largely absent. Existing literature predominantly focuses on either FNNs with less common activations, such as sigmoid, tanh, or rectified power unit (RePU) activations [21]–[23] or residual architectures with ReLU activations [24]–[26]. To enable the uniform convergence analysis of smooth DNNs, Section 4 develops a comprehensive set of theoretical tools. Specifically, we establish an upper bound on the pseudo-dimension in Theorem 3 and derive approximation error bounds for both Sobolev functions and hierarchical composition models in Theorems 4 and 5, respectively. To our knowledge, these results provide the first theoretical foundation for establishing uniform convergence guarantees for smooth DNN estimators.

• Uniform convergence and robustness guarantees. Building on our theoretical framework, we establish uniform convergence guarantees for smooth DNN estimators in Huber, least-squares, quantile, and logistic regression. In particular, Theorem 8 shows that the smooth DNN Huber estimator is non-asymptotically robust in the uniform norm, extending existing robustness results for ReLU FNN Huber estimators in the \(L^2(P)\) norm [4], [5]. Theorem 9 then gives the corresponding uniform convergence rate for the smooth DNN least-squares estimator as a special case of Huber regression, and shows that it can overcome the curse of dimensionality exhibited by ReLU least-squares regression in Theorem 1. Analogous uniform convergence and robustness guarantees for quantile and logistic regression are established in Theorems 10 and 11, respectively. Notably, the uniform convergence guarantee for logistic regression also provides rigorous theoretical support for downstream applications such as probability estimation for double Higgs boson production [27]; see Theorem 12. More broadly, across all these tasks, we prove that smooth DNN estimators can adapt to the low-dimensional hierarchical composition structure of the target function, thereby enjoying a clear theoretical advantage over ReLU FNNs in uniform convergence.

To support our methodological and theoretical results, Section 6 benchmarks the numerical \(L^2(P)\) and uniform estimation errors of smooth DNN estimators against those of ReLU FNN estimators in Huber regression. The substantial performance gains achieved by smooth DNN estimators provide empirical support for our methods and theory. Taken together, these results support smooth DNN estimators as effective alternatives to ReLU FNNs in practical applications where uniform convergence is required.

The rest of this paper is organized as follows. Section 2 establishes the lower and upper uniform convergence rates of the least squares ReLU FNN estimators. Section 3 presents the new methodology for smooth DNNs. In Section 4, we establish the key statistical properties of smooth DNNs. Building on these foundational elements, in Section 5, we derive uniform convergence rates of smooth DNN estimators for Huber, least-squares, quantile, and logistic regression. In Section 6, we present simulation studies and a real application. Section 7 includes a few concluding remarks. The complete proofs and additional technical or numerical results are presented in the Supplementary Materials.

Before proceeding, we introduce the following notation and terminology that will be used throughout this paper. We use the notation \(a\lesssim b\) to mean there is a constant \(C>0\) independent of \(a\) and \(b\), such that \(a\leq C b\). We say \(a\asymp b\) if both \(a\lesssim b\) and \(b\lesssim a\) hold. We say \(a\lesssim_{\log n}b\) if there exists a function \(C(\log n)\) such that \(a\leq C(\log n)\cdot b\). We may further use the notation \(\lesssim_{\log}\) to suppress any other logarithmic terms. \(a\lor b = \max(a,b)\) and \(a\land b = \min(a,b)\). For \(n\in\mathbb{N}_+\), denote \([n]=\{1,2,\cdots,n\}\). For \(x\in \mathbb{R}\), denote \(\lceil x\rceil=\min\{n\in\mathbb{N}: n\geq x\}\). For \(x\in\mathbb{R}^d\), denote its \(\ell_2\) norm as \(\|\cdot\|_2\) and \(\ell_\infty\) norm as \(\|\cdot\|_\infty\). For any \(M>0\), let \(L^\infty(M)\) denote the collection of real-valued functions uniformly bounded by \(M\). For \(\Omega\subset \mathbb{R}^d\), we denote the uniform norm of a function \(f\) on \(\Omega\) by \(\|f\|_{L^\infty(\Omega)}\). For \(\mathbf{x}=(x_1,\cdots,x_d)\in\mathbb{R}^d\) and \({\boldsymbol{\alpha}}=(\alpha_1,\cdots,\alpha_d)\in\mathbb{N}^d\), we denote \(\mathbf{x}^{\boldsymbol{\alpha}}=x_1^{\alpha_1}x_2^{\alpha_2}\cdots x_d^{\alpha_d}\), and \(|{\boldsymbol{\alpha}}|=|\alpha_1|+|\alpha_2|+\cdots+|\alpha_d|\). Let \(\mathrm{D}_i\) be the derivative operator to the \(i\)-th variable for \(i\in[d]\), the multi-index derivative of order \({\boldsymbol{\alpha}}\in\mathbb{N}^d\) is defined as \(\mathrm{D}^{\boldsymbol{\alpha}}:=\mathrm{D}_{1}^{\alpha_1}\cdots \mathrm{D}_{d}^{\alpha_d}\). Let \(\beta=r+s\) for some non-negative integer \(r\) and \(0<s\leq 1\), \(d\in\mathbb{N}_+\), and \(C>0\). Let \(\|\cdot\|_{C^\beta(\Omega)}\) denote the Hölder norm of order \(\beta\), and a \(d\)-variate function \(f\) is called \((\beta,C)\)-smooth on \(\Omega\subseteq\mathbb{R}^d\), if \[\|f\|_{C^\beta(\Omega)}:=\max_{0\leq |{\boldsymbol{\alpha}}|\leq r}\|\mathrm{D}^{\boldsymbol{\alpha}}f\|_{L^\infty(\Omega)}\lor\max_{{\boldsymbol{\alpha}}:|{\boldsymbol{\alpha}}|=r}\sup_{\mathbf{x},\mathbf{y}\in\Omega, \mathbf{x}\neq \mathbf{y}}\frac{|\mathrm{D}^{\boldsymbol{\alpha}}f(\mathbf{x})-\mathrm{D}^{\boldsymbol{\alpha}}f(\mathbf{y})|}{\|\mathbf{x}-\mathbf{y}\|_2^s}\leq C.\]

2 Uniform Convergence of ReLU FNN Estimators↩︎

In this section, we investigate the fundamental limitations of ReLU FNN estimators in achieving uniform convergence. Section 2.1 establishes a theoretical lower bound on the uniform convergence rates of least-squares ReLU estimators in Theorem 1, demonstrating that they inherently suffer from the curse of dimensionality, even when their \(L^2(P)\) convergence rate is minimax-optimal. To explain the mechanism behind this bottleneck, Section 2.2 reveals how the limited smoothness of ReLU FNNs leads to slow uniform convergence rates.

2.1 A Lower Bound on Uniform Convergence Rates↩︎

To establish that the curse of dimensionality in uniform convergence is an inherent limitation of ReLU FNN estimators, we focus on the foundational framework of nonparametric least-squares regression. Let \(\{(X_i,Y_i)\}_{i=1}^n\) be the i.i.d.observations generated from \[\label{eq:32definition32of32least32squares32regression} Y_i = f_0(X_i) + \xi_i,\tag{1}\] where \(f_0\) is the unknown regression function, \(\xi_i\) denotes the random noise, and \(P\) represents the distribution of \((X_i,Y_i)\).

To estimate \(f_0\), we consider empirical risk minimization over a class of ReLU FNNs. Specifically, given a network width \(W\in\mathbb{N}_+\) and depth \(D\in\mathbb{N}_+\), a ReLU FNN is defined as a function \(f:\mathbb{R}^d\to\mathbb{R}\) of the composition form \[\label{eq:32definition32of32feedforward32network} f(x) = \mathcal{L}_{D+1}\circ\sigma_\ast\circ\mathcal{L}_{D}\circ\sigma_\ast\circ\cdots\circ\mathcal{L}_{2}\circ\sigma_\ast\circ\mathcal{L}_{1}(x),\tag{2}\] where \(\sigma_\ast(x)=x\lor 0\) is the ReLU activation function applied entrywise, and \(\mathcal{L}_i:\mathbb{R}^{d_{i-1}}\to \mathbb{R}^{d_i}\) are affine transformation maps with the dimension vector \((d_0,\ldots,d_{D+1})=(d,W,\ldots,W,1)\).

To rule out pathological estimators that interpolate the noise through increasingly steep local oscillations, we restrict attention to uniformly bounded ReLU networks with uniformly bounded parameters. Accordingly, for some \(r\ge 2\) and \(M\ge 0\), we define the candidate class by \[\begin{align} \mathcal{F}_{\mathrm{ReLU}}&:=\bigl\{f\text{ of the form~\eqref{eq:32definition32of32feedforward32network}: }\max_{i\in[D+1]}\|A_i\|_{\infty}\lor \|b_i\|_{\infty}\leq r\bigr\}\cap L^\infty(M). \end{align}\] Let the corresponding least-squares ReLU estimator \(\hat{f}_n\) be defined by \[\label{eq:32definition32of32ReLU32least32squares32estimator} \hat{f}_n \in \mathop{\arg\min}_{f\in\mathcal{F}_{\mathrm{ReLU}}} \sum_{i=1}^n \bigl(f(X_i)-Y_i\bigr)^2,\tag{3}\] We further require that the estimator sequence is uniformly Lipschitz, in the sense that \(\sup_{n\ge 1}\|\hat{f}_n\|_{\text{Lip}}<\infty\). These restrictions are imposed to exclude the possibility that poor uniform convergence arises from pathological overfitting to the noise.

Meanwhile, to isolate the architectural limitation of ReLU FNNs from pathological distributions, we also impose standard regularity conditions on the data-generating process. Specifically, it is assumed in Theorem 1 below that the covariates \(X_i\) are uniformly distributed on \([0,1]^d\), the noises \(\xi_i\) are independent of \(X_i\) and uniformly bounded, and the target function \(f_0\) is Hölder-smooth and uniformly bounded. These assumptions describe a highly regular nonparametric estimation setting: the covariate distribution is benign, the regression function is smooth, and the noise is bounded. Therefore, any slow uniform convergence rate established below cannot be attributed to irregular features of the data-generating process.

Taken together, the restrictions on the estimator and on the data-generating distribution ensure that the lower bound in Theorem 1 reflects an intrinsic limitation of the ReLU architecture itself. In particular, the slow uniform convergence rate is not caused by heavy-tailed noise, pathological covariate design, lack of smoothness of \(f_0\), or overfitting to the observations, but rather by the intrinsic limitations of ReLU FNN architectures. Under these conditions, we obtain the following lower bound on the uniform convergence rate of \(\hat{f}_n\).

Theorem 1 explicitly formalizes the classical curse of dimensionality for the least-squares ReLU estimator \(\hat{f}_n\), since the exponent of its uniform convergence rate depends poorly on dimension \(d\). In contrast, [4] proved that the \(L^2(P)\) convergence rates of ReLU FNN estimators do not have this issue. Specifically, under these conditions, the \(L^2(P)\) estimation error of ReLU estimator \(\hat{f}_n\) satisfies \[\sup_{P\in\mathcal{P}}\mathbb{E}\|\hat{f}_n-f_0\|_{L^2(P)} \;\lesssim_{\log}\; \sqrt{\frac{(WD)^2}{n}} \;+\; (WD)^{-2\alpha/d}.\] In particular, provided the target function \(f_0\) is not excessively smooth (i.e., \(\alpha \le d(d-2)/6\)), one can choose a network architecture with \(W\asymp D \asymp n^{\frac{d}{4(d+2\alpha)}} \geq n^{\frac{3}{4(d+1)}}\). Under this configuration, \(\hat{f}_n\) simultaneously attains the minimax-optimal \(L^2(P)\) rate \(n^{-\frac{\alpha}{2\alpha+d}}\) while remaining trapped at the slow, dimension-dependent rate of \(n^{-\frac{1}{d+1}}\) in the uniform norm. This stark discrepancy highlights a critical intrinsic limitation that least-squares ReLU FNN estimators can be suboptimal for downstream applications requiring uniform convergence guarantees. We provide a more detailed discussion of Theorem 1 and its implications in Section S1 of the Supplementary Material.

2.2 An Interpolation Perspective on Uniform Convergence Rates↩︎

In Section 2.1, ReLU least squares estimators are shown to suffer from the curse of dimensionality under the uniform norm. In contrast, many classical nonparametric estimators, such as splines, wavelets, and Fourier series, typically avoid this issue. A key reason is their linear structure: these estimators lie in the linear span of basis functions and, conditioned on the design matrix, behave analogously to ordinary least squares (OLS) estimators.

To illustrate this structural contrast, consider observations \(\{(X_i,Y_i)\}_{i=1}^n\) generated from the parametric linear model \(Y_i = X_i^\top \beta_0 + \xi_i\), where \(X_i\) is supported on \([0,1]^d\), and let \(\hat{\beta}_n\) denote the OLS estimator. Define the true and estimated regression functions by \(f_0(x) = x^\top \beta_0\) and \(\hat{f}_n(x) = x^\top \hat{\beta}_n\), respectively. Assuming the population second-moment matrix \(\Sigma:=\mathbb{E}[XX^\top]\) is positive definite with the smallest eigenvalue \(\lambda_{\min}(\Sigma)>0\), the \(L^2(P)\) and \(L^\infty\) estimation errors are deterministically linked. Specifically, because \(\|x\|_2 \le \sqrt{d}\), Cauchy–Schwarz implies \(\|\hat{f}_n - f_0\|_{L^\infty([0,1]^d)} \le \sqrt{d} \|\hat{\beta}_n - \beta_0\|_2\). Concurrently, \(\|\hat{f}_n - f_0\|_{L^2(P)}^2 = (\hat{\beta}_n - \beta_0)^\top \Sigma (\hat{\beta}_n - \beta_0) \ge \lambda_{\min}(\Sigma) \|\hat{\beta}_n - \beta_0\|_2^2\). Combining these relations yields the inequality that implies that the \(L^\infty\) and \(L^2(P)\) estimation errors are equivalent for linear models: \[\sqrt{\frac{\lambda_{\min}(\Sigma)}{d}} \|\hat{f}_n - f_0\|_{L^\infty([0,1]^d)} \leq \|\hat{f}_n - f_0\|_{L^2(P)} \leq \|\hat{f}_n - f_0\|_{L^\infty([0,1]^d)}.\]

Beyond parametric linear models, analogous interpolation inequalities serve as powerful classical tools to establish uniform convergence rates for broader function classes by exploiting underlying smoothness properties. For example, the Gagliardo–Nirenberg inequality [28], [29] applies to Hölder functions, while the Brezis–Mironescu inequality [30] applies to (fractional) Sobolev functions.

This viewpoint provides an analytic explanation for the slow uniform convergence rates observed in ReLU FNN estimators. To see this, consider the standard Gagliardo–Nirenberg inequality on \([0,1]^d\): for any function \(f \in C^\alpha([0,1]^d)\), the uniform norm is bounded by \[\label{eq:32Interpolation32inequality32for32Hölder-smooth32functions} \|f\|_{L^\infty([0,1]^d)} \lesssim \|f\|_{C^\alpha([0,1]^d)}^{\frac{d}{2\alpha+d}} \cdot \|f\|_{L^2(P)}^{\frac{2\alpha}{2\alpha+d}}.\tag{4}\] This interpolation bound shows that when the candidate function class possesses substantial smoothness (\(\alpha \gg d\)), the \(L^\infty\) convergence rate closely tracks the \(L^2(P)\) rate. However, because ReLU FNNs are generally only piecewise affine and need not belong to \(C^1\), this interpolation argument can use at most Hölder exponents \(\alpha\in(0,1)\). Applying 4 together with the standard \(L^2(P)\) rate \(\|\hat{f}_n - f_0\|_{L^2(P)} \lesssim n^{-1/2}\), and then letting \(\alpha\to1\), yields \[\|\hat{f}_n - f_0\|_{L^\infty([0,1]^d)} \lesssim n^{-1/(d+2)}.\] When the Gagliardo–Nirenberg inequality is tight, this bound is sharp. This issue shows that the poor intrinsic smoothness of the ReLU activation function fundamentally prevents its network architecture from adapting to fast uniform convergence requirements, rendering it ill-suited for downstream tasks requiring strong uniform convergence guarantees.

3 Smoothly Activated Deep Neural Networks↩︎

This section introduces smoothly activated deep neural networks (smooth DNNs), which encompass both deep feedforward and deep residual architectures, as a powerful and theoretically rigorous alternative. Section 3.1 reviews smooth activation functions widely adopted in modern deep learning and then summarizes the foundational framework of residual networks. Leveraging these components, Section 3.2 presents the mathematical formulation and characterization of smooth DNNs.

3.1 Preliminaries↩︎

We first characterize the infinitely differentiable (i.e., \(C^\infty\)-smooth) activation functions and the structural configurations of residual networks.

Smooth Activation Functions. Nonlinear activation functions govern the expressiveness of neural networks. While the ReLU activation serves as a baseline for deep compositional mapping, its non-differentiability at the origin constrains its functional regularity. Non-smooth variants such as Leaky ReLU, PReLU, ELU, and SELU introduce non-zero components on the negative real line but remain non-differentiable, making them analytically ill-suited for establishing fast uniform convergence rates under the interpolation perspective described in Section 2.2.

To achieve higher-order functional smoothness, we consider a class of \(C^\infty\)-smooth activations, including SiLU (also known as Swish) [13], [14], GELU [15], and Mish [31]. As illustrated in Figure 1, these activations exhibit smooth transitions near the origin. Formally, these \(C^\infty\)-smooth activation functions admits a factorization of the form \(\sigma(x)=x\cdot \psi(x)\), where \(\psi: \mathbb{R}\to \mathbb{R}\) is a non-decreasing, uniformly bounded, and infinitely differentiable function. The explicit functional forms are given by: \[\begin{align} \text{SiLU}(x) =& x\cdot \frac{\exp(x)}{1+\exp(x)}; \tag{5}\\ \text{GELU}(x) =& x\cdot \frac{1}{2}\Bigl(1+\tanh\bigl(\sqrt{2/\pi}\,(x+0.044715\,x^3)\bigr)\Bigr); \tag{6}\\ \text{Mish}(x) =& x\cdot \tanh\!\bigl(\log(1+\exp(x))\bigr). \tag{7} \end{align}\] Here we use the hyperbolic tangent version of the GELU activation, and its associated function \(\psi\) approximates the distribution function of the standard normal distribution. Writing \(\mathrm{ReLU}(x)=x\cdot H(x)\) with \(H(x)=\mathbb{I}(x\geq 0)\), the corresponding functions \(\psi\) of SiLU, GELU, and Mish converge to \(H(x)\) exponentially fast as \(x\to\pm\infty\). As shown in Section 4.2, this property ensures that the resulting network classes enjoy strong approximation properties for functions in Sobolev spaces and hierarchical composition models.

Residual Networks. In the context of nonparametric function estimation, ReLU FNNs have been shown to possess strong approximation capabilities relative to classical nonparametric models, such as wavelets [1]–[3]. However, as compositional depth increases, standard deep feedforward architectures often exhibit ill-conditioned optimization landscapes, which can lead to higher empirical risk even under careful parameter initialization and normalization strategies [16].

To regularize the optimization landscape, residual networks (ResNets) incorporate identity skip connections around each affine–activation block [16], [17]. Rather than applying an isolated nonlinear transformation at each layer, the residual architecture models a functional perturbation of the preceding state. Formally, given the output of the \(\ell\)-th layer \(h_{\ell}(x)\), the \((\ell+1)\)-th layer is defined as \(h_{\ell+1}(x)\;=\;h_{\ell}(x)\;+\;\mathcal{F}_\ell\bigl(h_{\ell}(x)\bigr)\), where \(\mathcal{F}_\ell\) is typically implemented via a small stack of affine layers, normalizations, and nonlinear activations.

As illustrated in Figure 2, the identity skip connection creates a direct pathway \(h_{\ell}\!\to h_{\ell+1}\) for more stable information and gradient propagation across layers. By mitigating the structural degradation of gradients, residual architectures make DNNs easier to train. We will incorporate residual architectures into our theoretical framework as an optimization-oriented generalization of feedforward architectures.

3.2 Mathematical Formulation of Smooth DNNs↩︎

We now define the function classes of smooth DNNs under both feedforward and residual architectures. To parameterize these architectures simultaneously, we first introduce an activation operator that unifies coordinatewise mappings and Gated Linear Unit (GLU) formulations.

Next, we define smooth DNNs for both residual and feedforward architectures.

For a fixed activation function \(\sigma \in C^\infty(\mathbb{R})\), the functional class of smooth DNNs with depth \(D\), width \(W\), and mapping from \(\mathbb{R}^{d_1}\) to \(\mathbb{R}^{d_{D+2}}\) is given by: \[\mathcal{F}(d_1,d_{D+2},D,W,\sigma) = \bigl\{x\mapsto f(x): f \text{ is of the form \eqref{eq:32definition32of32smooth32rn} or \eqref{eq:32definition32of32smooth32fnn}}\bigr\}.\] Unlike ReLU networks, the compositional structure of smooth DNNs ensures that any \(f \in \mathcal{F}(d_1, d_{D+2}, D, W, \sigma)\) inherits the \(C^\infty\) regularity of its activation primitives, rendering the empirical process amenable to high-order metric entropy bounds. Apart from the architectural change, the estimation procedure remains the same as for ReLU networks. Specifically, given a loss function \(\ell: \mathbb{R}^{d_y} \times \mathbb{R}^{d_y} \to \mathbb{R}_+\) and i.i.d.observations \(\{X_i,Y_i\}\in \mathbb{R}^{d_x}\times \mathbb{R}^{d_y}\) for \(i\in[n]\), we define the smooth DNN estimator by \[\hat{f}_n \;\in\; \mathop{\arg\min}_{f\in\mathcal{F}(d_x,d_y,D,W,\sigma)} \frac{1}{n}\sum_{i=1}^n \ell\big(f(X_i),\, Y_i\big).\]

4 Theoretical Properties of Smooth DNNs↩︎

This section develops the core analytic properties of the smooth DNNs. Specifically, the pseudo–dimension (Pdim) is quantified in Section 4.1. The approximation error bounds for Sobolev functions and hierarchical composition models are established in Section 4.2. Hölder norm upper bounds for approximators are provided in Section 4.3.

4.1 Pseudo-Dimension↩︎

In statistical learning theory, the statistical (generalization) error of \(\hat{f}_n\) is governed by the complexity of the candidate class \(\mathcal{F}(d_1,d_{D+2},D,W,\sigma)\), with higher complexity typically yielding slower rates. Complexity analysis for neural networks began with binary threshold activations and was later extended to real-valued networks via the Pdim [32]. More recently, nearly tight Pdim bounds have been established for networks with piecewise-polynomial activations, most notably ReLU [33]. While these results largely tie complexity control to specific activation classes, [34] bounded complexity in terms of the operation number needed to compute the network output. Building on this idea, [35] showed that when the neural network is describable by Pfaffian functions, the Pdim admits polynomial upper bounds in the number of parameters. This framework is particularly convenient for establishing Pdim bounds of smooth DNNs. Hence, we first recall the definition of a Pfaffian function.

Indeed, the smooth activations of interest in modern deep learning practice are Pfaffian.

We next recall the definition of Pdim for real-valued function classes.

We now present an upper bound on the Pdim of smooth DNNs with Pfaffian activations.

4.2 Approximation Error↩︎

As the smooth DNN becomes deeper or wider, its increased Pdim may incur larger statistical error, while the increased capacity can reduce approximation error by representing more complex targets. Thus, non-asymptotic approximation bounds for smooth DNNs are required to identify architectural scalings that balance these two errors.

Building on approximation theory for ReLU FNNs [36], [37], recent studies have extended approximation guarantees to smooth FNNs. In particular, [38] showed that any ReLU FNN can be approximated arbitrarily well by FNNs for a broad class of smooth activations. More recently, [39] leveraged linear combinations of smooth activations to approximate monomials and derived non-asymptotic approximation error bounds for such networks. Following the strategy in [39], we establish non-asymptotic approximation error bounds for smooth DNNs when approximating Sobolev functions and hierarchical composition models. We begin by stating the analytic assumptions on the activation \(\sigma\).

This assumption ensures that smooth DNNs can approximate monomials and two specific functions that are central in neural network approximation theory. For the \(C^\infty\) activations of interest in this work, including SiLU, GELU, and Mish, the local \(C^\infty\) assumption is satisfied. Since these activations admit a uniformly bounded and increasing function \(\psi(x)=\sigma(x)/x\) with exponentially decaying derivatives, the \(\mathfrak{m}\)-th order quasi-decay assumption also holds.

We now present the approximation error bounds of smooth DNNs for Sobolev functions.

This result establishes approximation rates for smooth DNNs when approximating Sobolev functions. In what follows, we establish non-asymptotic approximation error bounds for smooth DNNs under hierarchical composition models.

Leveraging results in Theorem 4, we obtain the following approximation error results for smooth DNNs under hierarchical composition models.

Theorem 5 shows that when the unknown target \(f_0\) admits a hierarchical composition of several \((\beta,C)\)-smooth functions, each with either high smoothness or low input dimension, smooth DNNs can adapt to this intrinsic low-dimensional structure without explicit knowledge of the underlying composition. In particular, smooth DNNs achieve the same approximation bounds as ReLU FNNs under the hierarchical composition model [1]–[3].

4.3 Upper Bounds of Hölder Norms↩︎

As discussed in Section 2.2, uniform convergence rates are governed by the smoothness of the candidate function class. Consequently, to establish uniform convergence rates for smooth DNN estimators, we need to characterize their Hölder norms across different orders.

We first examine the analytic properties of the activations SiLU, GELU, and Mish.

Proposition 6 shows that the derivatives of the smooth activations of interest grow at a factorial rate. As these activations are the building blocks of smooth DNNs, this result provides a natural baseline for the order of the Hölder-norm bounds that one can expect for smooth DNNs. The following result further establishes factorial-order Hölder norm bounds for smooth DNN approximators on a subset of \(\Omega\) that is arbitrarily close to \(\Omega\).

This result shows that the smooth DNN approximator admits a factorial-order upper bound on the Hölder norm \(C^{\mathfrak{m}}\), as that for analytic activations in Proposition 6, even when \(\mathfrak{m}>\alpha\) exceeds the smoothness order of the target function. Moreover, this result does not require \(\sigma\) to be analytic, indicating that the factorial-order upper bound originates from the intrinsic structure of smooth DNNs rather than from the analyticity of \(\sigma\).

Rather than establishing a factorial-order upper bound on the entire domain \(\Omega=[0,1]^d\), we obtain such a bound only on a subset \(\Omega_{\delta}\), whose measure can be made arbitrarily close to that of \(\Omega\). This restriction arises because the smooth DNN approximator \(\phi\) involves certain sub-networks whose Sobolev norm control is difficult to establish on small regions of the domain. That said, we believe that the present result is sufficient to explain the well-behaved Hölder norms of smooth DNN approximators.

5 Uniform Convergence of Smooth DNN Estimators↩︎

This section derives uniform convergence rates for smooth DNN estimators under several widely used models. Specifically, in Section 5.1, we study smooth DNN estimators under the Huber regression model, establish their robustness to heavy-tailed noise, and consider least-squares regression as a special case of Huber regression. We then investigate quantile regression with smooth DNN estimators in Section 5.2. Uniform convergence rates for smooth DNN estimators under logistic regression are established in Section 5.3.

Before turning to these models, we assume the covariate \(X_i\) is supported on \(\Omega:=[0,1]^d\) and admits a bounded density throughout this section. In particular, there exist constants \(0<c<C<\infty\), such as the density \(p_X\) of \(X_i\) satisfies \[0<c<p_X(x)\leq C<\infty,\quad\text{for all }x\in[0,1]^d.\] We also impose the following assumption on the activation function \(\sigma\) to invoke the theoretical properties of smooth DNNs established in Section 4. SiLU, GELU, and Mish are key examples satisfying these conditions.

To leverage interpolation inequalities when deriving uniform convergence rates, we impose the following assumption on the Hölder norms of the candidate function class. As discussed in Section 4.3, the rationale is that both analytic activations and smooth DNN approximators to Sobolev functions satisfy factorial-type upper bounds on their Hölder norms. Thus, this condition can be interpreted as restricting attention to a bounded subset of smooth DNNs in a neighborhood of the approximator.

5.1 Huber Regression and Least Squares↩︎

In this section, we study the uniform convergence of smooth DNN estimators in Huber regression and least squares. Specifically, let \(\{(X_i,Y_i)\}_{i=1}^n\sim P\) be i.i.d.observations satisfying \[\label{eq:32definition32of32regression32data32generating32process} Y_i = f_0(X_i) + \xi_i, \qquad \text{where} \quad \mathbb{E}[\xi_i \mid X_i] = 0,\tag{10}\] where \(f_0\) denotes the target function of interest.

Motivated by the prevalence of heavy-tailed data in real-world applications, a growing line of work has investigated the \(L^2\) estimation error for nonparametric function estimators under heavy-tailed noise [4], [5], [40], [41]. In contrast to light-tailed noise, heavy-tailed noise produces outliers in \(Y_i\) more frequently, making classical least-squares estimators unreliable. To address this issue, the Huber loss was introduced as a robust alternative, defined as, for the Huber parameter \(\tau\geq 0\), \(\ell_\tau(x) = \frac{1}{2}x^2\mathbb{I}(|x|\leq\tau) + (\tau|x| - \frac{1}{2}\tau^2)\mathbb{I}(|x|>\tau).\) Then, for a function class \(\mathcal{F}_n\), the Huber estimator is defined as \[\label{eq:32definition32of32deep32huber32estimator} \hat{f}_n(\tau)\in\mathop{\arg\min}_{f\in\mathcal{F}_n}\frac{1}{n}\sum_{i=1}^n \ell_\tau\big(Y_i-f(X_i)\big).\tag{11}\]

Recent advances have characterized the robustness of the Huber estimator \(\hat{f}_n(\tau)\) for different function classes with non-asymptotic analysis. When \(\mathcal{F}_n\) consists of linear functions, [42] showed that the adaptive Huber estimator achieves a sub-Gaussian concentration bound in the presence of heavy-tailed noise. More recently, analogous robustness guarantees for ReLU FNN estimators were established in prior works [4], [5]. However, these results are formulated in the \(L^2(P)\) norm rather than in \(L^\infty\), leaving open whether DNN estimators remain robust in the uniform norm under heavy-tailed regimes.

The assumption below is imposed to show that smooth DNN Huber estimators are robust in the uniform norm under heavy-tailed noise.

We now establish robustness guarantees for the estimator \(\hat{f}_n(\tau)\) in the \(L^\infty(\Omega)\) norm.

In this result, the uniform estimation error of the smooth DNN Huber estimator \(\hat{f}_n(\tau)\) is decomposed into four components: (i) statistical error \(\delta_\mathfrak{s}\), (ii) approximation error \((WD)^{-2\gamma^\ast}\), (iii) Huberization bias \(v_m/\tau^{m-1}\), and (iv) a sub-Gaussian deviation term. These components reveal two natural trade-offs. First, increasing \(\tau\) reduces the Huber bias but enlarges both the deviation term and \(\delta_\mathfrak{s}\), since the Huber loss approaches the least-squares loss. Second, enlarging the network architecture (in terms of \(D\) and \(W\)) decreases the approximation error but reduces the effective sample size \(\tilde{n}\), thereby increasing the statistical error. Moreover, this \(L^\infty\) estimation error decomposition parallels that of ReLU Huber estimators in the \(L^2\) norm [4], [5], with the only distinction arising from the pseudo-dimensions of the two function classes. This comparison highlights a concrete advantage of smooth DNN estimators over existing ReLU FNNs in terms of uniform convergence guarantees.

We now turn to least-squares regression as a special case of Huber regression. As shown in Section 2, ReLU FNN least-squares estimators suffer from the curse of dimensionality even under idealized settings, in sharp contrast to their favorable \(L^2\) convergence behavior established in [4], [5]. Since the least-squares loss can be viewed as the limiting case of the Huber loss \(\ell_\tau\) as \(\tau\to\infty\), Theorem 8 suggests the corresponding uniform convergence behavior of smooth DNN least-squares estimators under heavy-tailed noise. However, one cannot directly let \(\tau\to\infty\) in Theorem 8, because the final term, namely the sub-Gaussian deviation bound, diverges in this limit. This is precisely because least-squares regression is not robust to heavy-tailed noise, and therefore cannot satisfy a non-asymptotic heavy-tailed robustness guarantee of sub-Gaussian type. For this reason, in Theorem 9 below, we formulate the uniform estimation error of smooth DNN least-squares estimators in terms of convergence in expectation rather than a non-asymptotic deviation bound.

While Theorem 1 shows that the convergence rate of ReLU FNN least-squares estimators can be as slow as \(n^{-\frac{1}{d+1}}\) regardless of the smoothness of \(f_0\), smooth DNN least-squares estimators can adapt to the hierarchical composition structure of \(f_0\), so that their convergence rates depend on the intrinsic dimensionality of its low-dimensional representation. When the noise is light-tailed, as in Theorem 1 (in the sense that \(m\to\infty\)), the resulting rate in Theorem 9 scales as \(n^{-\frac{0.99\,\gamma^\ast}{2\gamma^\ast+2}}\). Thus, smooth DNN estimators can mitigate the curse of dimensionality for least-squares regression in the uniform norm by adapting to the low-dimensional hierarchical composition structure of \(f_0\).

5.2 Quantile Regression↩︎

As another robust estimation task, quantile regression estimates conditional quantile functions [43]. For a fixed quantile level \(\tau \in (0,1)\), suppose we observe i.i.d.samples \(\{(X_i, Y_i)\}_{i=1}^n\) such that the \(\tau\)-th conditional quantile of \(Y_i\) given \(X_i = x\) is \(f_0(x)\). Then \(f_0\) minimizes the population quantile regression loss: for the check loss \(\rho_\tau(u) = u(\tau - \mathbb{I}(u < 0))\), \[f_0 \in \mathop{\arg\min}\mathbb{E}[\rho_\tau(Y_1 - f(X_1))].\] Recent work has established non-asymptotic \(L^2(P)\) convergence rates and non-asymptotic robustness guarantees for ReLU FNN quantile regression estimators [5], [44]. In the following, we establish robustness guarantees for smooth DNN quantile estimators in the \(L^\infty\) norm. Specifically, we introduce the following regularity condition to ensure well-behaved conditional densities in a neighborhood of the target quantile:

Leveraging the theoretical results developed in Section 4, we establish robustness guarantees for the smooth DNN quantile regression estimator in the \(L^\infty\) norm below.

Similar to Theorem 8, the exponent \(0.99\) is flexible and can be replaced by any value strictly smaller than \(1\). Moreover, by optimizing over the depth and width of smooth DNNs in terms of \((D,W)\), the convergence rate of the smooth DNN quantile estimator scales as \(\mathbb{E}\|\hat{f}_n-f_0\|_{L^\infty(\Omega)}\lesssim_{\log} n^{-\frac{0.99\,\gamma^\ast}{2+2\gamma^\ast}}.\) Notably, this convergence exponent does not depend on the ambient dimension \(d\), but only on the adjusted model-complexity parameter \(\gamma^\ast\). This justifies that the smooth DNN estimator can mitigate the curse of dimensionality by adapting to the low-dimensional hierarchical composition structure of \(f_0\).

5.3 Logistic Regression↩︎

In this subsection, we first establish uniform estimation error bounds for smooth DNN logistic regression estimators in Theorem 11. Building on this result, Theorem 12 further establishes uniform convergence guarantees for the density-ratio-based distribution estimator studied in [27]. For i.i.d.observations \(\{(X_i,Y_i)\}_{i=1}^n\), define the smooth DNN logistic regression estimator by \[\hat{f}_n\in \mathop{\arg\min}_{f\in\mathcal{F}_n}\frac{1}{n}\sum_{i=1}^n -Y_i\cdot f(X_i)+\log(1 + \exp(f(X_i)))=:\frac{1}{n}\sum_{i=1}^n\ell(f;X_i,y_i),\] where \(\mathcal{F}_n = \mathcal{F}(d,1,D,W,\sigma;A)\cap L^\infty(M)\). Let the target function \(f_0\) satisfy \[\mathbb{E}[Y_i|X_i]= \frac{\exp(f_0(X_i))}{1 + \exp(f_0(X_i))}.\]

As in the cases of Huber and quantile regression, convergence rates for the smooth DNN logistic regression estimator can be established below by leveraging the theoretical results developed in Section 4.

Theorem 11 shows that the smooth DNN logistic regression estimator enjoys an estimation error bound analogous to that of the smooth DNN quantile regression estimator established in Theorem 10. This further supports the reliability of smooth DNN estimators across a range of regression and classification models.

The uniform convergence guarantee established in Theorem 11 is also useful for supporting the practical application of smooth DNNs to estimating the probability of Higgs boson pair production decaying into four bottom quarks (denoted by \(HH\to 4b\)) [27]. We briefly reformulate this problem within the nonparametric framework considered in [27] and then state the corresponding uniform convergence guarantee.

In this application, the event \(HH\to 4b\) is the signal of interest, while other physical processes producing four bottom quarks form the background. Accurate estimation of both signal and background distributions is needed to distinguish signal from background and conduct subsequent hypothesis testing. While the signal distribution can be approximated using physical simulation, simulating the background distribution is computationally prohibitive. Consequently, the main problem in [27] is to estimate the background distribution from observed data contaminated by signal events.

To this end, [27] introduced an auxiliary distribution together with density ratio estimation. Let \(\Omega\subseteq\mathbb{R}^d\) denote the state space of the observed physical quantities, and let \(P\) and \(Q\) denote the probability measures of the background and auxiliary distributions on \(\Omega\), respectively. In that work, \(Q\) corresponds to events with three observed bottom quarks, which share similar physical characteristics with background events. Density ratio estimation is reduced to a binary classification problem: let \(Y\) be binary and let \(E\) satisfy \[E\mid Y=0\sim Q,\qquad E\mid Y=1\sim P.\] Define \[\psi(x):=\mathbb{P}(Y=1\mid E=x),\qquad x\in\Omega.\] Then Bayes’ rule gives, when \(\mathbb{P}(Y=1)=\mathbb{P}(Y=0)=0.5\), \[P(A)=\int_A \frac{\psi(y)}{1-\psi(y)}\,Q(\text{d}y),\qquad \text{for any measurable }A\subseteq\Omega.\] Hence, if \(\hat{\psi}\) is an estimator of \(\psi\) and \(\mathbb{Q}_N\) is an empirical measure of \(Q\), then a plug-in estimator of \(P\) is \[\hat{P}(A)=\int_A \frac{\hat{\psi}(y)}{1-\hat{\psi}(y)}\,\mathbb{Q}_N(\text{d}y),\qquad \text{for any measurable }A\subseteq\Omega.\] Similarly, if \(p\) and \(q\) are the densities of \(P\) and \(Q\), and \(\hat{q}\) is an estimator of \(q\), then \[\hat{p}(x)=\frac{\hat{\psi}(x)}{1-\hat{\psi}(x)}\,\hat{q}(x),\qquad x\in\Omega.\] In [27], the density ratio \(\frac{\psi(y)}{1-\psi(y)}\) is estimated using a SiLU ResNet tailored to the structure and symmetries of the collider data.

A key difficulty is that the observed four-bottom-quark data are contaminated by signal events, so one cannot directly fit a logistic regression model on the full domain. Following [27], suppose \[\Omega=\mathcal{X}\cup\mathcal{Y},\qquad \mathcal{X}\cap\mathcal{Y}=\emptyset,\] where signal events occur only on \(\mathcal{X}\). Then contamination is present only on \(\mathcal{X}\), while the observed four-bottom quark events on \(\mathcal{Y}\) are purely background. Therefore, one may estimate the classifier between background and three-bottom quark events in \(\mathcal{Y}\), and then extrapolate it to the full domain \(\Omega\). Although this extrapolation assumption is statistically strong, it is reasonable in this physical setting and can be cross-checked against an alternative optimal-transport-based approach. To establish theoretical guarantees for [27], we adopt this framework in a slightly relaxed form.

Compared with Assumption 1 in [27], which requires the binary classifier on \(\mathcal{Y}\) to coincide exactly with the binary classifier on the whole domain \(\Omega\), the above assumption allows the two population risk minimizers to differ and quantifies the resulting discrepancy between the corresponding population loss functions through \(\varepsilon_n\). In particular, when \(\widetilde{P}\) and \(\widetilde{P}_\mathcal{Y}\) exhibit similar physical patterns for distinguishing four-bottom-quark events from three-bottom-quark events, one expects \(\varepsilon_n\) to be small.

Leveraging Assumption 6 together with the uniform convergence guarantee in Theorem 11, we establish the following result for the density-ratio-based distribution estimators.

Theorem 12 transfers the uniform convergence rate of the smooth DNN logistic regression estimator to uniform convergence guarantees for the density-ratio-based estimators \(\hat{P}\) and \(\hat{p}\) of the background distribution. This provides a theoretical justification for the density-ratio-based method in [27] under the relaxed extrapolation condition in Assumption 6.

In particular, the theorem yields non-asymptotic estimation error bounds for both \(\hat{P}\) and \(\hat{p}\). Each bound consists of four components: the statistical error term \(\sqrt{\log(10/\delta)/\tilde{n}}\), the approximation error term \((WD)^{-2\gamma^\ast}\), the extrapolation error \(\varepsilon_n\), and the estimation error of the auxiliary distribution, represented by \(|(Q-\mathbb{Q}_N)\exp(f_0(\cdot))\mathbb{I}(A)|\) and \(\|\hat{q}-q\|_{L^\infty(\Omega)}\), respectively. The first two terms arise from smooth DNN logistic regression, while the term \(\varepsilon_n\) quantifies the discrepancy introduced by extrapolation, which is expected to be small in this physical setting. Moreover, when the auxiliary distribution can be simulated, the last term is typically negligible in this setting. Besides, if a collection of measurable sets \(\mathcal{A}\) is sufficiently regular so that \[\sup_{A\in\mathcal{A}}\Bigl|(Q-\mathbb{Q}_N)\exp(f_0(\cdot))\mathbb{I}(A)\Bigr|\to 0\qquad\text{as }N\to\infty,\] and if \(\|\hat{q}-q\|_{L^\infty(\Omega)}\to 0\) (e.g., permitted by kernel density estimation [45]), the auxiliary-distribution estimation error can be made asymptotically negligible for sufficiently large \(N\).

6 Numerical Experiments↩︎

Section 6.1 benchmarks SiLU FNN and ResNet Huber estimators against ReLU FNN estimators through simulation studies. In Section 6.2, we apply SiLU ResNet quantile estimators to analyze the temperature sensitivity of ozone concentrations across regions. All experiments were conducted on a computer running Ubuntu 24.04.2 LTS with an AMD EPYC 7532 CPU, 512 GB RAM, and an NVIDIA RTX 3090 GPU.¹

6.1 Simulation Studies↩︎

In this subsection, we compare the numerical performance of SiLU DNN and ReLU FNN estimators for Huber regression. We consider i.i.d.samples \(\{(X_i,Y_i)\}_{i=1}^n\) generated from \(Y_i = f_0(X_i) + g(X_i)\,\varepsilon_i,\) where the covariate \(X_i\in\mathbb{R}^8\) is sampled uniformly from \([0,1]^8\), and \(\varepsilon_i\) is drawn from a Student’s \(t\) distribution with either 2 or 4 degrees of freedom (\(t_2\) or \(t_4\)). Following [44] to induce noise, the target and scale functions are defined as \[\begin{align} f_0(x)=&\cos(2\pi x_1)+\frac{1}{1+\exp(-x_2-x_3)}+\frac{1}{(1+x_4+x_5)^3}+\frac{1}{x_6+\exp(x_7x_8)},\\ g(x)=&\sin\!\Big(\frac{\pi(x_1+x_2)}{2}\Big)+\log\!\big(1+x_3^2x_4^2x_5^2\big)+\frac{x_8}{1+\exp(-x_6-x_7)}. \end{align}\]

For the Huber estimator in 11 , we set the Huber parameter to \(\tau\in\{0.25, 0.5, 1, 2, 4\}\). The depth and width of the ReLU FNNs, SiLU FNNs, and SiLU ResNets are all set to \(4\). For the SiLU ResNets, we take \(\mathcal{L}_1\) and \(\mathcal{L}_3\) in Definition 2 to be identity maps so that all architectures have matched parameter sizes.

We consider sample sizes \(n\in\{512, 1024, 2048\}\) and replicate each setting 100 times. For evaluation, given an estimator \(\hat{f}_n\), we independently draw \(\{X_j\}_{j=1}^{10000}\) uniformly from \([0,1]^8\) and compute the pointwise errors \(\{|\hat{f}_n(X_j)-f_0(X_j)|\}_{j=1}^{10000}\). Based on these errors, we report the estimated \(L^2\) and \(L^\infty\) errors in Tables 1 and 2 under \(t_2\) and \(t_4\) noise, respectively. The corresponding side-by-side boxplots for the \(L^2\) and \(L^\infty\) metrics are presented in Figures 3 and 4, respectively. As shown in Tables 1–2 and Figures 3–4, the SiLU FNN and ResNet estimators consistently outperform the ReLU FNN estimator across all settings. Moreover, SiLU ResNet estimators exhibit smaller estimation error compared to SiLU FNN estimators, which is consistent with the optimization advantages of residual architectures. To further assess the statistical estimation performance, we additionally report the estimated convergence exponents. Specifically, under the scaling law \(n^{-b}\) for the estimation error, we estimate \(b\) via log–log regression and report the results in Tables 1 and 2, together with bootstrap standard deviations based on 1000 resamples. Overall, while all estimators yield comparable exponents for the \(L^2\) norm, ReLU FNNs have much smaller uniform convergence exponents. This empirical pattern is consistent with our theory: ReLU FNNs are not well suited to tasks requiring uniform convergence.

6.2 Real Data Analysis: Temperature Sensitivity of Ozone↩︎

Ozone is widely recognized as an air pollutant that affects human health [46]. The IPCC Sixth Assessment Report assesses that climate-driven changes in meteorology can influence ozone concentration, and that shifts in the geographical distribution of emissions have contributed to strong regional variations in ozone trends [47]. These assessments motivate us to study the relationship between temperature and ozone concentrations over the contiguous United States (CONUS) using quantile regression. From a methodological perspective, uniform convergence is important here because it supports reliable estimation uniformly over the spatial domain, rather than only in an average sense.

Our response is the daily maximum 8-hour average ozone concentration, measured in parts per billion (ppb), obtained from the U.S.EPA Air Quality System (AQS). The covariates include latitude, longitude, year, day-of-year, and daily maximum and minimum near-surface temperatures (tmmx, tmmn) in \(^\circ\)C at AQS monitoring sites. The temperature covariates are obtained from the gridMET dataset [48] and mapped to each monitoring site via nearest-neighbor assignment to the closest gridMET grid cell. We construct the dataset using daily observations from the warm season (May–September) over the years 2015–2024, when photochemical ozone formation is most active. The resulting dataset contains approximately 1.84 million observations from about 1,500 monitoring sites.

For this quantile regression task, we consider a SiLU ResNet and a ReLU FNN, both with depth 6 and width 128, and train them in PyTorch using the Adam optimizer for 200 epochs. As in Section 6.1, we take \(\mathcal{L}_1\) and \(\mathcal{L}_3\) for the SiLU RNs in Definition 2 to be identity maps so that the two architectures have comparable parameter sizes. We fit conditional quantiles at levels \(\tau\in\{0.10,0.30,0.50,0.87\}\), following the practice in ozone studies [49], since these quantiles admit meaningful practical interpretations. Specifically, the \(\tau=0.10\) quantile reflects relatively clean-condition days, the \(\tau=0.30\) quantile serves as a robust proxy for the photochemical baseline, the \(\tau=0.50\) quantile summarizes the typical ozone level, and the \(\tau=0.87\) quantile provides an upper-tail summary of broadly elevated ozone days during the warm season.

We then perform a counterfactual perturbation analysis using the trained models. Specifically, using temperatures from the warm season of 2025, we construct perturbed covariates by adding a constant warming increment \(\Delta\in\{0.5,1.0,1.5\}\) to both tmmx and tmmn, while holding the remaining covariates fixed. Using the original and perturbed covariates, we compute predicted daily ozone concentrations and report their differences over the CONUS at 20km resolution. Averaging the predictions over the warm season, Figure 5 presents the counterfactual differences in ozone concentrations across quantiles predicted by the SiLU ResNets over the CONUS.

Figure 5 reveals a spatially heterogeneous temperature sensitivity of ozone across the CONUS. Across all four quantile levels, warming induces broadly positive changes over the northern tier of the United States and California. These patterns are consistent with prior studies based on classical models [49], [50]. In contrast, ozone decreases over parts of the South and the Gulf Coast, consistent with prior evidence that, over tropical and subtropical oceans, increases in temperature and water vapor can reduce photochemical ozone production [51]. Taken together, these counterfactual patterns suggest that the SiLU ResNet estimators capture climate-driven and geographically heterogeneous temperature–ozone relationships.

To further compare the counterfactual prediction performance of the ReLU FNN and SiLU ResNet estimators, we complement the visual inspection with two quantitative diagnostics. The is Spatial Roughness (SR), measures the local irregularity of a rasterized counterfactual ozone-change surface through first-order finite differences; larger values indicate less spatial smoothness and more abrupt local variation. The second is Extreme Value Fraction (EVF), defined as the fraction of grid cells with extremely large values. This diagnostic measures the tendency of a fitted model to generate excessive local extremes. For each model, Table 3 reports the average diagnostic value over all combinations of quantile level \(\tau\) and warming increment \(\Delta\) within each month, as well as over the full warm season for the overall summary.

More precisely, let \(D_{\tau,\Delta,t}(s)\) denote the daily counterfactual ozone change at grid cell \(s\), for quantile level \(\tau\), warming increment \(\Delta\), and date \(t\). For each month \(m\), and for the overall warm season, we first average \(D_{\tau,\Delta,t}(s)\) over time to obtain a mean counterfactual field \(Z(s)\). Let \(Z\) be a rasterized counterfactual ozone-change field on a two-dimensional grid, represented by the array \((Z_{i,j})_{i,j}\). Let \(\mathcal{G}\subseteq\mathbb{N}_+\times\mathbb{N}_+\) denote the set of valid grid-cell indices, and let \(\mathcal{H},\mathcal{V}\subseteq\mathbb{N}_+\times\mathbb{N}_+\) denote the sets of valid grid cells whose horizontal and vertical neighbors are also valid, respectively, namely, \[(i,j)\in \mathcal{H}\text{ if } (i,j)\in\mathcal{G}\text{ and } (i,j+1)\in\mathcal{G};\qquad (i,j)\in \mathcal{V}\text{ if } (i,j)\in\mathcal{G}\text{ and } (i+1,j)\in\mathcal{G}.\] We define the spatial roughness of \(Z\) by \[\mathrm{SR}(Z)= \frac{\sum_{(i,j)\in \mathcal{H}}|Z_{i,j+1}-Z_{i,j}|+\sum_{(i,j)\in \mathcal{V}}|Z_{i+1,j}-Z_{i,j}|}{|\mathcal{H}|+|\mathcal{V}|},\] that is, the mean absolute difference over all valid horizontal and vertical neighboring pairs, and \[\mathrm{EVF}(Z;u)= \frac{1}{|\mathcal{G}|}\sum_{s\in\mathcal{G}} \mathbf{1}\{|Z(s)|>u\},\] with threshold \(u=3.0\) ppb. For each model and each month, we average these diagnostics over all quantile levels \(\tau\in\{0.10,0.30,0.50,0.87\}\) and warming increments \(\Delta\in\{0.5,1.0,1.5\}\). Standard deviations are computed using a month-stratified paired bootstrap over dates with 1000 resamples.

As shown in Table 3, the SiLU ResNet yields uniformly smaller SR values than the ReLU FNN, indicating more spatially coherent counterfactual predictions. For EVF, the SiLU ResNet also attains smaller values in most months and in the overall summary. In addition, the SiLU ResNet exhibits smaller bootstrap standard deviations for both diagnostics in all reported cases. Taken together, these results suggest that the SiLU ResNet produces more stable and spatially coherent quantile predictions, while generating fewer potentially spurious extreme local responses in this application.

7 Conclusion↩︎

In this work, we investigate the slow uniform convergence of ReLU FNN estimators and develop a statistical theory for smooth DNNs, an alternative class that aligns with modern deep learning practice while enjoying uniform convergence guarantees in many statistical learning tasks. We first show that, despite minimax-optimal \(L^2(P)\) convergence, least-squares ReLU FNN estimators can suffer from the curse of dimensionality in uniform convergence, owing to the low smoothness of the ReLU activation. Motivated by the prevalence of residual architectures and \(C^\infty\) activations in state-of-the-art deep learning models, we introduce smooth DNNs and establish their core statistical properties, including pseudo-dimension bounds, approximation error bounds for Sobolev functions and hierarchical composition models, and factorial-type Hölder norm bounds for approximators. Building on these foundations, we derive non-asymptotic \(L^\infty\) convergence rates for smooth DNN estimators in several representative tasks: Huber, least-squares, quantile, and logistic regression, showing that smooth DNNs can provide reliable uniform convergence guarantees while adapting to low-dimensional hierarchical structure. Together, these results position smooth DNNs as a theoretically grounded alternative to ReLU FNN estimators in settings where uniform convergence is essential.

References↩︎

---

1. The Python code used to reproduce both the simulation and real-data experiments is available in an anonymized GitHub repository at https://anonymous.4open.science/r/Uniform_convergence-DF22/.↩︎