Provable Identifiability of Two-Layer ReLU Neural Networks via LASSO
  Regularization
LASSO regularization is a popular regression tool to enhance the prediction
accuracy of statistical models by performing variable selection through the
$\ell_1$ penalty, initially formulated for the linear model and its variants.
In this paper, the territory of LASSO is extended to two-layer ReLU neural
networks, a fashionable and powerful nonlinear regression model. Specifically,
given a neural network whose output $y$ depends only on a small subset of input
$\boldsymbol{x}$, denoted by $\mathcal{S}^{\star}$, we prove that the LASSO
estimator can stably reconstruct the neural network and identify
$\mathcal{S}^{\star}$ when the number of samples scales logarithmically with
the input dimension. This challenging regime has been well understood for
linear models while barely studied for neural networks. Our theory lies in an
extended Restricted Isometry Property (RIP)-based analysis framework for
two-layer ReLU neural networks, which may be of independent interest to other
LASSO or neural network settings. Based on the result, we advocate a neural
network-based variable selection method. Experiments on simulated and
real-world datasets show promising performance of the variable selection
approach compared with existing techniques.
Gen Li, Ganghua Wang, Jie Ding
# Provable Identifiability of Two-Layer ReLU Neural Networks via LASSO Regularization

###### Abstract

LASSO regularization is a popular regression tool to enhance the prediction accuracy of statistical models by performing variable selection through the \(\ell_{1}\) penalty, initially formulated for the linear model and its variants. In this paper, the territory of LASSO is extended to two-layer ReLU neural networks, a fashionable and powerful nonlinear regression model. Specifically, given a neural network whose output \(y\) depends only on a small subset of input \(x\), denoted by \(\mathcal{S}^{*}\), we prove that the LASSO estimator can stably reconstruct the neural network and identify \(\mathcal{S}^{*}\) when the number of samples scales logarithmically with the input dimension. This challenging regime has been well understood for linear models while barely studied for neural networks. Our theory lies in an extended Restricted Isometry Property (RIP)-based analysis framework for two-layer ReLU neural networks, which may be of independent interest to other LASSO or neural network settings. Based on the result, we advocate a neural network-based variable selection method. Experiments on simulated and real-world datasets show promising performance of the variable selection approach compared with existing techniques.

 Lasso, Identifiability, Neural network, Nonlinear regression, Variable selection.

## I Introduction

Given \(n\) observations \((y_{i},\mathbf{x}_{i})\), \(i=1,\ldots,n\), we often model them with the regression form of \(y_{i}=f(\mathbf{x}_{i})+\xi_{i}\), with an unknown function \(f\), \(\mathbf{x}_{i}\in\mathbb{R}^{p}\) being the input variables, and \(\xi_{i}\) representing statistical errors. A general goal is to estimate a regression function \(\widehat{f}_{n}\) close to \(f\) for prediction or interpretation. This is a challenging problem when the input dimension \(p\) is comparable or even much larger than the data size \(n\). For linear regressions, namely \(f(\mathbf{x})=\mathbf{w}^{\top}\mathbf{x}\), the least absolute shrinkage and selection operator (LASSO) [1] regularization has been established as a standard tool to estimate \(f\). The LASSO has also been successfully used and studied in many nonlinear models such as generalized linear models [2], proportional hazards models [3], and neural networks [4]. For LASSO-regularized neural networks, existing works have studied different properties, such as convergence of training [5], model pruning [6, 7], and feature selection [8, 9]. The LASSO regularization has also been added into the standard deep learning toolbox of many open-source libraries, e.g., Tensorflow [10] and Pytorch [11].

Despite the practical success of LASSO in improving the generalizability and sparsification of neural networks, whether one can use LASSO for identifying significant variables is underexplored. For linear models, the variable selection problem is also known as support recovery or feature selection in different literature. Selection consistency requires that the probability of \(\text{supp}(\widehat{\mathbf{w}})=\text{supp}(\mathbf{w})\) converges to one as \(n\to\infty\). The standard approach to selecting a parsimonious sub-model is to either solve a penalized regression problem or iteratively pick up significant variables [12]. The existing methods differ in how they incorporate unique domain knowledge (e.g., sparsity, multicollinearity, group behavior) or what desired properties (e.g., consistency in coefficient estimation, consistency in variable selection) to achieve [13]. For instance, consistency of the LASSO method [1] in estimating the significant variables has been extensively studied under various technical conditions, including sparsity, mutual coherence [14], restricted isometry [15], irrepresentable condition [16], and restricted eigenvalue [17].

Many theoretical studies of neural networks have focused on the generalizability. For example, a universal approximation theorem was established that shows any continuous multivariate function can be represented precisely by a polynomialized two-layer network [18]. It was later shown that any continuous function could be approximated arbitrarily well by a two-layer perceptron with sigmoid activation functions [19], and an approximation error bound of using two-layer neural networks to fit arbitrary smooth functions has been established [20, 21]. Statistically, generalization error bounds for two-layer neural networks [21, 22] and multi-layer networks [23, 24, 25] have been developed. From an optimization perspective, the parameter estimation of neural networks was cast into a tensor decomposition problem where a provably global optimum can be obtained [26, 27, 28]. Very recently, a dimension-free Rademacher complexity to bound the generalization error for deep ReLU neural networks was developed to avoid the curse of dimensionality [29]. It was proved that certain deep neural networks with few nonzero network parameters could achieve minimax rates of convergence [30]. A tight error bound free from the input dimension was developed by assuming that the data is generated from a generalized hierarchical interaction model [31].

This work theoretically studies the identifiability of neural networks and uses it for variable selection. Specifically, suppose data observations are generated from a neural network with only a few nonzero coefficients. The identifiability concerns the possibility of identifying those coefficients, which may be further used to identify a sparse set of input variables that are genuinely relevant to the response. In this direction, LASSO and its variant Group LASSO [32] have been advocated to regularize neural-network for variable selection in
practice (see, e.g., [6, 8, 9, 33]).

In this paper, we consider the following class of two-layer ReLU neural networks for regression.

\[\mathcal{F}_{r}=\bigg{\{}f:\boldsymbol{x}\mapsto f(\boldsymbol{x})=\sum_{j=1}^{ r}a_{j}\mathrm{relu}(\boldsymbol{w}_{j}^{\top}\boldsymbol{x}+b_{j}),\]

\[\text{where }a_{j},b_{j}\in\mathbb{R},\boldsymbol{w}_{j}\in\mathbb{R}^{p} \bigg{\}}.\]

Here, \(p\) and \(r\) denote the input dimension and the number of neurons, respectively. We will assume that data are generated from a regression function in \(\mathcal{F}_{r}\) perturbed by a small term. We will study the following two questions.

First, if the underlying regression function \(f\) admits a parsimonious representation so that only a small set of input variables, \(\mathcal{S}^{\star}\), is relevant, can we identify them with high probability given possibly noisy measurements \((y_{i},\boldsymbol{x}_{i})\), for \(i=1,\ldots,n\)? Second, is such an \(\mathcal{S}^{\star}\) estimable, meaning that it can be solved from an optimization problem with high probability, even in small-\(n\) and large-\(p\) regimes?

To address the above questions, we will establish a theory for neural networks with the LASSO regularization by considering the problem \(\min_{\boldsymbol{W},\boldsymbol{a},\boldsymbol{b}}\|\boldsymbol{W}\|_{1}\) under the constraint of

\[\frac{1}{n}\sum_{i=1}^{n}\biggl{(}y_{i}-\sum_{j=1}^{r}a_{j}\mathrm{relu}( \boldsymbol{w}_{j}^{\top}\boldsymbol{x}_{i}+b_{j})\biggr{)}^{2}\leq\sigma^{2},\]

which is an alternative version of the \(\ell_{1}\)-regularization. More notational details will be introduced in Subsection II-B.

Our theory gives positive answers to the above questions. We theoretically show that the LASSO-type estimator can stably identify ReLU neural networks with sparse input signals, up to a permutation of hidden neurons. We only focus on the varying \(n\) and \(p\) and implicitly assume that the sparsity of \(\boldsymbol{W}^{\star}\) and the number of neurons \(r\) are fixed. While this does not address wide neural networks, we think it still corresponds to a nontrivial and interesting function class. For example, the class contains linear functions when input signals are bounded. Our result is rather general as it applies to noisy observations of \(y\) and dimension regimes where the sample size \(n\) is much smaller than the number of input variables \(p\). The theory was derived based on new concentration bounds and function analysis that may be interesting in their own right.

Inspired by the developed theory, we also propose a neural network-based variable selection method. The idea is to use the neural system as a vehicle to model nonlinearity and extract significant variables. Through various experimental studies, we show encouraging performance of the technique in identifying a sparse set of significant variables from large dimensional data, even if the underlying data are not generated from a neural network. Compared with popular approaches based on tree ensembles and linear-LASSO, the developed method is suitable for variable selection from nonlinear, large-dimensional, and low-noise systems.

The rest of the paper is outlined as follows. Section II introduces the main theoretical result and proposes an algorithm to perform variable selection. Section III uses simulated and real-world datasets to demonstrate the proposed theory and algorithm. Section IV concludes the paper.

## II Main result

### _Notation_

Let \(\boldsymbol{u}_{\mathcal{S}}\) denote the vector whose entries indexed in the set \(\mathcal{S}\) remain the same as those in \(\boldsymbol{u}\), and the remaining entries are zero. For any matrix \(\boldsymbol{W}\in\mathbb{R}^{p\times r}\), we define

\[\|\boldsymbol{W}\|_{1}=\sum_{1\leq k\leq p,1\leq j\leq r}|w_{kj}|,\,\| \boldsymbol{W}\|_{\mathrm{F}}=\biggl{(}\sum_{1\leq k\leq p,1\leq j\leq r}w_{ kj}^{2}\biggr{)}^{1/2}.\]

Similar notations apply to vectors. The inner product of two vectors is denoted as \(\langle\boldsymbol{u},\boldsymbol{v}\rangle\). Let \(\boldsymbol{w}_{j}\) denote the \(j\)-th column of \(\boldsymbol{W}\). The sparsity of a matrix \(\boldsymbol{W}\) refers to the number of nonzero entries in \(\boldsymbol{W}\). Let \(\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p})\) denote the standard \(p\)-dimensional Gaussian distribution, and \(1(\cdot)\) denote the indicator function. The rectified linear unit (ReLU) function is defined by \(\mathrm{relu}(v)=\max\{v,0\}\) for all \(v\in\mathbb{R}\).

### _Formulation_

Consider \(n\) independently and identically distributed (i.i.d.) observations \(\{\boldsymbol{x}_{i},y_{i}\}_{1\leq i\leq n}\) satisfying

\[y_{i}=\sum_{j=1}^{r}a_{j}^{\star}\cdot\mathrm{relu}(\boldsymbol{w}_{j}^{\star \top}\boldsymbol{x}_{i}+b_{j}^{\star})+\xi_{i} \tag{1}\]

