Nonparametric regression using deep neural networks with ReLU activation function

1. Introduction.

• Deep ReLU networks achieve near-minimax rates for arbitrarily smooth regression functions, despite producing piecewise-linear functions, under a hierarchical composition assumption that includes generalized additive models.This framework addresses the curse of dimensionality through structural assumptions on the regression function.Source: §1 · p. 2
• The theory indicates that ReLU network depth should increase with sample size, so larger samples naturally support additional hidden layers.The analysis is motivated by the growing depth of modern network architectures.Source: §1 · p. 2
• For networks with more potential parameters than observations, the number of nonzero parameters acts as the effective model dimension and must be carefully selected to control overfitting.The framework models large networks through sparsity rather than requiring the total parameter count to be below the sample size.Source: §1 · p. 3, §2 · pp. 6–7
• Network parameters can be restricted to absolute value at most one, matching the typically small learned weights and remaining compatible with projection-based training.The same bounded-parameter condition is incorporated into the formal network class.Source: §1 · p. 3, §2 · p. 6
• The estimator is minimax-rate optimal up to logarithmic factors if and only if its reconstruction method nearly minimizes empirical risk.The optimization gap measures the discrepancy between the estimator’s expected empirical risk and the global minimum over the network class.Source: §1 · p. 4, §2 · pp. 7–8
• Wavelet series estimators are suboptimal under the composition assumption, with lower bounds showing a polynomial loss in the sample size n relative to neural-network fitting.This identifies a function class where neural networks outperform wavelet estimators.Source: §1 · p. 4

2. Mathematical definition of multilayer neural networks.

• A network of depth L and width vector p alternates shifted ReLU transformations with matrix multiplications, with input dimension p0 = d and scalar output pL+1 = 1.The network parameters are the entries of the weight matrices and shift vectors.Source: §2 · pp. 5–6
• The formal class permits large depth and overparameterization, bounds every parameter by one, and models sparsity through only a small number s of nonzero parameters.Prediction performance is evaluated through expected squared-error risk, while empirical-risk minimization provides the optimization benchmark.Source: §2 · pp. 6–7

Main theorem

• Theorem 1 gives risk bounds for sparse network estimators over compositional regression classes, with rates governed by the effective smoothness and approximation error.The regression function is represented as a composition whose components have Hölder smoothness β_i and depend on at most t_i variables, allowing dimension-free exponents and including additive-model-type constraints.Source: §3 · pp. 8–10

Network architecture

• Choosing network depth L of order log^2 n minimizes the stated upper-rate expression, while sparsity s must be selected at the prescribed order.Network width can be chosen independently of smoothness under the stated condition, but the theorem requires sparse rather than fully connected architectures.Source: §3 · pp. 10–11

Oracle inequality

• Theorem 2 supplies an oracle inequality linking prediction risk to empirical risk and complexity, while increasing depth worsens the resulting upper bounds.The depth effect is consistent with observations that overly deep standard feedforward networks can degrade performance; residual networks are noted as a separate remedy requiring separate analysis.Source: §3 · pp. 11–12

ReLU properties

• ReLU networks support depth synchronization through identity propagation and bounded parameters, simplifying approximation constructions and avoiding the large weights required by standard smooth-activation square approximations.Identity propagation also underlies skip connections and residual-network behavior, whereas approximating the square function with other activations can require parameters diverging as sample size increases.Source: §3 · pp. 12–13

Minimax lower bounds

• Theorem 3 shows that φ_n is a minimax lower bound in the regime t_i ≤ min(d_0, ..., d_{i−1}), matching the upper rate up to logarithmic factors.The lower-bound analysis also shows that L2 and sup-norm minimax rates can differ by a polynomial power of n in some setups.Source: §3 · pp. 13–14

Approximation lower bounds

• Minimax lower bounds combined with the oracle inequality yield approximation-theoretic lower bounds on the number of weights required by sparse networks.This provides a statistical route to proving necessary network-size conditions for prescribed approximation accuracy, including refined bounds for L2 loss.Source: §3 · p. 14

