Scalable Subsampling Inference for Deep Neural Networks

In the last several years, machine learning (ML) methods have been developed rapidly fueled by ever-increasing amounts of data and computational power. Among different ML methods, a popular and widely-used technique is Neural Networks (NN) that models the relationship between inputs and outputs through layers of connected computational neurons. The idea of applying such a biology-analogous framework can be traced to the work of McCulloch and Pitts, (1943).

At the end of the 20th century, people focused on the feed-forward Shallow Neural Networks (SNN) with sigmoid-type activation functions. An SNN has only one hidden layer but is shown to possess the universal approximation property, i.e., it can be used to approximate any Borel measurable function from one finite dimensional space to another with any desired degree of accuracy—see Cybenko, (1989); Hornik et al., (1989) and references within. However, the SNN practical performance left much to be desired. In the last ten or so years, Deep Neural Networks (DNN) received increased attention due to their great empirical performance.

Although DNN have become a state-of-the-art model, their theoretical foundation is still in development. Notably, Yarotsky, (2018); Yarotsky and Zhevnerchuk, (2020) explored the approximation ability of DNN for any function $f$ that belongs to a Hölder Banach space; here, the sigmoid-type activation functions are now replaced by ReLU-type functions to avoid the gradient vanishing problem. The aforementioned work showed that the optimal error of the DNN estimator $f_{\text{DNN}}$ can be uniformly bounded, i.e.,

here, $\xi$ is some smoothness measurement of the target function $f:{\bf R}^{d}\to{\bf R}$ —see Section 4 for a formal definition; $W$ is the size of a neural network $f_{\text{DNN}}$, i.e., the total number of parameters; and $d$ is the dimension of the function inputs.

However, the bound (1) is not useful in practice. The reason is three-fold: (a) it requires a discontinuous weight assignment to build the desired DNN, so it is not feasible to train such DNN with usual gradient-based methods; (b) the structure of the DNN might not be the standard fully connected form so finding the satisfied specific structure becomes another difficult; most importantly, (c) this error bound is on the optimal estimation we can achieve from a finely designed DNN. It fails to tell us any story about the situation of applying the DNN estimator to solving real-world problems.

For example, what is the performance of the DNN to estimate a regression function with $n$ independent samples $\{(Y,\bm{X}_{i})\}_{i=1}^{n}$ generated from an underlying true model $f$? It is easy to see that the error $\varepsilon$ of $f_{\text{DNN}}$ in sup-norm can be arbitrarily small if we allow $W$ to be arbitrarily large based on Eq. 1. However, this optimal performance is hardly achievable and only represents the theoretically best estimation. What we attempt to do in this paper is to determine an empirically optimal $\widehat{f}_{\text{DNN}}$ with samples $\{(Y,\bm{X}_{i})\}_{i=1}^{n}$ and then explore its estimation and prediction inference. Guided by this spirit, people usually think $\widehat{f}_{\text{DNN}}$ as an $M$-estimator and set different loss functions for various purposes:

here $\mathscr{F}_{\text{DNN}}$ is a user-chosen space that contains all DNN candidates; $L(\cdot,\cdot)$ is the loss function, e.g., Mean Squared Errors loss for the regression problem with real-valued output, i.e., $L(f_{\theta}(\bm{x}_{i}),y_{i})=(f_{\theta}(\bm{x}_{i})-y_{i})^{2}/2$; $\{(y,\bm{x}_{i})\}_{i=1}^{n}$ are realizations of $\{(Y,\bm{X}_{i})\}_{i=1}^{n}$.

In the paper at hand, we consider DNN-based estimation and prediction inference in the data-generating model: $Y_{i}=f(\bm{X}_{i})+\epsilon_{i}$; here, the $\epsilon_{i}$ are independent, identically distributed (i.i.d.) from a distribution $F_{\epsilon}$ that has mean 0 and variance $\sigma^{2}$—we will use the shorthand $\epsilon_{i}\sim$ i.i.d. $(0,\sigma^{2}$). Consequently, $f(\bm{x}_{i})=\mathbb{E}(Y_{i}|\bm{X}_{i}=\bm{x}_{i})$. Furthermore, the regression function $f(\cdot)$ will be assumed to satisfy some smoothing condition which will be specified later. Note that the additive model with heteroscedastic error: $Y_{i}=f(\bm{X}_{i})+g(\bm{X}_{i})\cdot\epsilon_{i}$ can be analyzed similarly by applying two DNNs, one to estimate $f(\cdot)$ and one for $g(\cdot)$.