with \(\boldsymbol{x}_{i}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p})\), where \(r\) is the number of neurons, \(a_{j}^{\star}\in\{1,-1\}\), \(\boldsymbol{w}_{j}^{\star}\in\mathbb{R}^{p}\), \(b_{j}^{\star}\in\mathbb{R}\), and \(\xi_{i}\) denotes the random noise or approximation error obeying

\[\frac{1}{n}\sum_{i=1}^{n}\xi_{i}^{2}\leq\sigma^{2}. \tag{2}\]

In the above formulation, the assumption \(a_{j}^{\star}\in\{1,-1\}\) does not lose generality since \(a\cdot\mathrm{relu}(b)=ac\cdot\mathrm{relu}(b/c)\) for any \(c>0\). The setting of Inequality (2) is for simplicity. If \(\xi_{i}\)'s are unbounded random variables, the theoretical result to be introduced will still hold, with more explanations in the Appendix. The \(\xi_{i}\)'s are not necessarily i.i.d., and \(\sigma\) is allowed to be zero, which reduces to the noiseless scenario.

Let \(\boldsymbol{W}^{\star}=[\boldsymbol{w}_{1}^{\star},\ldots,\boldsymbol{w}_{r}^{ \star}]\in\mathbb{R}^{p\times r}\) denote the data-generating coefficients. The main problem to address is whether one can stably identify those nonzero elements, given that most entries in \(\boldsymbol{W}^{\star}\) are zero. The study of neural networks from an identifiability perspective is essential. Unlike the generalizability problem that studies the predictive performance of machine learning models, the identifiability may be used to interpret modeling results and help scientists make trustworthy decisions. To illustrate this point, we will propose to use neural networks for variable selection in Subsection II-E.

To answer the above questions, we propose to study the following LASSO-type optimization. Let \(\big{(}\widehat{\boldsymbol{W}},\widehat{\boldsymbol{a}},\widehat{ \boldsymbol{b}}\big{)}\) be a solution to the following optimization problem,

\[\min_{\boldsymbol{W},\boldsymbol{a},\boldsymbol{b}}\|\boldsymbol{W} \|_{1} \tag{3}\] \[\text{subject to }\frac{1}{n}\sum_{i=1}^{n}\biggl{(}y_{i}-\sum_{j=1}^{r}a_ {j}\cdot\mathrm{relu}(\boldsymbol{w}_{j}^{\top}\boldsymbol{x}_{i}+b_{j}) \biggr{)}^{2}\leq\sigma^{2},\]

within the feasible range \(\boldsymbol{a}\in\{1,-1\}^{r}\), \(\boldsymbol{W}\in\mathbb{R}^{p\times r}\), and \(\boldsymbol{b}\in\mathbb{R}^{r}\). Intuitively, the optimization operates under the
constraint that the training error is not too large, and the objective function tends to sparsify \(\mathbf{W}\). Under some regularity conditions, we will prove that the solution is indeed sparse and close to the data-generating process.

### _Main result_

We make the following technical assumption.

**Assumption 1**.: _For some constant \(B\geq 1\),_

\[1\leq\|\mathbf{w}_{j}^{\star}\|_{2}\leq B\quad\text{and}\quad|b_{j}^{\star}|\leq B \quad\forall 1\leq j\leq r. \tag{4}\]

_In addition, for some constant \(\omega>0\),_

\[\max_{j,k=1,\dots,r,j\neq k}\frac{\left|\langle\mathbf{w}_{j}^{\star},\mathbf{w}_{k}^ {\star}\rangle\right|}{\|\mathbf{w}_{j}^{\star}\|_{2}\|\mathbf{w}_{k}^{\star}\|_{2}}\leq \frac{1}{r^{\omega}}. \tag{5}\]

**Remark 1** (Discussion of Assumption 1).: _The condition in (4) is a normalization only for technical convenience, since we can re-scale \(\mathbf{w}_{j},b_{j},y_{i},\sigma\) proportionally without loss of generality. Though this condition implicitly requires \(\mathbf{w}_{j}^{\star}\neq\mathbf{0}\) for all \(j=1,\dots,r\), it is reasonable since it means the neuron \(j\) is not used/activated. The condition in (5) requires that the angle of any two different coefficient vectors is not too small. This condition is analogous to a bounded-eigenvalue condition often assumed for linear regression problems, where each \(w_{j}^{\star}\) is understood as a column in the design matrix. This condition is by no means mild or easy to verify in practice. Nevertheless, as our focused regime is large \(p,n\) but small \(r\), we think the condition in (5) is still reasonable. For example, when \(r=2\), this condition simply requires \(w_{1}^{\star}\neq w_{2}^{\star}\)._

Our main result shows that if \(\mathbf{W}^{\star}\) is sparse, one can stably reconstruct a neural network when the number of samples (\(n\)) scales logarithmically with the input dimension (\(p\)). A skeptical reader may ask how the constants exactly depend on the sparsity and \(r\). We will provide a more elaborated result in Subsection II-D and introduce the proof there.

**Theorem 1**.: _Under Assumption 1, there exist some constants \(c_{1},c_{2},c_{3}>0\) depending only (polynomially) on the sparsity of \(\mathbf{W}^{\star}\) such that for any \(\delta>0\), one has with probability at least \(1-\delta\),_

\[\widehat{\mathbf{a}}=\mathbf{\Pi}\mathbf{a}^{\star}\text{ and }\|\widehat{\mathbf{W}}-\mathbf{W}^{ \star}\mathbf{\Pi}^{\top}\|_{\mathrm{F}}+\|\widehat{\mathbf{b}}-\mathbf{\Pi}\mathbf{b}^{\star }\|_{2}\leq c_{1}\sigma \tag{6}\]

_for some permutation matrix \(\mathbf{\Pi}\), provided that_

\[n>c_{2}\log^{4}\frac{p}{\delta}\qquad\text{and}\qquad\sigma<c_{3}. \tag{7}\]

**Remark 2** (Interpretation of Theorem 1).: _The permutation matrix \(\mathbf{\Pi}\) is necessary since the considered neural networks produce identical predictive distributions (of \(y\) conditional \(\mathbf{x}\)) under any permutation of the hidden neurons. The result states that the underlying neural coefficients can be stably estimated even when the sample size \(n\) is much smaller than the number of variables \(p\). Also, the estimation error bound is at the order of \(\sigma\), the specified noise level in (2)._

_Suppose that we define the signal-to-noise ratio (SNR) to be \(\mathbb{E}\|\mathbf{x}\|^{2}/\sigma^{2}\). An alternative way to interpret the theorem is that a large SNR ensures the global minimizer to be close to the ground truth with high probability. One may wonder what if the \(\sigma<c_{3}\) condition is not met. We note that if \(\sigma\) is too large, the error bound in (6) would be loose, and it is not of much interest anyway. In other words, if the SNR is small, we may not be able to estimate parameters stably. This point will be demonstrated by experimental studies in Section III._

The estimation results in Theorem 1 can be translated into variable selection results as shown in the following Corollary 1. The connection is based on the fact that if \(i\)-th variable is redundant, the underlying coefficients associated with it should be zero. Let \(\mathbf{w}_{i,\cdot}^{\star}\) denote the \(i\)-th row of \(\mathbf{W}^{\star}\). Then,

\[\mathcal{S}^{\star}=\{1\leq i\leq p:\|\mathbf{w}_{i,\cdot}^{\star}\|_{2}>0\}\]

characterizes the "significant variables." Corollary 1 states that the set of variables with non-vanished coefficient estimates contains all the significant variables. The corollary also shows that with a suitable shrinkage of the coefficient estimates, one can achieve variable selection consistency.

**Corollary 1** (Variable selection).: _Let \(\widehat{\mathcal{S}}_{c_{1}\sigma}\subseteq\{1,\dots,p\}\) denote the sets of \(i\)'s such that \(\|\widehat{\mathbf{w}}_{i,\cdot}\|_{2}>c_{1}\sigma\). Under the same assumption as in Theorem 1, and \(\min_{1\leq i\leq r}\|\mathbf{w}_{i,\cdot}^{\star}\|_{2}>2c_{1}\sigma\), for any \(\delta>0\), one has_

\[\mathbb{P}(\mathcal{S}^{\star}=\widehat{\mathcal{S}}_{c_{1}\sigma})\geq 1-\delta,\]

_provided that \(n>c_{2}\log^{4}\frac{p}{\delta}\) and \(\sigma<c_{3}\)._

Considering the noiseless scenario \(\sigma=0\), Theorem 1 also implies the following corollary.

**Corollary 2** (Unique parsimonious representation).: _Under the same assumption as in Theorem 1, there exists a constant \(c_{2}>0\) depending only on the sparsity of \(\mathbf{W}^{\star}\) such that for any \(\delta>0\), one has with probability at least \(1-\delta\),_

\[\widehat{\mathbf{a}}=\mathbf{\Pi}\mathbf{a}^{\star},\quad\widehat{\mathbf{W}}=\mathbf{W}^{\star} \mathbf{\Pi}^{\top},\quad\text{and}\quad\widehat{\mathbf{b}}=\mathbf{\Pi}\mathbf{b}^{\star}\]

_for some permutation matrix \(\mathbf{\Pi}\), provided that \(n>c_{2}\log^{4}\frac{p}{\delta}\)._

Corollary 2 states that among all the possible representations \(\mathbf{W}\) in (1) (with \(\xi_{i}=0\)), the one(s) with the smallest \(\ell_{1}\)-norm must be identical to \(\mathbf{W}^{\star}\) up to a column permutation with high probability. In other words, the most parsimonious representation (in the sense of \(\ell_{1}\) norm) of two-layer ReLU neural networks is unique. This observation addresses the questions raised in Section I.

It is worth mentioning that since the weight matrix \(\mathbf{W}\) of the neural network is row-sparse, Group-LASSO is a suitable alternative to LASSO. We leave the analysis of Group-LASSO for future study.

### _Elaboration on the main result_

Suppose that \(\mathbf{W}^{\star}\) has at most \(s\) nonzero entries. The following theorem is a more elaborated version of Theorem 1.

**Theorem 2**.: _There exist some constants \(c_{1},c_{2},c_{3}>0\) such that for any \(\delta>0\), one has with probability at least \(1-\delta\),_

\[\widehat{\mathbf{a}}=\mathbf{\Pi}\mathbf{a}^{\star}\text{ and }\|\widehat{\mathbf{W}}-\mathbf{W}^{ \star}\mathbf{\Pi}^{\top}\|_{\mathrm{F}}+\|\widehat{\mathbf{b}}-\mathbf{\Pi}\mathbf{b}^{\star} \|_{2}\leq c_{1}\sigma \tag{8}\]
_for some permutation \(\mathbf{\Pi}\in\{0,1\}^{r\times r}\), provided that Assumption 1 holds and_