Additive models

• n^−2β/(2β+1) log^3 n + Δ(bfn, f0) bounds additive-model risk, matching the minimax rate up to the log^3 n factor.The additive decomposition separates univariate components from an infinitely smooth multivariate function, avoiding the curse of dimensionality.Source: §4 · p. 15

Generalized additive models

• n^−2β(γ∧1)/(2β(γ∧1)+1) + n^−2γ/(2γ+1) is also a lower bound for generalized additive models; when β = γ ≥2, the rate is n^−2β/(2β+1) [34].The model is a three-function composition with smooth component functions and an unknown link function.Source: §4 · pp. 15–16

Sparse tensor decomposition

• n^−2β/(2β+1) log^3 n + Δ(bfn, f0) is achieved for sparse tensor decompositions, with a rate exponent independent of input dimension d.Neural networks need only the sparse tensor representation; its univariate factors may be unknown and need not be orthogonal.Source: §4 · p. 16

Wavelet comparison

• Wavelet series estimators cannot exploit the underlying composition structure and can have much slower convergence rates in some setups.The composition assumption generalizes structural constraints such as additive models.Source: §4 · p. 17

5. Suboptimality of wavelet series estimators.

• n^−2α/(2α+1) log^3 n bounds the prediction error of deep neural-network reconstructions, whereas wavelet estimators cannot beat n^−2α/(2α+d) for the additive model f0(x)=h(x1+...+xd).Thus wavelet estimators can suffer the curse of dimensionality when d is large, while neural networks achieve faster rates.Source: §5 · p. 17
• The wavelet analysis considers uniform-design regression, empirical wavelet coefficients, truncation balancing squared bias and variance, and arbitrary coefficient index subsets.These ingredients support the lower-bound argument for wavelet series estimators.Source: §5 · pp. 17–18
• Theorem 4 establishes the wavelet lower bound R(bfn,f0) ≳ n^−2α/(2α+d) for compactly supported wavelets, arbitrary index sets, 0<α≤1, and Hölder radius K>0.The construction uses functions in the additive class with smoothness-controlled components.Source: §5 · pp. 18–19

6. A brief summary of related statistical theory for neural networks.

• The related-theory overview frames neural-network risk analysis around stochastic-error bounds from covering entropy and VC dimension and approximation theory for multivariate functions,.It first distinguishes shallow networks from multilayer architectures.Source: §6 · p. 19
• Shallow universal approximation guarantees arbitrary uniform approximation with enough hidden units, but rate results can be suboptimal because network representations require too many parameters.A cited rate is n^−2β/(2β+d+5), which loses a polynomial factor relative to the minimax rate.Source: §6 · p. 19
• Earlier multilayer results attain n^−2β/(2β+d) up to log n-factors for β≤1 with sigmoidal activations, extend to composition assumptions [40], but require activations at least as smooth as the signal.That smoothness requirement excludes ReLU activation.Source: §6 · p. 20
• For polynomial-growth activations, depth controls polynomial degree, but the quadratic-growth condition is crucial and the approach does not extend to ReLU, whose outputs are piecewise linear,.With σ(x)=x^2, L hidden layers generate polynomials of degree 2^L.Source: §6 · p. 20
• Recent ReLU approximation theory uses deep networks with few units that approximate the square function [72] [76], while this work proves good rates for every architecture satisfying Theorem 1 conditions.The approach also must handle bounded network parameters.Source: §6 · pp. 20–21

7.1. Embedding properties of network function classes.

• The network class supports enlargement, composition, depth synchronization, parallelization, and removal of inactive nodes, enabling construction and embedding of sparse ReLU networks.Composition joins compatible networks; identity layers synchronize depths, parallelization combines outputs, and inactive nodes can be removed without changing the computed function.Source: §7.1 · p. 21
• Approximation networks are built from multiplication modules and local Taylor expansions, then composed across normalized Hölder components and embedded into the target architecture.The proof uses m = ⌈log2 n⌉ and constructs component networks whose depth grows logarithmically with n before composition and embedding.Source: §7.1 · pp. 22–24