From a nonparametric regression view, it is well-known that the optimal convergence rate of the estimation for a $p$-times continuously differentiable regression function of a $d$-dimensional argument is $n^{2p/(2p+d)}$—see Stone, (1982). If we assume the regression function belongs to a more general Hölder Banach space, we can define a non-integer $\xi=p+s$ to represent the smoothness order of $f$; here $0<s\leq 1$ is the Hölder coefficient. The optimal rate of non-parametric estimation can also be extended to such non-integer smoothness order; see Condition $3^{{}^{\prime}}$ and Definition 2 of Kohler et al., (2023). Focusing on DNN estimation, the optimal and achievable error bound on the $L_{2}$ norm of $\widehat{f}_{\text{DNN}}$ is $O(n^{-\xi/(\xi+d)}\cdot\log^{8}(n))$ with a high probability; this bound is due to Farrell et al., (2021) but the rate appears slower than the optimal rate that we can attain. Besides, although $\widehat{f}_{\text{DNN}}$ will become more accurate as the sample size increases, training DNN becomes very time-consuming. Moreover, it is infeasible to load massive data into a PC or even a supercomputer since its node memory is also limited in the computation process as pointed out by Zou et al., (2021).

In this paper, we first give a small improvement on the bound of Farrell et al., (2021) using the latest results on the DNN approximation ability. Then, we resolve the computational issue involving huge data by applying the Scalable Subsampling technique of Politis, (2021) to create a set of subsamples and then build a so-called subagging DNN estimator $\overline{f}_{\text{DNN}}$. Under regularity conditions, we can show that the subagging DNN estimator $\overline{f}_{\text{DNN}}$ could possess a faster convergence rate than a single DNN estimator $\widehat{f}_{\text{DNN}}$ trained on the whole sample. Lastly, using the same set of subsamples, we can build a Confidence Interval (CI) for $f$ based on $\overline{f}_{\text{DNN}}$. Due to the prevalent undercoverage phenomenon of CIs with finite samples, we propose two ideas to improve the empirical coverage rate: (1) we enlarge the CI by maximizing the margin of errors in an appropriate way; (2) we take an iterated subsampling method to build a specifically designed CI which is a combination of pivot-CI and quantile-CI. Beyond estimation inference, we also perform predictive inference (with both point and interval predictions).

The paper is organized as follows. In Section 2, we give a short introduction to the structure of DNN. In Section 3, we describe the methodology of scalable subsampling. Subsequently, the performance of the subagging DNN estimator and its associated confidence/prediction intervals are analyzed in Section 4 and Section 5. Simulations and conclusions are provided in Section 6 and Section 7, respectively. Proofs are given in Appendix: A and additional simulations studies are presented in Appendix: B and Appendix: C.

In terms of notation, we will use the following norms: $\|g\|_{L_{2}(\bm{X})}:=\mathbb{E}[g(\bm{X})^{2}]^{1/2}$; $\|g\|_{\infty}:=\sup_{\bm{x}}|g(\bm{x})|$. Also, we employ the notation $a_{n}=\Theta\left(d_{n}\right)$ to denote “exact order”, i.e., that there exist two constants $c_{1},c_{2}$ satisfying $c_{1}\cdot c_{2}>0$, and $c_{1}d_{n}\leq a_{n}\leq c_{2}d_{n}$. We also use $\mathbb{E}_{n}[\cdot]$ to represent the sample average operator.

For completeness, we now give a brief introduction to the fully connected forward DNN with ReLU activation functions. Hereafter, we refer to DNN as the standard fully connected deep neural network with the ReLU activation function. In short, the DNN can be viewed as a parameterized family of functions. Its structure mainly depends on the input dimension $d$, depth $L\in\mathbb{N}$, width $\bm{H}\in\mathbb{N}^{L}$ and the output dimension. The depth $L$ describes how many hidden layers a DNN possesses; the width $\bm{H}=(H_{1},\ldots,H_{L})$ represents the number of neurons in each hidden layer. The forward property implies that the input, hidden neurons and output are connected in an acyclic relationship. The fully connected property indicates that each hidden neuron receives information from all hidden neurons at the previously hidden layer in a functional way.