\[n>c_{2}s^{3}r^{13}\log^{4}\frac{p}{\delta}\qquad\text{and}\qquad \sigma<\frac{c_{3}}{r}. \tag{9}\]

**Remark 3** (Sketch proof of Theorem 1).: _The proof of Theorem 1 is nontrivial and is included in the Appendix. Next, we briefly explain the sketch of the proof. First, we will define what we refer to as \(D_{1}\)-distance and \(D_{2}\)-distance between \((\mathbf{W},\mathbf{a},\mathbf{b})\) and \((\mathbf{W}^{\star},\mathbf{a}^{\star},\mathbf{b}^{\star})\). These distances can be regarded as the counterpart of the classical \(\ell_{1}\) and \(\ell_{2}\) distances between two vectors, but allowing the invariance under any permutation of neurons (Remark 2). Then, we let_

\[\Delta_{n}\coloneqq\frac{1}{n}\sum_{i=1}^{n} \bigg{[}\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x} _{i}+b_{j})\] \[-\sum_{j=1}^{r}a_{j}^{\star\top}\mathrm{relu}(\mathbf{w}_{j}^{\star \top}\mathbf{x}_{i}+b_{j}^{\star})\bigg{]}^{2},\]

_where \((\mathbf{W},\mathbf{a},\mathbf{b})\) is the solution of the problem in (3), and develop the following upper and lower bounds of it:_

\[\Delta_{n} \leq c_{6}\left(\frac{r}{S}+\frac{r\log^{3}\frac{p}{n\delta}}{n} \right)D_{1}^{2}+c_{6}\sigma^{2}\quad\text{and}\] \[\Delta_{n} \geq c_{4}\min\left\{\frac{1}{r},D_{2}^{2}\right\} \tag{10}\]

_hold with probability at least \(1-\delta\), provided that \(n\geq c_{5}S^{3}r^{4}\log^{4}\frac{p}{\delta}\), for some constants \(c_{4},c_{5},c_{6}\), and \(S\) to be specified. Here, the upper bound will be derived from a series of elementary inequalities. The lower bound is reminiscent of the Restricted Isometry Property (RIP) [15] for linear models. We will derive it from the lower bound of the population counterpart by concentration arguments, namely_

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top }\mathbf{x}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top }\mathbf{x}+b_{j}^{\star})\right]^{2}\] \[\geq c\min\left\{\frac{1}{r},D_{2}^{2}\right\},\]

_for some constant \(c>0\). The bounds in (10) imply that with high probability,_

\[c_{4}\min\left\{\frac{1}{r},D_{2}^{2}\right\}\leq c_{6}\left( \frac{r}{S}+\frac{r\log^{3}\frac{p}{n\delta}}{n}\right)D_{1}^{2}+c_{6}\sigma ^{2},\]

_Using this and an inequality connecting \(D_{1}\) and \(D_{2}\), we can prove the final result._

**Remark 4** (Alternative assumption and result).: _We provide an alternative to Theorem 2. Consider the following Assumption 1' as an alternative to Assumption 1._

**Assumption 1'**. For some constant \(B>0\),_

\[\|\mathbf{w}_{j}^{\star}\|_{2}\leq B\quad\text{and}\quad|b_{j}^{\star }|\leq B\quad\text{ for all }1\leq j\leq r.\]

_In addition,_

\[\mathbb{E}\Big{[}\langle\mathbf{a},\mathrm{relu}(\mathbf{W}^{\top}\mathbf{ x}+\mathbf{b})\rangle-\langle\mathbf{a}^{\star},\mathrm{relu}(\mathbf{W}^{\star\top}\mathbf{x}+\mathbf{b}^{ \star})\rangle\Big{]}^{2}\] \[\geq\psi D_{2}\left[(\mathbf{W},\mathbf{a},\mathbf{b}),(\mathbf{W}^{\star},\mathbf{ a}^{\star},\mathbf{b}^{\star})\right]^{2}, \tag{11}\]

_and_

\[n>\frac{c_{2}}{\psi}s^{3}r^{3}\log^{4}\frac{p}{\delta}. \tag{12}\]

_With Assumption 1' instead of Assumption 1, one can still derive the same result as Theorem 2. The proof of the above result is similar to that of Theorem 2, except that we insert Inequality (11) instead of Inequality (22) into (21) in Appendix A._

### _Variable selection_

To solve the optimization problem (3) in practice, we consider the following alternative problem,

\[\min_{\mathbf{W}\in\mathbb{R}^{p\times r},\mathbf{a}\in\mathbb{R}^{r}, \mathbf{b}\in\mathbb{R}^{r}} \bigg{\{}\frac{1}{n}\sum_{i=1}^{n}\biggl{(}y_{i}-\sum_{j=1}^{r}a_{j} \cdot\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})\biggr{)}^{2}\] \[+\lambda\|\mathbf{W}\|_{1}\bigg{\}}. \tag{13}\]

It has been empirically shown that algorithms such as the stochastic gradient descent can find a good approximate solution to the above optimization problem [34, 35]. Next, we will discuss some details regarding the variable selection using LASSO-regularized neural networks.

**Tuning parameters.** Given a labeled dataset in practice, we will need to tune hyper-parameters including the penalty term \(\lambda\), the number of neurons \(r\), learning rate, and the number of epochs. We suggest the usual approach that splits the available data into training and validation parts. The training data are used to estimate neural networks for a set of candidate hyper-parameters. The most suitable candidate will be identified based on the predictive performance on the validation data.

**Variable importance.** Inspired by Corollary 1, we interpret the \(\ell_{2}\)-norm of \(\widehat{\mathbf{w}}_{i,\cdot}\) as the importance of the \(i\)-th variable, for \(i=1,\ldots,p\). Corollary 1 indicates that we can accurately identify all the significant variables in \(\mathcal{S}^{\star}\) with high probability if we correctly set the cutoff value \(c_{1}\sigma\).

**Setting the cutoff value.** It is conceivable that variables with large importance are preferred over those with near-zero importance. This inspires us to cluster the variables into two groups based on their importance. Here, we suggest two possible approaches. The first is to use a data-driven approach such as \(k\)-means and Gaussian mixture model (GMM). The second is to manually set a threshold value according to domain knowledge on the number of important variables.

**Extension to deep neural networks.** Inspired by (13), we can intuitively generalize the proposed method to deep neural networks by penalizing the \(\ell_{1}\)-norm of the weight matrix in the input layer. Though we do not have a theoretical analysis for this broader setting, numerical studies show that it is still effective.

## III Experiments

We perform experimental studies to show the promising performance of the proposed variable selection method. We compare the variable selection accuracy and prediction performance of the proposed algorithm ('NN') with several
baseline methods, including LASSO ('LASSO'), orthogonal matching pursuit ('OMP'), random forest ('RF'), gradient boosting ('GB'), neural networks with group LASSO ('GLASSO') [5], group sparse regularization ('GSR') [8], and LNET ('LNET') [9]. The 'NN' hyperparameters to search over are the penalty term \(\lambda\in\{0.1,0.05,0.01,0.005,0.001\}\), the number of neurons \(r\in\{20,50,100\}\), the learning rate in \(\{0.01,0.005,0.001\}\), and the number of epochs in \(\{100,200,500\}\). Moreover, we extend 'NN' to a neural network that contains an additional hidden layer of ten neurons. We distinguish the proposed method with two-layer and three-layer neural networks by 'NN-2' and 'NN-3', respectively. Further experimental details are included in Appendix E.

### _Synthetic datasets_

#### Iv-A1 NN-generated dataset

The first experiment uses the data generated from Equation (1) with \(p=100\) variables, \(r=16\) neurons. The first \(10\) rows of neural coefficients \(\mathbf{W}\) are independently generated from the standard uniform distribution, and the remaining rows are zeros, representing \(10\) significant variables. The neural biases \(\mathbf{b}\) are also generated from the standard uniform distribution. The signs of neurons, \(\mathbf{a}\), follow an independent Bernoulli distribution. The training size is \(n=500\), and the test size is \(2000\). The noise is zero-mean Gaussian with standard deviation \(\sigma\) set to be \(0\), \(0.5\), \(1\), and \(5\). For each \(\sigma\), we evaluate its mean squared error on the test dataset and three quantities for variable selection: the number of correctly selected variables ('TP', the larger the better), wrongly selected variables ('FP', the smaller the better), and area-under-curve score ('AUC', the larger the better). Here, 'AUC' is evaluated based on the variable importance given by each method, which is detailed in Appendix E. The procedure is independently replicated \(20\) times.

The results are reported in Table I and Table II, which suggest that 'NN' has the best overall performance for both selection and prediction. In particular, 'NN-2' and 'NN-3' have almost the same performance among all situations, which empirically demonstrates that the proposed method also works for deeper neural networks. It is interesting to compare 'NN' with 'LNET': 'NN' has slightly higher test error than 'LNET' when the noise level is small, but a much smaller false positive rate and higher AUC score than 'LNET'. It indicates that 'NN' is more accurate for variable selection, while 'LNET' tends to over-select variables for better prediction accuracy. Also, all the methods perform worse as the noise level \(\sigma\) increases.
the quadratic factor \(x_{3}\). As for the prediction performance, neural network-based methods outperform other methods. In particular, 'NN' is better than 'GLASSO' and 'GSR', while 'LNET' exhibits better prediction and worse selection performance as seen in previous experiments.

### _BGSBoy dataset_

The BGSBoy dataset involves 66 boys from the Berkeley guidance study (BGS) of children born in 1928-29 in Berkeley, CA [36]. The dataset includes the height ('HT'), weight ('WT'), leg circumference ('LG'), strength ('ST') at different ages (2, 9, 18 years), and body mass index ('BMI18'). We choose 'BMI18' as the response, which is defined as follows.

\[\text{BMI18}=\text{WT18}/(\text{HT18}/100)^{2}, \tag{14}\]

where WT18 and HT18 denote the weight and height at the age of 18, respectively. In other words, 'WT18' and 'HT18' are sufficient for modeling the response among \(p=10\) variables. Other variables are correlated but redundant. The training size is \(n=44\) and the test size is 22. Other settings are the same as before. We compare the prediction performance and explore the three features which are most frequently selected by each method. The results are summarized in Table VII.

From the results, both linear and NN-based methods can identify 'WT18' and 'HT18' in all the replications. Meanwhile, tree-based methods may miss 'HT18' but select 'LG18' instead, which is only correlated with the response. Interestingly, we find that the linear methods still predict well in this experiment. A possible reason is that Equation (14) can be well-approximated by a first-order Taylor expansion on 'HT18' at the value of around \(180\), and the range of 'HT18' is within a small interval around \(180\).

### _UJIIndoorLoc dataset_