7.2. Proof of Theorem 2.

• Theorem 2 follows from a covering-entropy bound and an oracle inequality for bounded network classes under the nonparametric regression model.The argument applies the entropy result with δ = 1/n and assumes the regression function and candidate class are uniformly bounded, while design density bounds relate prediction risk to L2 distance.Source: §7.2 · pp. 25–26, §7.3 · p. 26

7.3. Proof of Theorem 3.

• Theorem 3 is proved by constructing separated regression functions from disjointly supported localized bumps and applying a Varshamov–Gilbert packing argument.The construction targets the composition component minimizing β_i + t_i, propagates the bump through the composition, and obtains pairwise L2 separation and controlled Kullback–Leibler divergence.Source: §7.3 · pp. 26–27

7.4. Proof of Lemma 1.

• The proof of Lemma 1 combines the minimax lower bound with Theorem 2 for empirical risk minimizers over bounded network classes.It sets the envelope to F = K + 1 and uses the lower-bound rate to control the relevant risk and complexity inequalities.Source: §7.4 · p. 28
• 7.4. Proof of Lemma 1 derives the result for constants depending on K, β, and d under a specified sparsity condition.The proof uses s ≤ c1ε^−d/β/(L log(1/ε)) and chooses c1 sufficiently small.Source: §7.4 · p. 28

7.5. Proofs for Section 5.

• 7.5. Proofs for Section 5 construct wavelet-based functions by selecting the first nonzero moment of ψ and treating the cases μ0 ≠ 0 and μ0 = 0 separately.The construction uses Lipschitz or Hölder-controlled functions and, in the μ0 = 0 case, the multinomial theorem to establish the lemma’s claim.Source: §7.5 · pp. 28–29
• Theorem 4 is proved by choosing an integer j* and applying Lemma 2 to construct a function of the form h(x1 + … + xd) with h in the stated Hölder class.The argument yields the displayed lower-bound scaling n^−2α/(2α+d) through the constructed wavelet coefficients.Source: §7.5 · pp. 29–30

1. Overparametrization and implicit regularization.

• The authors conjecture that implicit regularization is insufficient for consistency in overfitted additive-error regression and support this with a two-step argument involving shallow networks.The argument studies one-dimensional inputs with fixed first-layer parameters, then analyzes the resulting estimator under additive noise.Source: §1 · p. 31
• Under a weak overparameterization condition ensuring an interpolating shallow ReLU network exists, zero-initialized SGD with a suitable learning rate converges to a minimum-norm interpolant.The condition requires at least one knot in each interval determined by consecutive sample order statistics when m ≥ n.Source: §1 · p. 32
• The SGD limit is close to the natural cubic spline interpolant because its parameter norm corresponds to the L2 norm of the second derivative, which the natural spline minimizes.This links the optimization limit to a classical spline characterization.Source: §1 · p. 32
• The natural cubic spline interpolant is inconsistent with additive noise, so the corresponding SGD limit should also yield an inconsistent estimator.Noise-level deviations in the responses can produce a fixed lower bound on the spline’s average distance from the true regression function.Source: §1 · pp. 32–33
• Achieving rate-optimal nonparametric estimators therefore requires additional regularization that forces smoothing, with the article using network-weight sparsity as its complexity measure.The authors expect a related spline-interpolation phenomenon for deep networks, though not necessarily the natural cubic spline.Source: §1 · p. 33, §2 · p. 33

Pruning and topology

• Pruning trained fully connected networks can reduce parameters by about 90% without loss of efficiency, while the theory explicitly constructs a sparse topology attaining minimax rates up to log-factors.Fixing that topology and learning only nonzero weights may reduce flexibility to adapt to underlying data structure.Source: §2 · p. 33

Dynamic and structured sparsity

• Alternative sparse-network methods impose local sparsity constraints or evolve sparse topologies by removing weak connections, adding random connections, or retaining signal-channeling units.These approaches include predefined topologies, biological-network-inspired evolution, and connections between sparse networks and decision trees [61] [52] [1].Source: §2 · p. 34