Formally, if we let $\bm{u}_{l}=(u_{l,1},\ldots,u_{l,H_{l}})^{T}$ to represent all number of neurons at the $l$-th hidden layer for $l=0,\ldots,L+1$; here, $\bm{u}_{0}$ represents the input vector $(x_{1},\ldots,x_{d})^{T}$ and $\bm{u}_{L+1}$ is the output. Therefore, we can pretend that the input layer and the output layer are the $0$-th and $(L+1)$-th hidden layers, respectively. Then, $u_{l,i}=\sigma(\bm{u}_{l-1}^{T}\bm{w}_{l-1,i}+b_{l-1,i})$ for $l=1,\ldots,L$ and $i=1,\ldots,H_{l}$; here $\bm{w}_{l-1,i}\in\mathbb{R}^{H_{l-1}}$ is the weight vector which connects the $(l-1)$-th hidden layer and the neuron $u_{l,i}$; $b_{l-1,i}\in\mathbb{R}$ is the corresponding intercept term; $\sigma(\cdot)$ is the so-called activation function and we take the ReLU function in this paper. To get the output layer, we just take $u_{L+1,i}=\bm{u}_{L}^{T}\bm{w}_{L,i}+b_{L,i}$ for $i=1,\ldots,H_{L+1}$; here $H_{L+1}$ is equal to the output dimension. To express the functionality of the DNN in a more concise way, we can stack $\{\bm{w}_{l-1,i}^{T}\}_{i=1}^{H_{l}}$ by row to get $\bm{W}_{l-1}\in\mathbb{R}^{H_{l}}\times\mathbb{R}^{H_{l-1}}$ and collect $\{b_{l-1,i}\}_{i=1}^{H_{l}}$ to be a vector $\bm{b}_{l-1}$ for $l=1,\ldots,L+1$. Subsequently, we can treat the DNN as a function that takes the input $\bm{x}$ and returns output in the below way:

We can understand that the function $f_{\text{DNN}}(\bm{x})$ maps $\bm{x}$ to the $1$-st hidden layer and then map the $1$-st hidden to the $2$-nd hidden layer and so on iteratively with weights $\{\bm{W}_{l}\}_{l=0}^{L}$, $\{\bm{b}_{l}\}_{l=0}^{L}$ and the activation function $\sigma(\cdot)$. We can then compute the total number of parameters in a DNN by the formula $W=\sum_{i=0}^{L}(H_{i}\cdot H_{i+1}+H_{i+1})$. A simple DNN is presented in Fig. 1. It has a constant width of 4 and a depth of 2.

Scalable subsampling is one type of non-stochastic subsampling technique proposed by Politis, (2021). Assume that we observe the sample $\{\bm{Z}_{1},\ldots,\bm{Z}_{n}\}$; then, scalable subsampling relies on $q=\lfloor(n-b)/h\rfloor+1$ number of subsamples $B_{1},\ldots,B_{q}$ where $B_{j}=\{\bm{Z}_{(j-1)h+1},\ldots,$ $\bm{Z}_{(j-1)h+b}\}$; here, $\lfloor\cdot\rfloor$ denotes the floor function, and $h$ controls the amount of overlap (or separation) between $B_{j}$ and $B_{j+1}$. In general, the subsample size $b$ and the overlap $h$ are functions of $n$, but these dependencies will not be explicitly denoted, e.g.,

where $0<\beta<1$ and $a>0$. More importantly, tuning $b$ and $h$ can make scalable subsampling samples have different properties. For example, if $h=1$, the overlap is the maximum possible; if $h=0.2b$, there is 80$\%$ overlap between $B_{j}$ and $B_{j+1}$; if $h=b$, there is no overlap between $B_{j}$ and $B_{j+1}$ but these two blocks are adjacent; if $h=1.2b$, there is a block of about $0.2b$ data points that separate the blocks $B_{j}$ and $B_{j+1}$.

The bagging idea was initially proposed by Breiman, (1996), where the subsample is bootstrapped (sampling with replacement) with the same size as the original sample. As revealed by that work, the main benefit of taking this technique is that the mean-squared error (MSE) of the bagging estimator can decrease, especially for unstable estimators that may change a lot with different samples, e.g., neural networks and regression trees. There are ample works about combing the neural networks with the bagging technique to improve its generalization performance; e.g., see applications in the work of Ha et al., (2005); Khwaja et al., (2015) for references. However, the drawback of the original bagging method is that the estimation process needs to be performed with $n$-size bootstrap resamples many times which is infeasible with massive data. Bühlmann and Yu, (2002) proposed the subagging idea which is based on all subsamples as opposed to bootstrap resamples. However, even choosing a single random subsample could be computationally challenging when $n$ is large. As pointed out in Ting, (2021), drawing a random sample of size $b$ from $n$ items using the Sparse Fisher-Yates Sampler takes $O(b)$ time and space which corresponds to optimal time and space complexity for this problem.