The UJIINdoorLoc dataset aims to solve the indoor localization problem via WiFi fingerprinting and other variables such as the building and floor numbers. A detailed description
can be found in [37]. Specifically, we have 520 Wireless Access Points (WAPs) signals (which are continuous variables) and 'FLOOR', 'BUILDING', 'SPACEID', 'RELATIVEPOSITION', 'USERID', and 'PHONEID' as categorical variables. The response variable is a user's longitude ('Longitude'). The dataset has \(19937\) observations. We randomly sample \(3000\) observations and split them into \(n=2000\) for training and \(1000\) for test. As part of the pre-processing, we create binary dummy variables for the categorical variables, which results in \(p=681\) variables in total. We explore the ten features that are most frequently selected by each method. We set the cutoff value as the tenth-largest variable importance. The procedure is independently replicated \(100\) times. The results are reported in Table VIII.

Based on the results, the 'NN' achieves the best prediction performance and significantly outperforms other methods. As for variable selection, since 'BUILDING' greatly influences the location from our domain knowledge, it is non-surprisingly selected by all methods in every replication. However, except for 'BUILDING', different methods select different variables with some overlapping, e.g., 'PHONEID_14' selected by 'GLASSO' and 'GB', 'USERID_16' selected by 'NN' and 'LASSO', which indicate the potentially important variables. 'LNET' again selects more variables than other methods. There are nearly 60 variables selected by 'LNET' in every replication. Nevertheless, those methods do not achieve an agreement for variable selection. 'NN' implies that all the WAPs signals are weak while categorical variables provide more information about the user location. Given the very high missing rate of WAPs signals (\(97\%\) on average, as reported in [37]), the interpretation of 'NN' seems reasonable.

### _Summary_

The experiment results show the following points. First, 'NN' can stably identify the important variables and have competitive prediction performance compared with the baselines. Second, the increase of the noise level will hinder both the selection and prediction performance. Third, the LASSO regularization is crucial for 'NN' to avoid over-fitting, especially for small data. Using group LASSO or a mixed type of penalty has a similar performance as 'NN', while 'LNET' tends to over-select importance variables. Fourth, the selection and prediction performances are often positively associated for 'NN', but may not be the case for baseline methods.

## IV Concluding Remarks

We established a theory for the use of LASSO in two-layer ReLU neural networks. In particular, we showed that the LASSO estimator could stably reconstruct the neural network coefficients and identify the critical underlying variables under reasonable conditions. We also proposed a practical method to solve the optimization and perform variable selection. We briefly remark on some interesting further work. First, a limitation of the work is that we considered only a small \(r\). An interesting future problem is to study \(r\) that may grow fast with \(p\) and \(n\). Second, our experiments show that the algorithm can be extended to deeper neural networks. It will be exciting to generalize the main theorem to the multi-layer cases.

The Appendix includes proofs and experimental details.

## Acknowledgment

The authors would like to thank the anonymous reviewers and Associate Editor for their constructive review comments. The work was supported in part by the U.S. National Science Foundation under Grant Numbers DMS-2134148 and CNS-2220286.

## Appendix A Analysis: proof of Theorem 2

Let \(\mathcal{S}\) be the index set with cardinality \(S\) consisting of the support for \(\mathbf{W}^{\star}\) and top entries of \(\widehat{\mathbf{W}}\), where \(S\) will be specified momentarily. Define

\[\mathbf{W}\coloneqq\widehat{\mathbf{W}}_{\mathcal{S}}\in\mathbb{R}^{p\times r},\]

and \(a_{j}=\widehat{a}_{j}\), \(b_{j}=\widehat{b}_{j}\). Define

\[d_{1}(\mathbf{w}_{1},a_{1},b_{1},\mathbf{w}_{2},a_{2},b_{2})\] \[=\left\{\begin{array}{ll}\|\mathbf{w}_{1}-\mathbf{w}_{2}\|_{1}+|b_{1}- b_{2}|&\text{if }a_{1}=a_{2};\\ \|\mathbf{w}_{1}\|_{1}+\|\mathbf{w}_{2}\|_{1}+|b_{1}|+|b_{2}|&\text{if }a_{1}\neq a_{2}, \end{array}\right. \tag{15}\]

and

\[d_{2}(\mathbf{w}_{1},a_{1},b_{1},\mathbf{w}_{2},a_{2},b_{2})\] \[=\left\{\begin{array}{ll}\sqrt{\|\mathbf{w}_{1}-\mathbf{w}_{2}\|_{2}^{2 }+|b_{1}-b_{2}|^{2}}&\text{if }a_{1}=a_{2};\\ 1&\text{if }a_{1}\neq a_{2}.\end{array}\right. \tag{16}\]

In addition, for permutation \(\pi\) on \([r]\), let

\[D_{1} \coloneqq\min_{\pi}\sum_{j=1}^{r}d_{1}(\mathbf{w}_{\pi(j)},a_{\pi(j)}, b_{\pi(j)},\mathbf{w}_{j}^{\star},a_{j}^{\star},b_{j}^{\star}), \tag{17a}\] \[D_{2} \coloneqq\min_{\pi}\sqrt{\sum_{j=1}^{r}d_{2}(\mathbf{w}_{\pi(j)},a_{ \pi(j)},b_{\pi(j)},\mathbf{w}_{j}^{\star},a_{j}^{\star},b_{j}^{\star})^{2}} \tag{17b}\]
denote the \(D_{1}\)-distance and \(D_{2}\)-distance between \((\mathbf{W},\mathbf{a},\mathbf{b})\) and \((\mathbf{W}^{\star},\mathbf{a}^{\star},\mathbf{b}^{\star})\), respectively. Then, one has the following bounds.

**Lemma 1**.: _For any \(\mathbf{W}\in\mathbb{R}^{p\times r}\) with \(\left\|\mathbf{W}\right\|_{0}\leq S\), there exists some universal constants \(c_{4},c_{5}>0\) such that_

\[\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}( \mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\right]^{2}\] \[\geq c_{4}\min\left\{\frac{1}{r},D_{2}^{2}\right\} \tag{18}\]

_holds with probability at least \(1-\delta\) provided that_

\[n\geq c_{5}S^{3}r^{4}\log^{4}\frac{p}{\delta}. \tag{19}\]

**Lemma 2**.: _Then, there exists a universal constant \(c_{6}>0\) such that_

\[\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}( \mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\right]^{2}\] \[\leq c_{6}\left(\frac{r}{S}+\frac{r\log^{3}\frac{p}{n\delta}}{n }\right)D_{1}^{2}+c_{6}\sigma^{2}\]

_holds with probability at least \(1-\delta\)._

By comparing the bounds given in Lemma 1 and 2, one has

\[c_{4}\min\left\{\frac{1}{r},D_{2}^{2}\right\}\leq c_{6}\left(\frac{r}{S}+\frac {r\log^{3}\frac{p}{n\delta}}{n}\right)D_{1}^{2}+c_{6}\sigma^{2},\]

provided that

\[n>c_{5}S^{3}r^{4}\log^{4}\frac{p}{\delta}.\]

Let \(\widehat{\mathcal{S}}^{\star}\) be the index set with cardinality \(2s\) consisting of the support for \(\mathbf{W}^{\star}\) and top entries of \(\widehat{\mathbf{W}}\). In addition, let \(D_{1}^{\star}\) and \(D_{2}^{\star}\) denote the \(D_{1}\)-distance and \(D_{2}\)-distance between \(\left(\widehat{\mathbf{W}}_{\widehat{\mathcal{S}}},\widehat{\mathbf{a}},\widehat{\mathbf{ b}}\right)\) and \(\left(\mathbf{W}^{\star},\mathbf{a}^{\star},\mathbf{b}^{\star}\right)\) in a similar way as (17). Observing the fact that for \(S\geq 2s\), one has \(\mathcal{S}^{\star}\subset\widehat{\mathcal{S}}^{\star}\subset\mathcal{S}\), we have

\[\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_{2}\geq\|\mathbf{w}_{j,\widehat{\mathcal{S}}^{ \star}}-\mathbf{w}_{j}^{\star}\|_{2}=\|\widehat{\mathbf{w}}_{j}-\mathbf{w}_{j}^{\star}\|_ {2},\]

after some permutation, and then

\[D_{2}^{\star}\leq D_{2}.\]

In addition, after some permutation, we have \(D_{1}^{\star}\geq\left\|\widehat{\mathbf{W}}_{\widehat{\mathcal{S}}^{\star}}- \mathbf{W}^{\star}\right\|_{1}\geq\left\|\mathbf{W}^{\star}\right\|_{1}-\left\| \widehat{\mathbf{W}}_{\mathcal{S}^{\star}}\right\|_{1}\) and \(\left\|\mathbf{W}\right\|_{1}\leq\left\|\widehat{\mathbf{W}}\right\|_{1}\leq\left\| \mathbf{W}^{\star}\right\|_{1}\). Then,

Combined with Lemma 3 in Appendix D, the above results give

\[D_{2}^{\star}\leq\frac{2c_{6}}{c_{4}}\sigma,\]

provided that for some constant \(c_{7}>0\)

\[n\geq c_{5}S^{3}\log^{4}\frac{p}{\delta}\qquad\text{with }S\geq c_{7}sr,\]

such that

\[c_{6}\left(\frac{r}{S}+\frac{r\log^{3}\frac{p}{n\delta}}{n}\right)D_{1}^{\star 2 }\leq\frac{c_{4}}{8}D_{2}^{\star 2}.\]
Then, we can conclude the proof since after appropriate permutation,

\[\|\widehat{\mathbf{W}}-\mathbf{W}^{\star}\|_{\mathrm{F}}\leq 2\|\widehat{\mathbf{W}}_{ \widehat{\mathcal{S}}^{\star}}-\mathbf{W}^{\star}\|_{\mathrm{F}}.\]

## Appendix B Proof of Lemma 1 (Lower bound)

This can be seen from the following three properties.

* Consider the case that \[D_{1}\leq\epsilon=\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta}}}.\] With probability at least \(1-\delta\), \[\frac{1}{n}\sum_{i=1}^{n} \bigg{[}\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}_ {i}+b_{j})\] \[-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x}_{i}+b_{j}^{\star})\bigg{]}^{2}\] \[=\frac{D_{1}^{2}}{\epsilon^{2}}\frac{1}{n}\sum_{i=1}^{n} \bigg{[}\sum_{j=1}^{r}a_{j}\mathrm{relu}(\widetilde{\mathbf{w}}_{j}^{\top}\mathbf{x }_{i}+\widetilde{b}_{j})\] \[-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x}_{i}+b_{j}^{\star})\bigg{]}^{2},\] (20) where \(\widetilde{\mathbf{w}}_{j}=\mathbf{w}_{j}^{\star}+\frac{\epsilon}{D_{1}}\left(\mathbf{w} _{j}-\mathbf{w}_{j}^{\star}\right)\) and \(\widetilde{b}_{j}=b_{j}^{\star}+\frac{\epsilon}{D_{1}}\left(b_{j}-b_{j}^{ \star}\right)\).
* For any \(\epsilon>0\) and \[D_{1}\geq\frac{\epsilon}{\sqrt{\frac{S}{n}\log\frac{pr}{S}\log\frac{BS}{ \epsilon\delta}}},\] there exists some universal constant \(C_{1}>0\), such that with probability at least \(1-\delta\), \[\frac{1}{n}\sum_{i=1}^{n} \bigg{[}\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}_ {i}+b_{j})\] \[-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x}_{i}+b_{j}^{\star})\bigg{]}^{2}\] \[\geq\mathbb{E} \left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_ {j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}+b_ {j}^{\star})\right]^{2}\] \[-C_{1}D_{1}^{2}\log\frac{pn}{\delta}\sqrt{\frac{S}{n}\log\frac{pr} {S}\log\frac{BS}{\epsilon\delta}}.\] (21)
* For some universal constant \(C_{2}>0\) \[\mathbb{E} \left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_ {j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}+b_ {j}^{\star})\right]^{2}\] \[\geq C_{2}\min\left\{\frac{1}{r},D_{2}^{2}\right\}.\] (22)