Topology-only search

• Weight-agnostic neural networks [17] keep initialized weights fixed and learn only the topology through iterative candidate-model scoring and randomized generation around high-scoring models.No weight training is performed in this procedure.Source: §2 · p. 34

Extremely deep networks

• Extremely deep fully connected networks with polynomially increasing depth achieve optimal estimation rates without sparsity, but may require discontinuous weight assignments that raise stability concerns [KL, 77].The result concerns fixed-width networks and differs substantially from architectures commonly used in practice.Source: §2 · p. 34

Classification and nonparametric regression

• Well-separated classes and mostly correct labels can yield much faster convergence, making classification sample complexity lower than in nonparametric regression.With negligible mislabeling randomness, reconstruction resembles interpolation rather than denoising; this helps explain the gap observed for object recognition tasks [S].Source: §3 · p. 35
• Additive Gaussian perturbations increase the regression noise level but do not change the estimation rates, while deep-learning image classification is not robust to such perturbations [31].For the considered estimators, noise affects estimation-risk bounds through constants rather than rates.Source: §3 · p. 35
• The behavior of deep learning and reconstruction properties depend essentially on data structure, motivating estimation theory beyond nonparametric regression.The section identifies broader models as an open direction for future research.Source: §3 · p. 35
• Theorem 1 bounds the risk for any data-dependent estimator returning a sparsely connected neural network, rather than analyzing one fixed estimator [GMMM] and [S].The estimator-specific quantity Δ_n(bf_n, f_0) tightly controls estimation risk from above and below, unlike in zero-training-error interpolation settings where it is insufficient [S] and [78].Source: §4 · p. 35

4. Algorithms.

• Precise control of Δ_n(bf_n, f_0) remains unresolved because it may depend heavily on the learning rate, initialization, and energy landscape.The expectation defining Δ_n(bf_n, f_0) may include all randomness, including algorithmic randomization.Source: §4 · p. 35

5. High-dimensional input.

• For additive models, a parallelized network construction yields an upper-bound rate of d n^-2β/(2β+1) log^3 n when Δ_n(bf_n, f_0) is sufficiently small and d grows no faster than a power of n.Separate subnetworks approximate each component with O(N^-β) error using s ≍ dN log N parameters, producing total approximation error O(dN^-β).Source: §5 · p. 36
• High-dimensional inputs motivate studying networks with d = d_n ↑∞ and comparing their rates with other nonparametric procedures.The supplement also notes that basis expansions may require huge coefficients or exponentially many basis functions for functions generated by a single neuron [S].Source: §5 · p. 36

6. Function classes.

• Composition classes naturally reflect deep-network structure and include generalized additive models, while generic function spaces enable comparisons across estimation methods.The article reports slower wavelet convergence than sparsely connected deep ReLU networks for generalized additive models.Source: §6 · pp. 36–37
• More realistic network classes may need explicit weight-magnitude constraints reflecting SGD norm control and initialization neighborhoods rather than parameters bounded by one.Results transferring from ReLU to other activations may require parameters to become large, while training often keeps them near initialization [KL].Source: §6 · p. 37

7. Choice of the activation function.

• Deep networks can recover specialized constructions from the proof, but careful initialization may be necessary and reinitialization may help when optimization becomes trapped in a high-error local minimum.A simulation study found finite-sample deep-network performance near theoretical bounds, whereas another study reported difficulties with a smooth function using up to three hidden layers [GMMM].Source: §8 · p. 37

8. Real data.

• The real-data subsection contains acknowledgments for discussions about overfitting and stochastic gradient descent.No real-data results are supplied in the provided passages.Source: Acknowledgments · p. 37

SUPPLEMENT TO ”NONPARAMETRIC REGRESSION USING DEEP NEURAL NETWORKS WITH RELU ACTIVATION FUNCTION”