Facing such computational dilemmas, scalable subsampling and subagging as proposed by Politis, (2021) can be seen as an extension of the Divide-and-Conquer principle—see e.g. Jordan, (2013). Moreover, in addition to the computational savings, scalable subagging may yield an estimator that is not less (and sometimes more) accurate than the original; the following example illustrates such a case.

In the next section, we will introduce how to compute the scalable subsampling DNN estimator. Then, we will show that our aggregated DNN estimator could possess a smaller MSE than the optimal DNN estimator trained on the whole sample, under some conditions. We also discuss some specifically designed confidence intervals to measure the estimation accuracy via the approaches mentioned in Section 1.

Although the DNN has captured much attention in practice, its theoretical validation is still in development. Recently, Farrell et al., (2021) gave a non-asymptotic error bound to measure the performance of the DNN estimator under two regularity assumptions. To sync with the latest results on the approximation ability of DNNs, we replace their second assumption (Assumption 2 in their work) with the smooth function condition assumed by Yarotsky and Zhevnerchuk, (2020). We present these new assumptions below:

• •

  A1: The regression data are i.i.d. copies of $\bm{Z}=(Y,\bm{X})\in$ $\mathscr{Y}\times[-1,1]^{d}$, where $\bm{X}$ has a continuous distribution, and $\mathscr{Y}\subset[-M,M]$ for some positive constant $M$. Correspondingly, we set the space of all DNN candidate functions to be $\mathscr{F}_{\text{DNN}}=\{f_{\theta}:||f_{\theta}||_{\infty}\leq 2M\}$.

• •

  A2: The target regression function $f$ lies in the Hölder Banach space $\mathscr{C}^{k,\alpha}\left([-1,1]^{d}\right)$ which is the space of $k$ times continuously differentiable functions on $[-1,1]^{d}$ having a finite norm defined by

  $$\|f\|_{\mathscr{C}^{k,\alpha}\left([-1,1]^{d}\right)}=\max\left\{\max_{\mathbf{k}:|\mathbf{k}|\leq k}\max_{\mathbf{x}\in[-1,1]^{d}}\left|D^{\mathbf{k}}f(\mathbf{x})\right|,\max_{\mathbf{k}:|\mathbf{k}|=k}\sup_{\begin{subarray}{c}\mathbf{x},\mathbf{y}\in[-1,1]^{d}\\ \mathbf{x}\neq\mathbf{y}\end{subarray}}\frac{\left|D^{\mathbf{k}}f(\mathbf{x})-D^{\mathbf{k}}f(\mathbf{y})\right|}{\|\mathbf{x}-\mathbf{y}\|^{\alpha}}\right\},$$

  where the smoothness index is $\xi=k+\alpha$ with an integer $k\geq 0$ and $0<\alpha\leq 1$.

• •

  A3: The sample size $n$ is larger than $(2eM)^{2}\vee\text{Pdim}(f_{\text{DNN}})$ where $\text{Pdim}(f_{\text{DNN}})$ is the pseudo-dimension of $f_{\text{DNN}}$ which satisfies:

  $$c\cdot WL\log(W/L)\leq\operatorname{Pdim}(f_{\text{DNN}})\leq C\cdot WL\log W,$$

  with some universal constants $c,C>0$ and Euler’s number $e$; see Bartlett et al., (2019) for details.

As shown in Farrell et al., (2021), with $H_{1}=H_{2}=\cdots=H_{L}=\Theta(n^{\frac{d}{2(\xi+\alpha)}}\log^{2}n)$, $L=\Theta(\log n)$, the $L_{2}$ norm loss and empirical mean squared error of the deep fully connected ReLU network estimator from Eq. 2 can be bounded with probability at least $1-\exp\left(-n^{\frac{d}{\xi+d}}\log^{8}n\right)$, i.e.,

here $C_{1}>0$ is a constant which is independent of $n$ and may depend on $d,M$, and other fixed constants.

Obviously, the $L_{2}$ norm error bound in (3) is sub-optimal compared to the fastest convergence rate we can achieve for nonparametric function estimation. With the latest approximation theory on DNNs, we can improve the error bound in Eq. 3 by decreasing the power of the $\log(n)$ term. Meanwhile, this faster rate is satisfied with a narrower DNN. We give our first theorem about the convergence rate of $\widehat{f}_{\text{DNN}}$ below.