_Putting the above together._ Let

\[\epsilon=C_{3}\frac{\delta}{nr}\sqrt{\frac{S}{n}\log\frac{BnS}{\delta}},\]

for some universal constant \(C_{3}>0\) such that

\[\frac{\epsilon}{\sqrt{\frac{S}{n}\log\frac{pr}{S}\log\frac{BS}{\epsilon\delta} }}<\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta}}}.\]

Inserting (22) into (21) gives that

\[\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}( \mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\right]^{2}\] \[\geq C_{2}\min\left\{\frac{1}{r},D_{2}^{2}\right\}-C_{1}D_{1}^{2} \log\frac{pn}{\delta}\sqrt{\frac{S}{n}\log\frac{pr}{S}\log\frac{BS}{\epsilon \delta}}\] \[\geq\frac{C_{2}}{2}\min\left\{\frac{1}{r},D_{2}^{2}\right\}, \tag{23}\]

holds with probability at least \(1-\delta\) provided that for some constant \(C_{4}>0\)

\[n\geq C_{4}S^{3}r^{4}\log\frac{pr}{S}\log\frac{BS}{\epsilon \delta}\log^{2}\frac{pn}{\delta}\quad\text{and}\] \[D_{1}\geq\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta }}}.\]

Here, the last line holds due to Lemma 3 and we assume that \(\max\left\{\|\mathbf{W}\|_{\infty},\|\mathbf{b}\|_{\infty}\right\}\) is bounded by some constant. On the other hand, if

\[D_{1}<\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta}}},\]

it follows from (20) and (23) that

\[\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}( \mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\right]^{2}\] \[\geq\frac{D_{1}^{2}}{\epsilon^{2}}\frac{C_{2}}{2}\min\left\{ \frac{1}{r},\widetilde{D}_{2}^{2}\right\}=\frac{C_{2}}{2}D_{2}^{2}.\]

Summing up, we conclude the proof by verifying (20), (21), and (22) below.

### _Proof of (20)_

Since \(D_{1}\leq\epsilon=\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta}}}\), without loss of generality, we assume that \(a_{j}=a_{j}^{\star}\) for \(1\leq j\leq r\), and

\[D_{1}=\sum_{j=1}^{r}\left(\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_{1}+|b_{j}-b_{j}^{ \star}|\right)\leq\epsilon.\]

By taking union bound, with probability at least \(1-\frac{\delta}{2}\), one has for all \(1\leq i\leq n\) and \(1\leq j\leq r\),

\[\left|\mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star}\right|>\frac{\delta}{2nr} \sqrt{\frac{\pi}{2}},\]

since \(\|\mathbf{w}_{j}^{\star}\|_{2}\geq 1\) and \(\mathbf{x}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{I})\). In addition, for all \(1\leq i\leq n\) and \(1\leq j\leq r\),

\[\left|\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j}-\mathbf{w}_{j}^{\star\top}\mathbf{x}_ {i}-b_{j}^{\star}\right| \leq\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_{1}\|\mathbf{x}_{i}\|_{\infty}+|b_{j}-b_{j }^{\star}|\] \[\leq\epsilon\sqrt{2\log\frac{4pn}{\delta}}\]
holds with probability at least \(1-\frac{\delta}{2}\). Here, the last inequality comes from the fact that with probability at least \(1-\frac{\delta}{2}\),

\[\|\mathbf{x}_{i}\|_{\infty}\leq\sqrt{2\log\frac{4pn}{\delta}}\qquad\text{for all }1\leq i\leq n. \tag{24}\]

Putting together, we have with probability at least \(1-\delta\),

\[u(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})=u(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j}^{ \star}), \tag{25}\]

provided that

\[\epsilon\leq\frac{\delta}{4nr}\sqrt{\frac{\pi}{\log\frac{4pn}{\delta}}}.\]

Note that \(u(x)=1\) if \(x>0\), and \(u(x)\!=0\) if \(x\leq 0\). Then combining with the definition of \(\widetilde{\mathbf{w}}_{j}\) and \(\widetilde{b}_{j}\), the above property yields

\[\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}( \mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j}^{\star})\right]^{2}\] \[=\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}^{\star}u(\mathbf{ w}_{j}^{\top}\mathbf{x}_{i}+b_{j}^{\star})(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j}- \mathbf{w}_{j}^{\top}\mathbf{x}_{i}-b_{j}^{\star})\right]^{2}\] \[=\frac{D_{1}^{2}}{\epsilon^{2}}\frac{1}{n}\sum_{i=1}^{n}\biggl{[} \sum_{j=1}^{r}a_{j}^{\star}u(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j}^{\star})\] \[\qquad\qquad\qquad\times(\widetilde{\mathbf{w}}_{j}^{\top}\mathbf{x}_{i}+ \widetilde{b}_{j}-\mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}-b_{j}^{\star})\biggr{]}^{2}\] \[=\frac{D_{1}^{2}}{\epsilon^{2}}\frac{1}{n}\sum_{i=1}^{n}\biggl{[} \sum_{j=1}^{r}a_{j}\mathrm{relu}(\widetilde{\mathbf{w}}_{j}^{\top}\mathbf{x}_{i}+ \widetilde{b}_{j})\] \[\qquad\qquad\qquad-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{ w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\biggr{]}^{2},\]

and the claim is proved. Here, the last equality holds due to (25) and \(a_{j}=a_{j}^{\star}\) for \(j=1,\ldots,r\).

### _Proof of (21)_

Notice that

\[\left|a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_{j})-a_{j}^{ \star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star})\right|\] \[\leq\left\{\begin{array}{ll}\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_ {1}\|\mathbf{x}\|_{\infty}+|b_{j}-b_{j}^{\star}|&\text{if }a_{j}=a_{j}^{\star},\\ \left(\|\mathbf{w}_{j}\|_{1}+\|\mathbf{w}_{j}^{\star}\|_{1}\right)\|\mathbf{x}\|_{\infty}+ |b_{j}|+|b_{j}^{\star}|&\text{if }a_{j}\neq a_{j}^{\star},\end{array}\right.\]

which leads to

\[\left|\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b _{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}+b _{j}^{\star})\right|\] \[\leq D_{1}\max\left\{\|\mathbf{x}\|_{\infty},1\right\}. \tag{26}\]

For any fixed \((\mathbf{W},\mathbf{a},\mathbf{b})\), let

\[z_{i}:=\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})- \sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j }^{\star}),\]

and define the following event set

\[\mathcal{E}:=\left\{\|\mathbf{x}_{i}\|_{\infty}\leq\sqrt{2\log\frac{4pn}{\delta}} \qquad\text{for all }1\leq i\leq n\right\}.\]

Then, with probability at least \(1-\delta\),

\[\frac{1}{n}\sum_{i=1}^{n}\left(z_{i}^{2}-\mathbb{E}\left[z_{i}^{ 2}\right]\right) \tag{27}\] \[=\frac{1}{n}\sum_{i=1}^{n}\biggl{\{}z_{i}^{2}\mathbb{1}(\mathcal{ E})-\mathbb{E}\left[z_{i}^{2}\mathbb{1}(\mathcal{E})\right]-\mathbb{E}\left[z_{i}^{2} \mathbb{1}(\overline{\mathcal{E}})\right]\biggr{\}}\] \[\geq-4D_{1}^{2}\log\frac{4pn}{\delta}\sqrt{\frac{1}{n}\log\frac{ 2}{\delta}}-D_{1}^{2}\frac{\delta}{n}\] \[\geq-5D_{1}^{2}\log\frac{4pn}{\delta}\sqrt{\frac{1}{n}\log\frac{ 2}{\delta}}. \tag{28}\]

Here, the first line holds due to (24); the last line comes from Hoeffding's inequality, and the fact that

\[\left|\mathbb{E}\left[z_{i}^{2}\mathbb{1}(\overline{\mathcal{E}})\right]\right| \leq D_{1}^{2}\left|\mathbb{E}\left[\|\mathbf{x}_{i}\|_{\infty}^{2}1( \|\mathbf{x}_{i}\|_{\infty}>\sqrt{2\log\frac{4pn}{\delta}})\right]\right|\] \[\leq D_{1}^{2}\int_{\sqrt{2\log\frac{4pn}{\delta}}}^{\infty}x^{2}d \mathbb{P}(\|\mathbf{x}_{i}\|_{\infty}<x)\] \[\leq D_{1}^{2}\int_{\sqrt{2\log\frac{4pn}{\delta}}}^{\infty}4xp \exp(-\frac{x^{2}}{2})dx\leq D_{1}^{2}\frac{\delta}{n}.\]

In addition, consider the following \(\epsilon\)-net

\[\mathcal{N}_{\epsilon}\coloneqq\bigg{\{} (\mathbf{W},\mathbf{a},\mathbf{b}):|W_{ij}|\in\frac{\epsilon}{r+S}\left[\left. \frac{B(r+S)}{\epsilon}\right]\right\rfloor,\] \[\|\mathbf{W}\|_{0}\leq S,|b_{j}|\in\frac{\epsilon}{r+S}\left[\left. \frac{B(r+S)}{\epsilon}\right],|a_{j}|=1\right\}\!,\]

where \([n]:=\{1,2,\ldots,n-1\}\). Then, for all \((\mathbf{W},\mathbf{a},\mathbf{b})\) with \(\|\mathbf{W}\|_{1}\leq B\) and \(\|\mathbf{b}\|_{1}\leq B\), there exists some point, denoted by \(\left(\widetilde{\mathbf{W}},\widetilde{\mathbf{a}},\widetilde{\mathbf{b}}\right)\), in \(\mathcal{N}_{\epsilon}\) whose \(D_{1}\)-distance from \((\mathbf{W},\mathbf{a},\mathbf{b})\) is less than \(\epsilon\). For simplicity, define