• The supplement provides additional material for the article “Nonparametric regression using deep neural networks with ReLU activation function.”The supplied passage identifies the supplement’s scope but reports no additional substantive result.Source: SUPPLEMENT TO ”NONPARAMETRIC… · p. 38

• The triangle-wave construction approximates x(1−x) exponentially fast while keeping all network parameters bounded by one and providing an explicit error bound.The resulting functions are piecewise linear interpolants of x(1−x) on dyadic grids, based on an inductive triangle-wave characterization.Source: Appendix A · pp. 38–39
• A network Mult_m with m+4 hidden layers maps [0,1]^2 into [0,1] and approximates xy with error at most 2^-m.The construction combines the x(1−x) approximation with a polarization identity and uses ReLU truncation to keep the output in [0,1].Source: Appendix A · pp. 40–41
• Pairwise composition of Mult_m yields an r-input product network with q(m+5) hidden layers, q=⌈log2(r)⌉, parameters bounded by one, and error at most r2^-m.The construction pads the inputs with ones and repeatedly multiplies neighboring entries in a binary tree; it also returns zero whenever any input component is zero.Source: Appendix A · pp. 41–42
• For monomials x^α of degree |α|<γ, parallelized multiplication networks approximate every monomial with depth proportional to (m+5)⌈log2(γ∨1)⌉ and widths scaling with the monomial count C_{r,γ}.For |α|≥2, the approximation error is |α|^2 2^-m; linear and constant monomials are represented exactly by a shallow network.Source: Appendix A · p. 42

Proof of Theorem 5

• The proof begins with local Taylor approximations and rewrites them as linear combinations of monomials around grid points.This follows the classical function-approximation strategy cited in [76].Source: Appendix B · p. 43
• A regular grid D(M) contains (M + 1)^r points, enabling localized approximation over [0, 1]^r.The grid points have coordinates ℓ_j/M for ℓ_j ∈ {0, ..., M}.Source: Appendix B · p. 43
• Lemma B.2 constructs localized hat-function products with support control and at most 49r^2(M + 1)^r(1 + (m + 5)⌈log2 r⌉) nonzero parameters.The network uses ReLU units for one-dimensional hat functions and multiplication subnetworks for their r-fold products.Source: Appendix B · pp. 44–45
• The theorem proof combines polynomial and hat-function subnetworks, multiplies corresponding outputs, sums over grid points, and adds a bounded-parameter scaling and shifting network.The resulting construction uses depth determined by m and log2(β ∨ r), while its active-parameter bounds are controlled using (M + 1)^r ≤ N.Source: Appendix B · pp. 45–46
• Choosing M so that (M + 1)^r ≤ N and using Lemma B.1 yields the claimed result because M^−β ≤ N^−β/r 3^β.The proof also uses (M + 2)^r ≤ (3M)^r to relate the grid resolution to N.Source: Appendix B · p. 47

Proof of Lemma 4

• Lemma 4 reduces the analysis to log N_n ≤ n, since both bounds are trivial when log N_n ≥ n, and then proves the result through four parts.The proof combines empirical-to-population risk control, estimator comparison, and empirical-risk-minimizer arguments.Source: Appendix C · pp. 47–48
• The first two parts relate population and empirical risks using a minimal δ-cover, concentration inequalities, Cauchy–Schwarz, Bernstein’s inequality, and a union bound.The resulting control uses covering complexity through log N_n and handles approximation errors introduced by the δ-cover.Source: Appendix C · pp. 47–49
• The final two parts convert the preceding inequalities into the stated lower and upper risk bounds, including the empirical-risk-minimizer comparison.The argument applies the algebraic inequality from part (I), the estimator bound from part (II), and concludes parts (III) and (IV).Source: Appendix C · pp. 50–51

Proof of Lemma 5

• Lemma 5 bounds the covering number of sparse neural-network functions by controlling Lipschitz propagation across layers and discretizing the nonzero parameters.The proof uses the product of layerwise Lipschitz constants, counts possible supports of s nonzero parameters, and applies a parameter grid of size δ/(2(L + 1)V).Source: Appendix C · pp. 51–52