\[z_{i}\coloneqq\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top} \mathbf{x}_{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{ \star\top}\mathbf{x}_{i}+b_{j}^{\star}),\] \[\widetilde{z}_{i}\coloneqq\sum_{j=1}^{r}\widetilde{a}_{j}\mathrm{relu }(\widetilde{\mathbf{w}}_{j}^{\top}\mathbf{x}_{i}+\widetilde{b}_{j})-\sum_{j=1}^{r}a_{j}^ {\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}_{i}+b_{j}^{\star}).\]

Similar to (26), we can derive that

\[\left|\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_{j})-\sum_{j =1}^{r}\widetilde{a}_{j}\mathrm{relu}(\widetilde{\mathbf{w}}_{j}^{\top}\mathbf{x}+ \widetilde{b}_{j})\right|\] \[\leq\epsilon\max\left\{\|\mathbf{x}\|_{\infty},1\right\},\]

which implies

\[\left|z_{i}^{2}-\widehat{z}_{i}^{2}\right|\leq\epsilon\left(\epsilon+D_{1} \right)\max\left\{\|\mathbf{x}_{i}\|_{\infty}^{2},1\right\},\]

and then with probability at least
for some universal constant \(C_{5}>0\). Combining (28), (29), and (30) leads to

\[\frac{1}{n}\sum_{i=1}^{n}\left(z_{i}^{2}-\mathbb{E}\left[z_{i}^{2} \right]\right)\geq -5\left(\epsilon+D_{1}\right)^{2}\log\frac{4pn}{\delta}\sqrt{ \frac{1}{n}\log\frac{2\left|\mathcal{N}_{c}\right|}{\delta}}\] \[-4\epsilon\left(\epsilon+D_{1}\right)\log\frac{4pn}{\delta}.\]

It follows that (21) holds.

### _Proof of (22)_

We first consider a simple case that \(b_{j}=0\) and \(b_{j}^{\star}=0\) for \(1\leq j\leq r\), and show that for some small constant \(C_{6}>0\),

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{ \top}\mathbf{x})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x})\right]^{2}\] \[\geq C_{6}\min\left\{\frac{1}{r},D_{2}^{2}\right\}. \tag{31}\]

Next, we will assume that

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{ \top}\mathbf{x})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x})\right]^{2}\leq\frac{C_{6}}{r}.\]

Otherwise, Inequality (31) already holds. According to Lemma 4, one has for any constant \(k\geq 0\), there exists some constant \(\alpha_{k}>0\) such that

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{ \top}\mathbf{x})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x})\right]^{2}\] \[\geq\alpha_{k}\bigg{\|}\sum_{j=1}^{r}a_{j}\|\mathbf{w}_{j}\|_{2}\big{(} \frac{\mathbf{w}_{j}}{\|\mathbf{w}_{j}\|_{2}}\big{)}^{\otimes 2k}\] \[\qquad-\sum_{j=1}^{r}a_{j}^{\star}\|\mathbf{w}_{j}^{\star}\|_{2}\big{(} \frac{\mathbf{w}_{j}^{\star}}{\|\mathbf{w}_{j}^{\star}\|_{2}}\big{)}^{\otimes 2k} \bigg{\|}_{F}^{2}. \tag{32}\]

Assumption 1 tells us that for any integer \(k\geq\frac{2}{\omega}\),

\[\big{|}\langle\mathbf{v}_{j_{1}}^{\star},\mathbf{v}_{j_{2}}^{\star}\rangle\big{|}\leq \frac{1}{r^{2}}. \tag{33}\]

where

\[\mathbf{v}_{j}\coloneqq\mathrm{vec}\left(\big{(}\frac{\mathbf{w}_{j}}{\|\mathbf{w}_{j}\|_ {2}}\big{)}^{\otimes k}\right)\quad\text{with}\quad\beta_{j}\coloneqq a_{j}\| \mathbf{w}_{j}\|_{2},\]

and

\[\mathbf{v}_{j}^{\star}\coloneqq\mathrm{vec}\left(\big{(}\frac{\mathbf{w}_{j}^{\star}} {\|\mathbf{w}_{j}^{\star}\|_{2}}\big{)}^{\otimes k}\right)\quad\text{with}\quad \beta_{j}^{\star}\coloneqq a_{j}^{\star}\|\mathbf{w}_{j}^{\star}\|_{2}.\]

Then, (32) gives

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{ \top}\mathbf{x})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top} \mathbf{x})\right]^{2}\] \[\geq\alpha_{3k}\left\|\sum_{j=1}^{r}\beta_{j}\mathbf{v}_{j}^{\otimes 6 }-\sum_{j=1}^{r}\beta_{j}^{\star}\mathbf{v}_{j}^{\star\otimes 6}\right\|_{F}^{2}.\]

Define

\[\mathbb{S}_{+}\coloneqq\mathrm{span}\left\{\mathbf{v}_{j}\right\}_{j;\beta_{j}>0} \qquad\mathbb{S}_{-}\coloneqq\mathrm{span}\left\{\mathbf{v}_{j}\right\}_{j:\beta_ {j}<0},\]

and

\[\mathbb{S}_{+}^{\star}\coloneqq\mathrm{span}\left\{\mathbf{v}_{j}^{\star}\right\}_{j :\beta_{j}^{\star}>0}\qquad\mathbb{S}_{-}^{\star}\coloneqq\mathrm{span} \left\{\mathbf{v}_{j}^{\star}\right\}_{j:\beta_{j}^{\star}<0}.\]

Let \(\mathbf{P}_{\mathbb{S}}\) and \(\mathbf{P}_{\mathbb{S}}^{\perp}\) denote the projection onto and perpendicular to the subspace \(\mathbb{S}\), respectively. By noticing that \(\mathbf{P}_{\mathbb{S}}^{\perp}\mathbf{v}_{j}=\mathbf{0}\) for \(j\) obeying \(\beta_{j}<0\), and \(\mathbf{P}_{\mathbb{S}_{+}^{\star}}^{\perp}\mathbf{v}_{j}^{\star}=\mathbf{0}\) for \(j\) obeying \(\beta_{j}^{\star}>0\), one has

\[\left\|\sum_{j=1}^{r}\beta_{j}\mathbf{v}_{j}^{\otimes 6}-\sum_{j=1}^{r} \beta_{j}^{\star}\mathbf{v}_{j}^{\star\otimes 6}\right\|_{F}^{2}\] \[\geq\left\|\sum_{j:\beta_{j}>0}\beta_{j}\big{(}\mathbf{P}_{\mathbb{S} _{-}}^{\perp}\mathbf{v}_{j}\big{)}^{\otimes 2}\otimes\big{(}\mathbf{P}_{\mathbb{S}_{+}}^{ \perp}\mathbf{v}_{j}\big{)}^{\otimes 4}\right.\] \[\qquad-\sum_{j:\beta_{j}^{\star}<0}\beta_{j}^{\star}\big{(}\mathbf{P}_{ \mathbb{S}_{-}}^{\perp}\mathbf{v}_{j}^{\star}\big{)}^{\otimes 2}\otimes\big{(}\mathbf{P}_{\mathbb{S}_{+}}^{ \perp}\mathbf{v}_{j}^{\star}\big{)}^{\otimes 4}\bigg{\|}_{F}^{2}\] \[\geq\sum_{j:\beta_{j}^{\star}<0}\left\|\mathbf{P}_{\mathbb{S}_{-}}^{ \perp}\mathbf{v}_{j}^{\star}\right\|_{2}^{4},\]

where the penultimate inequality holds since the inner product between every pair of terms is positive, and the last inequality comes from the facts that \(|\beta_{j}^{\star}|\geq 1\) and (33).

Moreover, by means of AM-GM inequality and (33), one can see that

\[\sum_{j:\beta_{j}^{\star}<0}\left\|\mathbf{P}_{\mathbb{S}_{-}}^{\perp }\mathbf{v}_{j}^{\star}\right\|_{2}^{4} \geq\frac{1}{r}\Big{(}\sum_{j:\beta_{j}^{\star}<0}\left\|\mathbf{P}_{ \mathbb{S}_{-}}^{\perp}\mathbf{v}_{j}^{\star}\right\|_{2}^{2}\Big{)}^{2}\] \[=\frac{1}{r}\left\|\mathbf{P}_{\mathbb{S}_{-}}^{\perp}\big{[}\mathbf{v}_{j }^{\star}\big{]}_{j:\beta_{j}^{\star}<0}\right\|_{F}^{4}\] \[\geq\frac{1}{2r}\left\|\mathbf{P}_{\mathbb{S}_{-}}^{\perp}\mathbf{P}_{ \mathbb{S}_{-}}^{\perp}\right\|_{F}^{4}.\]

Then combining with (31), the above result and the counterpart for \(\beta_{j}^{\star}>0\) lead to

\[\dim(\mathbb{S}_{-})\geq\dim(\mathbb{S}_{-}^{\star})\quad\text{and}\quad\dim( \mathbb{S}_{+})\geq\dim(\mathbb{S}_{+}^{\star}),\]

which gives

\[\dim(\mathbb{S}_{-})=\dim(\mathbb{S}_{-}^{\star})\quad\text{and}\quad\dim( \mathbb{S}_{+})=\dim(\mathbb{S}_{+}^{\star}).\]

Furthermore, for some small constant \(C_{6}>0\), we have

\[\mathrm{dist}(\mathbb{S}_{-},\mathbb{S}_{-}^{\star})\leq C_{6}\quad\text{and} \quad\mathrm{dist}(\mathbb{S}_{+},\mathbb{S}_{+}^{\star})\leq C_{6}.\]

Let \(\mathbf{P}_{i}^{\perp}\) denote the projection perpendicular to

\[\mathrm{span}\left\{\mathbf{v}_{j}^{\star}\right\}_{j\neq i:\beta_{j}^{\star}>0},\]

and

\[\gamma_{j}\coloneqq\frac{\beta_{j}\langle\mathbf{P}_{\mathbb{S}_{-}}^{\perp}\mathbf{v}_{j },P_{\mathbb{S}_{-}}^{\perp}\mathbf{v}_{i}^{\star}\rangle^{2}\langle\mathbf{P}_{i}^{ \perp}\mathbf{v}_{i},P_{\mathbb{S}_{-}}^{\perp}\mathbf{v}_{i}^{\star}\rangle^{2}}{ \left\|\mathbf{P}_{\mathbb{S}
Then for any \(i\),

\[\left\|\sum_{j=1}^{r}\beta_{j}\mathbf{v}_{j}^{\otimes 6}-\sum_{j=1}^{r} \beta_{j}^{*}\mathbf{v}_{j}^{*\otimes 6}\right\|_{F}^{2}\] \[\geq\left\|\sum_{j;\beta_{j}>0}\beta_{j}\big{(}\mathbf{P}_{\mathbb{S} \_}}^{\perp}\mathbf{v}_{j}\big{)}^{\otimes 2}\otimes\mathbf{v}_{j}^{\otimes 4}\right.\] \[\left.\qquad-\sum_{j=1}^{r}\beta_{j}^{*}\big{(}\mathbf{P}_{\mathbb{S} \_}}^{\perp}\mathbf{v}_{j}^{*}\big{)}^{\otimes 2}\otimes\mathbf{v}_{j}^{*\otimes 4}\right\|_{F}^{2}\] \[\geq\frac{1}{2}\bigg{\|}\sum_{j;\beta_{j}>0}\beta_{j}\big{(}\mathbf{P} _{\mathbb{S}\_}}^{\perp}\mathbf{v}_{j}\big{)}^{\otimes 2}\otimes\mathbf{v}_{j}^{ \otimes 4}\] \[\geq\frac{1}{2}\bigg{\|}\sum_{j;\beta_{j}>0}\beta_{j}\big{(}\mathbf{P} _{\mathbb{S}\_}}^{\perp}\mathbf{v}_{j}\big{)}^{\otimes 2}\otimes\mathbf{(}\mathbf{P}_{i}^{ \perp}\mathbf{v}_{i}\big{)}^{\otimes 2}\otimes\mathbf{v}_{j}^{\otimes 2}\] \[\qquad\qquad-\beta_{i}^{*}\big{(}\mathbf{P}_{\mathbb{S}\_}}^{\perp} \mathbf{v}_{i}\big{)}^{\otimes 2}\otimes\mathbf{(}\mathbf{P}_{i}^{\perp}\mathbf{v}_{i}\big{)}^{ \otimes 2}\otimes\mathbf{v}_{i}^{*\otimes 2}\bigg{\|}_{F}^{2}\] \[\geq\frac{1}{2}\left\|\sum_{j;\beta_{j}>0}\gamma_{j}\mathbf{v}_{j}^{ \otimes 2}-\beta_{i}^{*}\big{\|}\mathbf{P}_{\mathbb{S}\_}}^{\perp}\mathbf{v}_{i}^{*} \big{\|}_{2}^{2}\big{\|}\mathbf{P}_{i}^{*}\mathbf{v}_{i}^{*}\right\|^{2}\mathbf{v}_{i}^{* \otimes 2}\Bigg{\|}_{F}^{2}\,,\]

which, together with (31), implies that there exists some \(j\) such that

\[\|\sqrt{\beta_{j}}\mathbf{v}_{j}-\sqrt{\beta_{i}^{*}}\mathbf{v}_{i}^{*}\|_{2}^{2}\leq \frac{1}{r}.\]

Without loss of generality, assume that

\[\|\sqrt{\beta_{j}}\mathbf{v}_{j}-\sqrt{\beta_{j}^{*}}\mathbf{v}_{j}^{*}\|_{2}^{2}\leq \frac{1}{r},\qquad\text{for all }1\leq j\leq r. \tag{34}\]

Then

\[\mathbb{E}\big{[}\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{ \top}\mathbf{x})-\sum_{j=1}^{r}a_{j}^{*}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}) \big{]}^{2}\] \[\geq\alpha_{k}\left\|\sum_{j=1}^{r}\beta_{j}\mathbf{v}_{j}\mathbf{v}_{j}^{ \top}-\sum_{j=1}^{r}\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{*\top}\right\|_{F}^ {2}\] \[\geq\alpha_{k}\sum_{j=1}^{r}\big{\|}\beta_{j}\mathbf{v}_{j}\mathbf{v}_{j}^ {\top}-\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{*\top}\big{\|}_{F}^{2}\] \[\quad-\frac{\alpha_{k}}{2r}\left(\sum_{j=1}^{r}\big{\|}\beta_{j} \mathbf{v}_{j}\mathbf{v}_{j}^{\top}-\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{*\top} \big{\|}_{F}^{2}\right)^{2}\] \[\geq\frac{\alpha_{k}}{2}\sum_{j=1}^{r}\big{\|}\beta_{j}\mathbf{v}_{j} \mathbf{v}_{j}^{\top}-\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{*\top}\big{\|}_{F}^{2}\,.\]

Here, the first line comes from (32); the second line holds through the following claim

\[\leq\frac{1}{2r}\|\beta_{j_{1}}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}-\beta_{ j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{*\top}\|_{2}\|\beta_{j_{2}}\mathbf{v}_{j_{2}}\mathbf{v}_{j_{2}} ^{\top}-\beta_{j_{2}}^{*}\mathbf{v}_{j_{2}}^{*}\mathbf{v}_{j_{2}}^{*\top}\|_{2}\]

since for \(\mathbf{\delta}_{j}\coloneqq\sqrt{\beta_{j}}\mathbf{v}_{j}-\sqrt{\beta_{j}^{*}}\mathbf{v}_{j }^{*}\),

\[\beta_{j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}-\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{j}^{* \top}=\mathbf{\delta}_{j}\mathbf{\delta}_{j}^{\top}+\sqrt{\beta_{j}^{*}}\mathbf{\delta}_{j} \mathbf{v}_{j}^{*\top}+\sqrt{\beta_{j}^{*}}\mathbf{\delta}_{j}\mathbf{\delta}_{j}^{*\top}.\]

Then the conclusion is obvious by noticing that

\[\big{\|}\beta_{j}\mathbf{v}_{j}\mathbf{v}_{j}^{\top}-\beta_{j}^{*}\mathbf{v}_{j}^{*}\mathbf{v}_{ j}^{*\top}\big{\|}_{F}\geq\|\mathbf{w}_{j}-\mathbf{w}_{j}^{*}\|_{2}.\]

Finally, we analyze the general case with \(b_{j},b_{j}^{*}\neq 0\), which is similar to the above argument. For simplicity, we only explain the different parts here. According to Lemma 4, one has for any constant \(k\geq 0\), there exists some constant \(\alpha_{k}>0\) and some function \(f_{k}:\mathbb{R}\to\mathbb{R}\) such that

\[\mathbb{E}\left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top }\mathbf{x}+b_{j})-\sum_{j=1}^{r}a_{j}^{*}\mathrm{relu}(\mathbf{w}_{j}^{*\top}\mathbf{x}+b_ {j}^{*})\right]^{2}\] \[\geq\sum_{k\geq\frac{1}{2}\omega}^{\infty}\left\|\sum_{j=1}^{r}a_{j }f_{k}(\frac{b_{j}}{\|\mathbf{w}_{j}\|_{2}})\|\mathbf{w}_{j}\|_{2}\big{(}\frac{\mathbf{w}_{j} }{\|\mathbf{w}_{j}\|_{2}})^{\otimes k}\right.\] \[\left.\qquad\qquad-\sum_{j=1}^{r}a_{j}^{*}f_{k}(\frac{b_{j}^{*}}{ \|\mathbf{w}_{j}^{*}\|_{2}})\|\mathbf{w}_{j}^{*}\|_{2}\big{(}\frac{\mathbf{w}_{j}^{*}}{\|\mathbf{w} _{j}^{*}\|_{2}})^{\otimes k}\right\|_{F}^{2}\] \[\gtrsim\sum_{j=1}^{r}\sum_{k\geq\frac{12}{r}}^{\infty}\left\|a_{j }f_{k}(\frac{b_{j}}{\|\mathbf{w}_{j}\|_{2}})\mathbf{w}_{j}-a_{j}^{*}f_{k}(\frac{b_{j}^{*}}{ \|\mathbf{w}_{j}^{*}\|_{2}})\mathbf{w}_{j}^{*}\right\|_{F}^{2}\] \[\gtrsim\sum_{j=1}^{r}\inf_{R_{l}(\mathbf{x})}\mathbb{E}\Big{[}a_{j} \mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_{j}) \tag{35}\] \[\left.\qquad\qquad-a_{j}^{*}\mathrm{relu}(\mathbf{w}_{j}^{*\top}\mathbf{x} +b_{j}^{*})-R_{l}(\mathbf{x})\right]^{2}\] \[\gtrsim\sum_{j=1}^{r}\big{(}\|\mathbf{w}_{j}-\mathbf{w}_{j}^{*}\|_{2}^{2}+|b _{j}-b_{j}^{*}|^{2}\big{)}\,.\]

Here, \(l=\big{\lceil}\frac{12}{\omega}\big{\rceil}\), and the second inequality holds in a similar way to above analysis. Then the general conclusion is handy.

## Appendix C Proof of Lemma 2 (upper bound)

For simplicity, let

\[z_{i} \coloneqq\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x} _{i}+b_{j})-\sum_{j=1}^{r}a_{j}^{*}\mathrm{relu}(\mathbf{w}_{j}^{*\top}\mathbf{x}_{i}+b_ {j}^{*}),\] \[\widehat{z}_{i}
We can bound the first term in the right hand side by

\[\frac{1}{n}\sum_{i=1}^{n}\widehat{z}_{i}^{2}\] \[=\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}a_{j}\left(\mathrm{relu }(\mathbf{w}_{j}^{\top}\mathbf{x}_{i}+b_{j})-\mathrm{relu}(\widehat{\mathbf{w}}_{j}^{\top} \mathbf{x}_{i}+\widehat{b}_{j})\right)\right]^{2}\] \[\leq\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{j=1}^{r}\left|(\mathbf{w}_{j }-\widehat{\mathbf{w}}_{j})^{\top}\mathbf{x}_{i}\right|\right]^{2}\] \[\leq\frac{r}{n}\sum_{i=1}^{n}\sum_{j=1}^{r}\left|(\mathbf{w}_{j}- \widehat{\mathbf{w}}_{j})^{\top}\mathbf{x}_{i}\right|^{2},\]

where the second line holds due to the contraction property of ReLu function, and the last line comes from the AM-GM inequality. Lemma 5 further gives for some constant \(C_{7}>0\),

\[\sum_{j=1}^{r}\frac{1}{n}\sum_{i=1}^{n}\left|(\mathbf{w}_{j}-\widehat {\mathbf{w}}_{j})^{\top}\mathbf{x}_{i}\right|^{2}\leq C_{7}\sum_{j=1}^{r}\left\|\mathbf{w}_{j}-\widehat{\mathbf{w}}_{j} \right\|_{1}^{2}\] \[+C_{7}\frac{\log^{3}\frac{p}{n\delta}}{n}\sum_{j=1}^{r}\left\| \mathbf{w}_{j}-\widehat{\mathbf{w}}_{j}\right\|_{1}^{2}\]

holds with probability at least \(1-\delta\). In addition,

\[\sum_{j=1}^{r}\left\|\mathbf{w}_{j}-\widehat{\mathbf{w}}_{j}\right\|_{1}^ {2} \leq\left\|\mathbf{W}-\widehat{\mathbf{W}}\right\|_{1}^{2}\] \[\leq\left(\left\|\mathbf{W}^{\star}\right\|_{1}-\left\|\widehat{\mathbf{ W}}\right\|_{1}\right)^{2}\leq D_{1}^{2},\]

and

\[\sum_{j=1}^{r}\left\|\mathbf{w}_{j}-\widehat{\mathbf{w}}_{j}\right\|_{2}^ {2} =\left\|\mathbf{W}-\widehat{\mathbf{W}}\right\|_{1}^{2}\leq\left\|\mathbf{W}- \widehat{\mathbf{W}}\right\|_{1}\left\|\mathbf{W}-\widehat{\mathbf{W}}\right\|_{\infty}\] \[\leq\frac{\left(\left\|\mathbf{W}^{\star}\right\|_{1}-\left\|\widehat {\mathbf{W}}\right\|_{1}\right)\left(\left\|\mathbf{W}^{\star}\right\|_{1}-\left\| \widehat{\mathbf{W}}^{\star}\right\|_{1}\right)}{S/2}\] \[\leq\frac{4}{S}D_{1}^{2}.\]

Here, \(\widehat{\mathbf{W}}^{\star}\) denote the entries of \(\widehat{\mathbf{W}}\) on the support set for \(\mathbf{W}^{\star}\), and we make use of the fact that \(\|\widehat{\mathbf{W}}\|_{1}\leq\|\mathbf{W}^{\star}\|_{1}\) and

\[\left\|\mathbf{W}-\widehat{\mathbf{W}}\right\|_{\infty}\leq\frac{\|\widehat{\mathbf{W}}^ {\star}-\widehat{\mathbf{W}}\|_{1}}{S-s}\leq\frac{\|\mathbf{W}^{\star}\|_{1}-\| \widehat{\mathbf{W}}^{\star}\|_{1}}{S/2}.\]

Putting everything together gives the desired result.

## Appendix D Technical lemmas

**Lemma 3**.: _For any \((\mathbf{W},\mathbf{a},\mathbf{b})\) with \(\left\|\mathbf{W}\right\|_{0}+\left\|\mathbf{b}\right\|_{0}+\left\|\mathbf{W}^{\star} \right\|_{0}+\left\|\mathbf{b}^{\star}\right\|_{0}\leq S\). Assume that \(\left\|\mathbf{W}\right\|_{1}+\left\|\mathbf{b}\right\|_{1}\leq\left\|\mathbf{W}^{\star} \right\|_{1}+\left\|\mathbf{b}^{\star}\right\|_{1}\) and \(\left\|\mathbf{w}_{j}^{\star}\right\|_{2}^{2}+|b_{j}^{\star}|^{2}\leq 1\). Then one has_

\[D_{1}\leq 2\sqrt{S}D_{2}, \tag{37}\]

_where \(D_{1},D_{2}\) are defined in (17)._

Proof.: For simplicity, assume that

\[D_{2}^{2}=\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_{2}^{2 }+|b_{j}-b_{j}^{\star}|^{2}\right)+\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}^{ \star}\|_{2}^{2}+|b_{j}^{\star}|^{2}\right).\]

Here, \(j\in\mathcal{J}\) means that \(a_{j}=a_{j}^{\star}\) and

\[\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_{2}^{2}+|b_{j}-b_{j}^{\star}|^{2}\leq\|\mathbf{w} _{j}^{\star}\|_{2}^{2}+|b_{j}^{\star}|^{2}.\]

Then, according to the AM-GM inequality, one has

\[\sqrt{S}D_{2} \geq\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star} \|_{1}+|b_{j}-b_{j}^{\star}|\right)\] \[\quad+\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}^{\star}\|_{1}+|b_{ j}^{\star}|\right)\] \[\geq\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}^{\star}\|_{1}-\|\mathbf{ w}_{j}\|_{1}+|b_{j}^{\star}|-|b_{j}|\right)+\|\mathbf{W}^{\star}\|_{1}\] \[\quad+\|\mathbf{b}^{\star}\|_{1}-\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_ {j}^{\star}\|_{1}+|b_{j}^{\star}|\right)\] \[\geq\sum_{j\notin\mathcal{J}}\biggl{(}\|\mathbf{w}_{j}\|_{1}+|b_{j}| \biggr{)},\]

which implies that

\[2\sqrt{S}D_{2} \geq\sum_{j\in\mathcal{J}}\left(\|\mathbf{w}_{j}-\mathbf{w}_{j}^{\star}\|_ {1}+|b_{j}-b_{j}^{\star}|\right)\] \[\quad+\sum_{j\notin\mathcal{J}}\left(\|\mathbf{w}_{j}^{\star}\|_{1}+|b_ {j}^{\star}|+\|\mathbf{w}_{j}\|_{1}+|b_{j}|\right).\]

Thus we conclude the proof. 

**Lemma 4** (Theorem 2.1 [27]).: _For any constant \(k\geq 0\), there exists some universal function \(f_{k}:\mathbb{R}\to\mathbb{R}\) such that_

\[\mathbb{E} \left[\sum_{j=1}^{r}a_{j}\mathrm{relu}(\mathbf{w}_{j}^{\top}\mathbf{x}+b_ {j})-\sum_{j=1}^{r}a_{j}^{\star}\mathrm{relu}(\mathbf{w}_{j}^{\star\top}\mathbf{x}+b_{j }^{\star})\right]^{2}\] \[=\sum_{k=0}^{\infty}\biggl{\|}\sum_{j=1}^{r}a_{j}f_{k}\biggl{(}\frac{ b_{j}}{\|\mathbf{w}_{j}\|_{2}}\biggr{)}\|\mathbf{w}_{j}\|_{2}\biggl{(}\frac{\mathbf{w}_{j}}{\|\mathbf{w}_{j}\|_{2}} \biggr{)}^{\otimes k}\] \[\qquad-\sum_{j=1}^{r}a_{j}^{\star}f_{k}\biggl{(}\frac{b_{j}^{\star}}{ \|\mathbf{w}_{j}^{\star}\|_{2}}\biggr{)}\|\mathbf{w}_{j}^{\star}\|_{2}\biggl{(}\frac{\mathbf{ w}_{j}^{\star}}{\|\mathbf{w}_{j}^{\star}\|_{2}}\biggr{)}^{\otimes k}\biggr{\|}_{F}^{2}, \tag{38}\]

_with_

\[\alpha_{k}\coloneqq f_{2k}(0)>0,\qquad\text{for all }k>0. \tag{39}\]

_In addition, we have_

\[\inf_{R_{l}}\mathbb{E}\left[a\mathrm{relu}(\mathbf{w}^{\top}\mathbf{x}+b )-\sum_{j=1}^{r}a^{\star}\mathrm{relu}(\mathbf{w}^{\star\top}\mathbf{x}+b^{\star})-R_{l }(\mathbf{x})\right]^{2}\] \[=\sum_{k>l}^{\infty}\biggl{\|}af_{k}\biggl{(}\frac{b}{\|\mathbf{w}\|_{2} }\biggr{)}\|\mathbf{w}\|_{2}\biggl{(}\frac{\mathbf{w}}{\|\mathbf{w}\|_{2}}\biggr{)}^{ \otimes k}\] \[\qquad-a^{\star}f_{k}\biggl{(}\frac{b^{\star}}{\|\mathbf{w}^{\star}\|_{2}} \biggr{)}\|\mathbf{w}^{\star}
holds with probability at least \(1-\delta\)._

Proof.: Before proceeding, we introduce some useful techniques about Restricted Isometry Property (RIP). Let \(\mathbf{X}\coloneqq\frac{1}{\sqrt{n}}[\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{n}]\). For some constant \(c_{0}>0\), if \(n\geq c_{0}\left(s\log\frac{p}{s}+\log\frac{1}{s}\right)\), then with probability at least \(1-\delta\),

\[\big{\|}\mathbf{X}^{\top}\mathbf{w}\big{\|}_{2}^{2}\leq 2\|\mathbf{w}\|_{2}^{2} \tag{42}\]

holds for all \(\mathbf{w}\) satisfying \(\|\mathbf{w}\|_{0}\leq s\).

We divide the entries of \(\mathbf{w}\) into several groups \(\mathcal{S}_{1}\cup\mathcal{S}_{2}\cup\dots\cup\mathcal{S}_{L}\) with equal size \(s\) (except for \(\mathcal{S}_{L}\)), such that the entries in \(\mathcal{S}_{j}\) are no less than \(\mathcal{S}_{k}\) for any \(j<k\). Then, according (42), one has

\[\frac{1}{n}\sum_{i=1}^{n}(\mathbf{w}^{\top}\mathbf{x}_{i})^{2} =\mathbf{w}^{\top}\mathbf{X}\mathbf{X}^{\top}\mathbf{w}=\sum_{j,k}\mathbf{w}_{\mathcal{ S}_{j}}^{\top}\mathbf{X}\mathbf{X}^{\top}\mathbf{w}_{\mathcal{S}_{k}}\] \[\leq 2\sum_{j,k}\|\mathbf{w}_{\mathcal{S}_{j}}\|_{2}\|\mathbf{w}_{ \mathcal{S}_{k}}\|_{2}=2\Big{(}\sum_{l=1}^{L}\|\mathbf{w}_{\mathcal{S}_{l}}\|_{2 }\Big{)}^{2}.\]

In addition, the order of \(\mathbf{w}_{\mathcal{S}_{l}}\) yields for \(l>1\),

\[\|\mathbf{w}_{\mathcal{S}_{l}}\|_{2} \leq\sqrt{s}\|\mathbf{w}_{\mathcal{S}_{l}}\|_{\infty}\leq\frac{1}{(l -1)\sqrt{s}}\|\mathbf{w}\|_{1},\]

which leads to

\[\Big{(}\sum_{l=1}^{L}\|\mathbf{w}_{\mathcal{S}_{l}}\|_{2}\Big{)}^{2} \leq 2\|\mathbf{w}_{\mathcal{S}_{l}}\|_{2}^{2}+2\Big{(}\sum_{l=2}^{L} \frac{1}{(l-1)\sqrt{s}}\|\mathbf{w}\|_{1}\Big{)}^{2}\] \[\leq 2\|\mathbf{w}\|_{2}^{2}+\frac{2\log^{2}L}{s}\|\mathbf{w}\|_{1}^{2}.\]

We conclude the proof by combining the above inequalities. 

## Appendix E Further experiments details

The hyper-parameters used in Section III are summarized in Table IX.

We briefly explain the variable selection procedure. We first obtain a vector of the variables' importance. For 'LASSO' and 'OMP', we use the absolute value of the estimated coefficient as the variable importance; for 'NN', 'GLASSO', and 'GSR', we obtain the importance by applying row-wise \(\ell_{2}\)-norm to the weight matrix in the input layer of the neural network; for 'RF', 'GB', and 'LNET', we use the importance produced by those methods. Once we have the importance vector, we can obtain the receiver operating characteristic (ROC) curve for synthetic datasets by varying the cut-off thresholds and calculate the AUC score. As for variable selection, we apply GMM of two mixtures to the importance vector for the synthetic datasets. The variables in the cluster with higher importance are considered significant. Then, we calculate the correctly or wrongly selected variables accordingly. For BGSBoy and UJIIndoorLoc datasets, the variables with the three- and ten-largest importance are selected, respectively.
