Sparse-Input Neural Network using Group Concave Regularization
Simultaneous feature selection and non-linear function estimation are
challenging, especially in high-dimensional settings where the number of
variables exceeds the available sample size in modeling. In this article, we
investigate the problem of feature selection in neural networks. Although the
group LASSO has been utilized to select variables for learning with neural
networks, it tends to select unimportant variables into the model to compensate
for its over-shrinkage. To overcome this limitation, we propose a framework of
sparse-input neural networks using group concave regularization for feature
selection in both low-dimensional and high-dimensional settings. The main idea
is to apply a proper concave penalty to the $l_2$ norm of weights from all
outgoing connections of each input node, and thus obtain a neural net that only
uses a small subset of the original variables. In addition, we develop an
effective algorithm based on backward path-wise optimization to yield stable
solution paths, in order to tackle the challenge of complex optimization
landscapes. Our extensive simulation studies and real data examples demonstrate
satisfactory finite sample performances of the proposed estimator, in feature
selection and prediction for modeling continuous, binary, and time-to-event
outcomes.
Bin Luo, Susan Halabi
# Sparse-Input Neural Network using Group Concave Regularization

###### Abstract

Simultaneous feature selection and non-linear function estimation are challenging, especially in high-dimensional settings where the number of variables exceeds the available sample size in modeling. In this article, we investigate the problem of feature selection in neural networks. Although the group LASSO has been utilized to select variables for learning with neural networks, it tends to select unimportant variables into the model to compensate for its over-shrinkage. To overcome this limitation, we propose a framework of sparse-input neural networks using group concave regularization for feature selection in both low-dimensional and high-dimensional settings. The main idea is to apply a proper concave penalty to the \(l_{2}\) norm of weights from all outgoing connections of each input node, and thus obtain a neural net that only uses a small subset of the original variables. In addition, we develop an effective algorithm based on backward path-wise optimization to yield stable solution paths, in order to tackle the challenge of complex optimization landscapes. Our extensive simulation studies and real data examples demonstrate satisfactory finite sample performances of the proposed estimator, in feature selection and prediction for modeling continuous, binary, and time-to-event outcomes.

Concave Regularized Neural Network

Neural networks, Feature selection, High dimensionality, LASSO, Concave penalty

## 1 Introduction

In the past decade, advancements in molecular, imaging and other laboratory tests have led to a growing interest in high-dimensional data analysis (HDDA). High-dimensional data refers to a dataset that contains a large number of observed variables relative to the small sample size, which presents a significant challenge in building accurate and interpretable models. For example, in bioinformatics, hundreds of thousands of RNA expressions, Genome-Wide Association Study (GWAS) data, and microarray data are used to understand the biology of disease, with only hundreds of patients involved (Visscher et al., 2012; Hertz et al., 2016; Kim and Halabi, 2016; Beltran et al., 2017). To address the curse of dimensionality, feature selection has become a critical step in HDDA. By identifying the most
representative features to characterize the biology of the disease or the outcome, feature selection approaches can increase the model interpretability and improve the generalization of the model.

There are various methods for feature selection, including filter methods (Koller and Sahami, 1996; Guyon and Elisseeff, 2003; Gu et al., 2012), wrapper methods (Kohavi and John, 1997; Inza et al., 2004; Tang et al., 2014), and embedded methods (Tibshirani, 1996; Zou, 2006; Fan and Li, 2001; Zhang et al., 2010). Among them, penalized regression methods have become very popular in HDDA since the introduction of the least absolute shrinkage and selection operator (LASSO) (Tibshirani, 1996). Penalized regression method can perform simultaneous parameter estimation and feature selection by shrinking some of the parameter coefficients to exact zeros. While LASSO has been widely used to obtain sparse estimations in machine learning and statistics, it tends to select unimportant variables to compensate for the over-shrinkage for relevant variables (Zou, 2006). To address the bias and inconsistent feature selection of LASSO, several methods have been proposed, including adaptive LASSO (Zou, 2006), the minimax concave penalty (MCP) (Zhang et al., 2010), and the smoothly clipped absolute deviation (SCAD) (Fan and Li, 2001).

However, most of these penalized methods assume linearity in the relationship between the variables and the outcomes, while the actual functional form of the relationship may not be available in many applications. Some additive non-parametric extensions have been proposed to resolve this problem (Lin and Zhang, 2006; Ravikumar et al., 2009; Meier et al., 2009), but their models rely on sums of univariate or low-dimensional functions and may not be able to capture the complex interactions between multiple covariates. Yamada et al. (2014) propose the HSIC-LASSO approach that leverages kernel learning for feature selection while uncovering non-linear feature interactions. However, it suffers from quadratic scaling in computational complexity with respect to the number of observations.

Neural networks are powerful tools for modeling complex relationships in a wide range of applications, from image (Krizhevsky et al., 2017; He et al., 2016) and speech recognition (Graves et al., 2013; Chan et al., 2016) to natural language processing (Young et al., 2018; Devlin et al., 2018) and financial forecasting (Fischer and Krauss, 2018). Their state-of-the-art performance has been achieved through powerful computational resources and the use of large sample sizes. Despite that, high-dimensional data can still lead to overfitting and poor generalization performance for neural networks (Liu et al., 2017). Recently, there have been novel developments in using regularized neural networks for feature selection or HDDA. A line of research focuses on utilizing the regularized neural networks, specifically employing the group LASSO technique to promote sparsity among input nodes (Liu et al., 2017; Scardapane et al., 2017; Feng and Simon, 2017). These methods consider all outgoing connections from a single input neuron as a group and apply the LASSO penalty on the \(l_{2}\) norm of weight vectors of each group. Other LASSO-regularized neural networks in feature selection can be found in the work of Li et al. (2016) and Lemhadri et al. (2021). However, regularized neural networks incorporating LASSO suffer from a tendency to over-shrink the non-zero weight of relevant variables and include many false positives in the selected model. The adaptive LASSO was employed to alleviate this problem (Dinh and Ho, 2020), yet their results are limited to continuous outcomes and assume that the conditional mean function is exactly a neural network. The work in Yamada et al. (2020) bypassed the \(l_{1}\) regularization by introducing stochastic gates to the input layer of neural networks. They considered
\(l_{0}\)-like regularization based on a continuous relaxation of the Bernoulli distribution. Their method, however, requires a cutoff value for selecting variables with weak signals, and the stochastic gate is unable to completely exclude the non-selected variables during model training and prediction stages.

In this paper, we propose a novel framework for sparse-input neural networks using group concave regularization to overcome the limitations of existing feature selection methods. Although concave penalties like MCP and SCAD have been shown to perform well in both theoretical and numerical settings for feature selection and prediction, they have not received the same level of attention as LASSO in the context of machine learning. Our proposed framework aims to draw attention to the underutilized potential of the concave penalty for feature selection in neural networks, by providing a comprehensive approach for simultaneous feature selection and function estimation in both low-dimensional and high-dimensional settings. In particular, our proposed method considers all outgoing connections from a single input neuron as a group and applies a proper concave penalty to the \(l_{2}\) norm of weights for each group. By shrinking all the weights of certain groups to exact zeros, it obtains a neural net that uses only a small subset of variables. In addition, we developed an effective algorithm based on the backward path-wise optimization that yields stable solution paths, to tackle the challenge of complex optimization landscapes. Our simulation studies and real data examples demonstrate the satisfactory finite sample performance of the group concave regularization, which outperforms existing methods in terms of feature selection and prediction accuracy for modeling continuous, binary, and time-to-event outcomes.

The rest of this article is organized as follows. In Section 2, we formulate the problem of feature selection for a generic non-parametric model and introduce our proposed method. The implementation of the method, including the composite gradient descent algorithm and the backward path-wise optimization, is presented in Section 3. In Section 4, we conduct extensive simulation studies to demonstrate the performance of the proposed method. The application of the method to real-world datasets is presented in Section 5. Lastly, in Section 6, we discuss the results and their implications.

## 2 Method

### Problem setup

Let \(X\in\mathbb{R}^{d}\) be a \(d\)-dimensional random vector and \(Y\) be a response variable. We assume the conditional distribution \(P_{Y|X}\) depends on a form of \(f(X_{S})\) with a function \(f\in F\) and a subset of variables \(S\subseteq\{1,\cdots,d\}\). We are interested in identifying the true set \(S\) for significant variables and estimating function \(f\) so that we can predict \(Y\) based on selected variable \(X_{S}\).

At the population level, we aim to minimize the loss

\[\min_{f\in F,S}\mathbb{E}_{X,Y}\ell(f(X_{S}),Y)\]

where \(\ell\) is a loss function tailored to a specific problem. In practical settings, the distribution of \((X,Y)\) is often unknown, and instead only an independent and identically distributed (i.i.d.) random sample of size \(n\) is available, consisting of pairs of observations \(\left(X_{i},Y_{i}\right)_{i=1}^{n}\). Additionally, if the number of variables \(d\) is large, an exhaustive search over all possible
subsets \(S\) becomes computationally infeasible. Furthermore, we do not assume any specific form of the unknown function \(f\) and aim to approximate \(f\) nonparametrically using neural networks. Thus, our goal is to develop an efficient method that can simultaneously select a variable subset \(S\) and approximate the solution \(f\) for any given class of functions using a sparse-input neural network.

### Proposed framework

We consider function estimators based on feedforward neural networks. Let \(\mathcal{F}_{n}\) be a class of feed forward neural networks \(f_{\mathbf{w}}:\mathbb{R}^{d}\mapsto\mathbb{R}\) with parameter \(\mathbf{w}\). The architecture of a multi-layer perceptron (MLP) can be expressed as a composition of a series of functions

\[f_{\mathbf{w}}(x)=L_{D}\circ\sigma\circ L_{D-1}\circ\sigma\circ\cdots\circ \sigma\circ L_{1}\circ\sigma\circ L_{0}(x),x\in\mathbb{R}^{d},\]

where \(\circ\) denotes function composition and \(\sigma(x)\) is an activation function defined for each component of \(x\). Additionally,

\[L_{i}(x)=\mathbf{W}_{i}x+b_{i},i=0,1,\ldots,\mathcal{D},\]

where \(\mathbf{W}_{i}\in\mathbb{R}^{d_{i+1}\times d_{i}}\) is a weight matrix, \(D\) is the number of hidden layers, \(d_{i}\) is the width defined as the number of neurons of the \(i\)-th layer with \(d_{0}=d\), and \(b_{i}\in\mathbb{R}^{d_{i+1}}\) is the bias vector in the \(i\)-th linear transformation \(L_{i}\). Note that the vector \(\mathbf{w}\in\mathbb{R}^{P}\) is the column-vector concatenation of all parameters in \(\{\mathbf{W}_{i},b_{i}:i=0,1,\ldots,\mathcal{D}\}\). We define the empirical loss of \(f_{\mathbf{w}}\) as

\[\mathcal{L}_{n}(\mathbf{w})=\frac{1}{n}\sum_{i=1}^{n}\ell(f_{\mathbf{w}}(X_{i} ),Y_{i}).\]

The ideal scenario is to have a sparse-input neural network \(f_{\mathbf{w}}\) that only takes signals from the important variables, meaning that \(\mathbf{W}_{0,j}=\mathbf{0}\) for \(j\notin S\), where \(\mathbf{W}_{0,j}\) denotes the \(j\)th column vector of \(\mathbf{W}_{0}\). In order to minimize the empirical loss \(\mathcal{L}_{n}(\mathbf{w})\) while inducing sparsity in \(\mathbf{W}_{0}\), we propose to train the neural network by minimizing the following group regularized empirical loss

\[\hat{\mathbf{w}}=\operatorname*{argmin}_{\mathbf{w}\in\mathbb{R}^{P}}\left\{ \mathcal{L}_{n}(\mathbf{w})+\sum_{j=1}^{d}\rho_{\lambda}(\|\mathbf{W}_{0,j}\| _{2})+\alpha\|\mathbf{w}\|_{2}^{2}\right\}, \tag{1}\]

where \(\|\cdot\|_{2}\) denotes the Euclidean norm of a vector.

The objective function in Eq. (1) comprises three components:

1. \(\mathcal{L}_{n}(\mathbf{w})\) is the empirical loss function such as the mean squared error loss for regression tasks, the cross-entropy loss for classification tasks, and the negative log partial likelihood for proportional hazards models. Further details can be found in Appendix A.
2. \(\rho_{\lambda}\) is a concave penalty function parameterized by \(\lambda\geq 0\). To simultaneously select variables and learn the neural network, we group the outgoing connections from each single input neuron that corresponds to each variable. The concave penalty function \(\rho_{\lambda}\) is designed to shrink the weight vectors of specific groups to exact zeros, resulting in a neural network that utilizes only a small subset of the original variables.
3. \(\alpha\|\mathbf{w}\|_{2}^{2}\), where \(\alpha>0\), represents the ridge regularization term used to prevent overfitting in neural networks. Note that feature selection, employing \(\rho_{\lambda}\), depends exclusively on the magnitudes of weights in the input layer. However, it is possible to diminish the influence of \(\rho_{\lambda}\) by reducing all weights in the input layer while simultaneously allowing larger weights in other layers, without affecting the network's output. The ridge regularization addresses this issue by promoting smaller, well-balanced weights, thereby improving model stability and mitigating overfitting.

Note that when the number of hidden layers \(D=0\), the function \(f_{\mathbf{w}}\) reduces to a linear function, and the optimization problem in Eq. (1) becomes the framework of elastic net (Zou and Hastie, 2005), SCAD-\(L_{2}\)(Zeng and Xie, 2014), and Mnet (Huang et al., 2016), with the choice of \(\rho_{\lambda}\) to be LASSO, SCAD, and MCP, respectively.

### Concave regularization

There are several commonly used penalty functions that encourage sparsity in the solution, such as LASSO (Tibshirani, 1996), SCAD (Fan and Li, 2001), and MCP (Zhang et al., 2010). When applied to the \(l_{2}\)-norm of the coefficients associated with each group of variables, these penalty functions give rise to group regularization methods, including group LASSO (GLASSO) (Yuan and Lin, 2006), group SCAD (GSCAD) (Guo et al., 2015), and group MCP (GMCP) (Huang et al., 2012). Specifically, LASSO, SCAD, and MCP are defined as follows.

* **LASSO** \[\rho_{\lambda}(t)=\lambda|t|.\]
* **SCAD** \[\rho_{\lambda}(t)=\begin{cases}\lambda|t|&\text{for }|t|\leq\lambda,\\ -\frac{t^{2}-2a\lambda|t|+\lambda^{2}}{2(a-1)}&\text{for }\lambda<|t|\leq a \lambda,\\ \frac{(a+1)\lambda^{2}}{2}&\text{for }|t|>a\lambda,\end{cases}\] where \(a>2\) is fixed.
* **MCP** \[\rho_{\lambda}(t)=\operatorname{sign}(t)\lambda\int_{0}^{|t|}\left(1-\frac{z} {\lambda a}\right)_{+}dz,\] where \(a>0\) is fixed.

It has been demonstrated, both theoretically and numerically, that the concave regularization methods of SCAD and MCP exhibit strong performance in terms of feature selection and prediction (Fan and Li, 2001; Zhang et al., 2010). Unlike the convex penalty LASSO, which tends to over-regularize large terms and provide inconsistent feature selection, concave regularization can reduce LASSO's bias and improve model selection accuracy. The rationale behind the concave penalty lies in the behavior of its derivatives. Specifically, SCAD and MCP initially apply the same level of penalization as LASSO, but gradually reduce the penalization rate until it drops to zero when \(t>a\lambda\). Given the benefits of the concave penalization, we propose using the group concave regularization in our framework for simultaneous feature selection and function estimation.
## 3 Implementation

### Composite gradient descent

Note that the optimization in Eq. (1) is not a convex optimization problem since both empirical loss function \(\mathcal{L}_{n}(\mathbf{w})\) and the penalty function \(\rho_{\lambda}\) can be non-convex. To obtain the stationary point, we use the composite gradient descent algorithm (Nesterov, 2013). This algorithm is also incorporated in Feng and Simon (2017); Lembadri et al. (2021) for sparse-input neural networks based on the LASSO regularization.

Denote \(\bar{\mathcal{L}}_{n,\alpha}(\mathbf{w})=\mathcal{L}_{n}(\mathbf{w})+\alpha\| \mathbf{w}\|_{2}^{2}\) as the smooth component of the objective function in Eq. (1). The composition gradient iteration for epoch \(t\) is given by

\[\mathbf{w}^{t+1}=\operatorname*{argmin}_{\mathbf{w}}\left\{\frac{1}{2}\| \mathbf{w}-\tilde{\mathbf{w}}^{t+1}\|_{2}^{2}+\sum_{j=1}^{d}\rho_{\lambda}(\| \mathbf{W}_{0,j}\|_{2})\right\}, \tag{2}\]

where \(\tilde{\mathbf{w}}^{t+1}=\mathbf{w}^{t}-\gamma\nabla\bar{\mathcal{L}}_{n, \alpha}(\mathbf{w}^{t})\) is the gradient update only for the smooth component \(\bar{\mathcal{L}}_{n,\alpha}(\mathbf{w}^{t})\) that can be computed using the standard back-propagation algorithm. Here \(\gamma>0\) is the learning rate for the update and can be set as a fixed value or determined by employing the backtracking line search method, as described in Nesterov (2013). Let \(A_{j}\) represent the index set of \(\mathbf{W}_{0,j}\) within \(\mathbf{w}\). We define \(A\) as the index set that includes all weights in the input layer, given by \(A=\{\bigcup_{j=1}^{d}A_{j}\}\). By solving Eq. (2), we obtain the iteration form \(\mathbf{w}_{A^{c}}^{t+1}=\tilde{\mathbf{w}}_{A^{c}}^{t+1}\) and

\[\mathbf{w}_{A_{j}}^{t+1}=h(\tilde{\mathbf{w}}_{A_{j}}^{t+1},\lambda),\quad \text{for }j=1,\cdots,d. \tag{3}\]

Here, \(A^{c}\) refers to the complement of the set \(A\), and the function \(h\) represents the thresholding operator, which can be determined by the specific penalty \(\rho_{\lambda}\). By taking \(\rho_{\lambda}\) to be the LASSO, MCP, and SCAD penalty, it can be verified that the GLASSO, GSCAD, and GMCP solutions for the iteration in Eq. (3) have the following form:

* GLASSO \[h_{\text{GLASSO}}(z,\lambda)=S_{\text{g}}(z,\lambda).\]
* GSCAD \[h_{\text{GSCAD}}(z,\lambda)=\begin{cases}S_{\text{g}}(z,\lambda),&\text{if } \|z\|_{2}\leq 2\lambda,\\ \frac{a-1}{a-2}S_{\text{g}}(z,\frac{a\lambda}{a-1}),&\text{if }2\lambda<\|z\|_{2} \leq a\lambda,\\ z,&\text{if }\|z\|_{2}>a\lambda.\end{cases}\]
* GMCP \[h_{\text{GMCP}}(z,\lambda)=\begin{cases}\frac{a}{a-1}S_{\text{g}}(z,\lambda),& \text{if }\|z\|_{2}\leq a\lambda,\\ z,&\text{if }\|z\|_{2}>a\lambda,\end{cases}\]

where \(S_{\text{g}}(z;\lambda)\) is the group soft-thresholding operator defined as

\[S_{\text{g}}(z;\lambda)=\left(1-\frac{\lambda}{\|z\|_{2}}\right)_{+}z.\]

Therefore, we can efficiently implement the composite gradient descent by integrating an additional thresholding operation into the input layer. This operation follows the gradient
descent step using the smooth component \(\bar{\mathcal{L}}_{n,\alpha}(\mathbf{w})\). The calculation for epoch t can be summarized as follows:

1. compute gradient of \(\bar{\mathcal{L}}_{n,\alpha}(\mathbf{w}^{t})\) using back-propagation,
2. update \(\tilde{\mathbf{w}}^{t+1}\leftarrow\mathbf{w}^{t}-\gamma\nabla\bar{\mathcal{L} }_{n,\alpha}(\mathbf{w}^{t})\),
3. update \(\mathbf{w}^{t+1}_{A^{c}}\leftarrow\tilde{\mathbf{w}}^{t+1}_{A^{c}}\) and \(\mathbf{w}^{t+1}_{A_{j}}\gets h(\tilde{\mathbf{w}}^{t+1}_{A_{j}},\lambda)\), for \(j=1,\cdots,d\).

The final index set of the selected variables is \(\hat{S}=\{j:\hat{\mathbf{w}}_{A_{j}}\neq\mathbf{0}\}\).

### Backward path-wise optimization

We are interested in learning neural networks not only for a specific value of \(\lambda\), but also for a range of \(\lambda\) where the networks vary by the number of included variables. Specifically, we consider a range of \(\lambda\) from \(\lambda_{min}\), where the networks include all or an excessively large number of variables, up to \(\lambda_{max}\), where all variables are excluded and \(|W_{0}|\) becomes a zero matrix. Since the objective function is not convex and has multiple local minima, the solution of Eq. (1) with random initialization may not vary continuously for \(\lambda\in[\lambda_{min},\lambda_{max}]\), resulting in a highly unstable path of solutions that are regularized by \(\lambda\).

To address this issue, we consider a path-wise optimization strategy by varying the regularization parameter along a path. In this approach, we use the solution of a particular value of \(\lambda\) as a warm start for the next problem. Regularized linear regression methods (Friedman et al., 2007, 2010; Breheny and Huang, 2011) typically adopt a forward path-wise optimization, starting from a null model with all variables excluded at \(\lambda_{max}\) and working forward with decreasing \(\lambda\)s. However, our numerical studies for sparse-input neural networks showed that starting from a sparse solution as an initial model does not produce a larger model along the path until jumping to the full model at a sufficiently small \(\lambda\). To tackle this problem, we implement a backward path-wise optimization approach, starting from a dense model at the minimum value of \(\lambda_{min}\) and solving toward sparse models up to \(\lambda_{max}\) with all variables excluded from the network. This dense-to-sparse warm start approach is also employed in (Lembadri et al., 2021) using LASSO regularization.

To further illustrate the importance of using backward path-wise optimization in regularized neural networks, we investigate variables selection and function estimation of a regression model \(Y=f(X)+\epsilon\), where \(f(X)=\log(|X_{1}|+0.1)+X_{1}X_{2}+X_{2}+\exp(X_{3}+X_{4})\) with 4 informative and 16 nuisance variables, and each \(X_{i}\) and \(\epsilon\) follow the standard normal distribution. More details of the simulations are presented in Section 4. Figure 1 shows the solution paths of GMCP and GLASSO based on different types of optimization. It is observed that non-pathwise optimization leads to fluctuations or variations in the solution path, whereas forward path-wise optimization tends to remain in the same sparse model until transitioning to a full model with a sufficiently small \(\lambda\). In contrast, backward path-wise optimization using GMCP and GLASSO produces relatively smooth solution paths. It should be noted that GLASSO has a tendency to over-shrink the weight vectors of informative variables and include more variables in the model. In contrast, GMCP is designed to prevent over-shrinkage and offers a smooth transition from the full model to the null model.

In addition to providing stable and smooth solution paths, backward path-wise optimization is advantageous computationally. In particular:
* The consecutive estimates of weights in the path are close, which reduces the rounds of gradient descent needed for each iteration. Therefore, the bulk of the computational cost occurs at \(\lambda_{min}\), and a lower number of iterations for the remaining \(\lambda\)s results in low computational costs.
* We observe that the excluded variables from previous solutions are rarely included in the following solutions. By pruning the inputs of the neural network along the solution path, further reduction in computation complexity can be achieved as the model becomes sparse. Since the computational cost scales with the number of input features, this approach can significantly speed up computation, particularly for high-dimensional data.

### Tuning Parameter Selection

Two tuning parameters are required in our proposed framework: the group penalty coefficient \(\lambda\) and the ridge penalty coefficient \(\alpha\). The former controls the number of selected variables and yields sparser models for larger values of \(\lambda\), while the latter imposes a penalty on the size of the network weights to prevent overfitting.

In all numerical studies presented in this paper, we adopted a 20% holdout validation set from the training data. The model was trained using the remaining data, and the optimal values for \(\lambda\) and \(\alpha\) were selected from a fine grid of values based on their performance on the validation set.

Python code and examples for the proposed group concave regularized neural networks are available at [https://github.com/r08in/GCRNN](https://github.com/r08in/GCRNN).

## 4 Simulation Studies

We assess the performance of the proposed regularized neural networks in feature selection and prediction through several simulation settings. The data are generated through the following function:

\[f(X)=\log(|X_{1}|+0.1)+X_{1}X_{2}+X_{2}+\exp(X_{3}+X_{4}),\]

where each component of the covariate vector \(X=(X_{1},\cdots,X_{d})^{T}\in\mathbb{R}^{d}\) are generated from independent standard normal distribution. Here \(d>4\) and function \(f(X)\) is sparse that only the first four variables are relevant to the outcome. We generate \(n\) i.i.d. random samples with continuous outcomes, binary outcomes, and time-to-event outcomes in the following three examples, respectively.

**Example 1** (Regression Model): _The continuous response \(Y\) is generated from a standard regression model with an additive error as follows_

\[Y=f(X)+\epsilon,\]

_where \(\epsilon\) follows a standard normal distribution._
Figure 1: **Solution path of \(l_{2}\) norm of the weight vector associated with each input node \(\|\mathbf{W}_{0j}\|_{2}\). Top left: Non-pathwise optimization using GMCP. All the neural network weights are initialized by drawing from \(N(0,0.1)\) for each \(\lambda\). Top right: forward path-wise optimization using GMCP. It starts from the null model and computes the solution with decreasing \(\lambda\). Random initialization is used before the selection of the first set of variables. Bottom left: backward path-wise optimization using GMCP. Bottom right: backward path-wise optimization using GLASSO.**
**Example 2**: **(Classification Model)** _The binary response \(Y\in\{0,1\}\) is generated from a Bernoulli distribution with the following conditional probability_

\[P(Y=1|X)=\frac{1}{1+\exp{(-f(X))}}.\]

**Example 3**: **(Proportional Hazards Model)**_The survival time \(T\) follows the proportional hazards model with a hazard function of the form_

\[h(t|X)=h_{0}(t)\exp{(f(X))}, \tag{4}\]

_where \(h_{0}(t)\) is the baseline hazard function. Thus, \(T=H_{0}^{-1}\left(-\log(U)\exp{(f(X))}\right)\), where \(U\) is a uniform random variable in \([0,1]\), and \(H_{0}\) is the baseline cumulative hazard function defined as \(H_{0}(t)=\int_{0}^{t}h_{0}(u)du\). We considered a Weibull hazard function for \(H_{0}\), with the scale parameter \(=2\) and the shape parameter \(=2\). Among \(n\) samples, \(\mathcal{C}\times n\) of them are randomly chosen as censoring observations with the censoring indicator \(\delta_{i}=0\) and otherwise \(\delta_{i}=1\) for event observations. The censoring rate \(\mathcal{C}=0\), \(0.2\) and \(0.4\) in our simulation. We define the observed time_

\[Y_{i}=\begin{cases}T_{i}&\text{if }\delta_{i}=1,\\ C_{i}&\text{if }\delta_{i}=0,\end{cases}\]

_where censoring time \(C_{i}\) is drawn from a uniform distribution \((0,T_{i})\)._

For each example, we consider the low and high dimensional settings in the following scenarios:

1. Low dimension (LD): \(d=20\) and \(n=300\) and \(500\).
2. High dimension (HD): \(d=1000\) and \(n=500\).

We perform \(200\) simulations for each scenario. The performance of the trained model in prediction and feature selection are evaluated on independently generated \(n\) random samples by the following measures:

1. Prediction score, which is defined as the \(R^{2}\) score, classification accuracy, and C-index for the regression, classification, and proportional hazards model, respectively.
2. Model size (MS), is the average number of selected covariates.
3. False positives rate (FPR), is the percent of selected but unimportant covariates: \[FPR=\frac{|\hat{S}\bigcap S^{c}|}{|S^{c}|}\times 100\%.\]
4. False negatives rate (FNR), is the percent of non-selected but important covariates: \[FNR=\frac{|\hat{S}^{c}\bigcap S|}{|S|}\times 100\%.\]
Recall that \(S\) represents the true index sets of important variables and \(\hat{S}=\{j:\|\mathbf{W}_{0,j}\|_{2}\neq 0\}\) denote the index sets of selected variables.

In our numerical studies, we consider the concave regularization GMCP and GSCAD for our proposed framework. We name the method of regularized neural networks using GLASSO, GMCP, and GSCAD as GLASSONet, GMCPNet, and GSCADNet, respectively. We compare the proposed group concave regularized estimator GMCPNet and GSCADNet with GLASSONet, neural network (NN) without feature selection (\(\lambda=0\)), random (survival) forest (RF), and the STG method proposed in Yamada et al. (2020). We also include the oracle version of NN and RF (Oracle-NN and Oracle-RF) as benchmarks, where true relevant variables are known in advance and used directly in the model fitting process. See Appendix D for the implementation details.

### Results

Table 1 presents a summary of the feature selection performance of the four approaches, namely STG, GLASSONet, GMCPNet, and GSCADNet, across all simulation scenarios. We exclude the results of the STG method for Example 3 as it either selects all variables or none of them for the survival outcome. For both LD and HD settings, GMCPNet and GSCADNet consistently outperform the STG and GLASSONet in terms of feature selection. These models exhibit superior performance, achieving model sizes that closely matched the true model, along with low false positive rates (FPR) and false negative rates (FNR) for most scenarios. While STG performs well in certain LD settings, it tends to over select variables in HD scenarios with a large variability in the model size. On the other hand, GLASSONet is prone to selecting more variables, leading to larger model sizes in both LD and HD settings, which aligns with the inherent nature of the LASSO penalty.

Figure 2 displays the distribution of testing prediction scores for the regression, classification, and proportional hazards models (PHM) with a censoring rate of \(\mathcal{C}=0.2\). The complete results of the PHM can be found in Appendix B. GMCPNet and GSCADNet demonstrate comparable performance in both LD and HD settings, achieving similar results to the Oracle-NN and outperforming NN, RF, and even Oracle-RF in most of scenarios. STG performs similarly to Oracle-NN in the LD setting of the regression model, but its performance deteriorates in the HD setting and other models. Conversely, while GLASSONet outperforms or is comparable to the Oracle-RF method in the LD settings, it suffers from overfitting in the HD settings by including a large number of false positives in the final model.

It is worth pointing out that the Oracle-NN outperforms the Oracle-RF in every scenario, indicating that neural network-based methods can serve as a viable alternative to tree-based methods when the sample size is sufficiently large relative to the number of predictors. Additionally, NN without feature selection performs the worst across all the simulation scenarios, highlighting the importance of feature selection, especially in the high-dimensional space.

Overall, the simulation results demonstrate the superior performance of the concave penalty in terms of feature selection and prediction. The proposed GMCPNet and GSCADNet methods exhibit remarkable capabilities in selecting important variables with low FPR and low FNR, while achieving accurate predictions across various models. These meth
ods show promise for tackling the challenges of feature selection and prediction in high-dimensional data.

## 5 Real Data Example

### Survival Analysis on CALGB-90401 dataset

We utilize the data from the CALGB-90401 study, a double-blinded phase III clinical trial that compares docetaxel and prednisone with or without bevacizumab in men with metastatic castration-resistant prostate cancer (mCRPC) to illustrate the performance of our proposed method. The CALGB-90401 data consists of 498,801 single-nucleotide polymorphisms (SNPs) that are processed from blood samples from patients. We assume a dominant model for SNPs and thus each of the SNPs is considered as a binary variable. Since most SNPs are irrelevant for predicting patient survival, we only consider 181 SNPs that are associated with DNA damage-repair genes, and 444 prioritized SNPs based on an updated literature search (Mateo et al., 2015; Wyatt et al., 2016; Beltran et al., 2011;

\begin{table}
\begin{tabular}{|l|l|l l|l l|l l|} \hline \multirow{2}{*}{Model} & \multirow{2}{*}{Method} & \multicolumn{2}{l|}{\(n=300\), \(d=20\)} & \multicolumn{2}{l|}{\(n=500\), \(d=20\)} & \multicolumn{2}{l|}{\(n=500\), \(d=1000\)} \\  & & FPR, FNR & MS (SD) & FPR, FNR & MS (SD) & FPR, FNR & MS (SD) \\ \hline \multirow{4}{*}{Regression} & STG & 7.8, 5.4 & 5.0 (2.0) & 7.2, 2.1 & 5.1 (1.7) & 1.6, 12.1 & 19.2 (28.0) \\  & GLASSONet & 86.7, 4.4 & 17.7 (4.7) & 96.0, 0.6 & 19.3 (2.2) & 24.3, 29.2 & 245.0 (98.7) \\  & GMCPNet & 2.2, 4.5 & 4.2 (1.0) & 2.1, 4.2 & 4.2 (1.0) & 0.0, 5.8 & 4.1 (0.9) \\  & GSCADNet & 2.4, 5.0 & 4.2 (1.1) & 2.0, 3.2 & 4.2 (0.9) & 0.0, 7.1 & 4.1 (1.0) \\ \hline \multirow{4}{*}{Classification} & STG & 25.3, 16.5 & 7.4 (6.9) & 10.1, 11.0 & 5.2 (4.8) & 3.8, 15.6 & 40.9 (183.4) \\  & GLASSONet & 89.2, 1.0 & 18.2 (2.8) & 94.7, 0.2 & 19.1 (2.0) & 16.3, 21.5 & 165.4 (92.9) \\  & GMCPNet & 14.4, 3.9 & 6.2 (3.6) & 9.3, 0.8 & 5.5 (2.6) & 0.3, 16.2 & 6.5 (4.2) \\  & GSCADNet & 11.6, 5.8 & 5.6 (2.9) & 7.0, 1.0 & 5.1 (1.9) & 0.3, 16.8 & 6.8 (5.9) \\ \hline \multirow{4}{*}{Survival (\(\mathcal{C}=0\))} & GLASSONet & 97.2, 0.0 & 19.5 (1.0) & 99.2, 0.0 & 19.9 (0.5) & 18.2, 20.0 & 184.8 (56.2) \\  & GMCPNet & 1.6, 0.4 & 4.2 (0.6) & 0.8, 0.0 & 4.1 (0.4) & 0.0, 1.5 & 4.1 (0.5) \\  & GSCADNet & 1.9, 0.2 & 4.3 (0.6) & 1.2, 0.0 & 4.2 (0.5) & 0.0, 1.6 & 4.1 (0.7) \\ \hline \multirow{4}{*}{Survival (\(\mathcal{C}=0.2\))} & GLASSONet & 98.0, 0.1 & 19.7 (0.8) & 99.6, 0.0 & 19.9 (0.3) & 16.8, 18.0 & 170.6 (49.0) \\  & GMCPNet & 1.9, 0.4 & 4.3 (0.9) & 1.7, 0.0 & 4.3 (1.0) & 0.0, 2.6 & 4.2 (0.9) \\  & GSCADNet & 1.8, 0.2 & 4.3 (0.9) & 1.7, 0.1 & 4.3 (0.8) & 0.0, 3.5 & 4.1 (0.7) \\ \hline \multirow{4}{*}{Survival (\(\mathcal{C}=0.4\))} & GLASSONet & 95.0, 0.0 & 19.2 (1.7) & 98.8, 0.0 & 19.8 (0.5) & 15.2, 19.9 & 154.6 (48.4) \\  & GMCPNet & 5.8, 8.1 & 4.6 (1.5) & 1.2, 0.1 & 4.2 (0.5) & 0.0, 4.2 & 4.1 (1.0) \\ \cline{1-1}  & GSCADNet & 4.8, 7.5 & 4.5 (1.3) & 1.7, 0.0 & 4.3 (0.7) & 0.0, 4.9 & 4.2 (1.0) \\ \hline \end{tabular}
\end{table}
Table 1: **Feature selection results of STG, GLASSONet, GMCPNet, and GSCADNet under the regression, classification, and proportional hazards models.** The false positives rate (FPR %), false negatives rate (FNR %), and model size (MS) with standard deviation (SD) in parentheses are displayed.
Figure 2: **Top row**: \(R^{2}\) score of the proposed methods for the regression model outlined in Example 1. **Middle row**: Accuracy of the proposed methods for the classification model outlined in Example 2. **Bottom row**: C-Index of the proposed methods for the survival model outlined in Example 3. The dashed lines represent the median score of the Oracle-NN, used as a benchmark for comparison.
Mosquera et al., 2013; Robinson et al., 2015; Abida et al., 2019; De Laere et al., 2017). We also include the eight clinical variables that have been identified as prognostic markers of overall survival in patients with mCRPC (Halabi et al., 2014): opioid analgesic use (PAIN), ECOG performance status, albumin (ALB), disease site (defined as lymph node only, bone metastases with no visceral involvement, or any visceral metastases), LDH greater than the upper limit of normal (LDH.High), hemoglobin (HGB), PSA, and alkaline phosphatase (ALKPHOS). The final dataset has \(d=635\) variables with a number of patients \(n=631\) and censoring rate \(C=6.8\%\).

We consider the proportional hazard model in the form of Eq. (4) for our proposed methods to identify clinical variables or SNPs that can predict the primary outcome of overall survival in these patients. To evaluate the feature selection and prediction performance of the methods, we randomly split the dataset 100 times into training sets (n=526) and testing sets (n=105) using a 5:1 allocation ratio. We apply the methods to each of the training sets and calculate the time-dependent area under the receiver operating characteristic curve (tAUC) on the corresponding testing sets. The tAUC assesses the discriminative ability of the predicted model and is computed using the Uno method (Uno et al., 2007). The results of the 100 random splits are presented in Figure 3. Our proposed method, GSCADNet, outperforms the others in survival prediction (left panel). It is worth noting that the NN method, which lacks feature selection, tends to overfit in high-dimensional data and performs poorly. Although these three regularized methods of sparse-input neural networks perform similarly in survival prediction, GLASSONet has a tendency to over-select variables and the proposed GMCPNet and GSCADNet select a relatively smaller set of variables without compromising prediction performance (middle panel). The right panel of Figure 3 demonstrates that GSCADNet successfully selects most of the significant clinical variables and detects some of the important SNPs in predicting overall survival.

Figure 3: **Left**: Boxplots of tAUC from testing set over 100 random splits. **Middle**: the number of selected variables for GLASSONet, GMCPNet, and GSCADNet. **Right**: Variables selected by GSCADNet with selection proportion\(\geq 10\%\) over 100 random splits.
### Classification on MNIST Dataset

We aim to visualize the selection of variables by considering the classification problem on the MNIST dataset. The MNIST dataset is a well-known benchmark dataset in computer vision, consisting of grayscale images of handwritten digits from 0 to 9. In this study, we focus on the binary classification problem of distinguishing between 7s and 8s in the MNIST dataset. We evaluate our proposed methods GMCPNet and GSCADNet, along with existing methods GLASSONet, STG, NN, and RF, based on their feature selection and classification accuracy.

The MNIST dataset consists of grayscale images with 28x28 pixels, which gives 784 variables. We randomly select 250 pictures of 7s and 8s from the MNIST dataset, respectively, to form a high-dimensional training dataset with \(d=784\) and \(n=500\). Note that the class labels depend only on the pixels of the central area of the images, and thus a good method for feature selection should identify the relevant pixels and classify the images of 7s and 8s. We also corrupt the images with i.i.d. random noise from a standard normal distribution so that the input features are not sparse. The trained models are evaluated on the testing dataset with 2002 images. We repeated the process of random sampling and model fitting 100 times, and the feature (pixel) selection and classification results are shown in Figure 4. We observe that GLASSONet, GMCPNet, GSCADNet all achieve median accuracies greater than 91% and outperform the other methods. While the heatmaps of feature selection show that GLASSONet, GMCPNet, GSCADNet consistently select relevant pixels in high frequencies, GLASSONet tends to over select variables and GMCPNet and GSCADNet choose irrelevant pixels in much lower frequencies (indicated by dark red colors).

## 6 Discussion

Among the plethora of feature selection methods, penalized regression has gained significant popularity. However, many of these methods rely on the assumption and application of linear theory, which may not capture the complex relationships between covariates and the outcome of interest. In biomedical research, for instance, researchers often normalize data and employ penalized techniques under a linear model for feature selection. However, relying solely on data transformation risks overlooking intricate biological relationships and fails to address the dynamic nature of on-treatment biomarkers. Moreover, advancements in molecular and imaging technologies have introduced challenges in understanding the non-linear relationships between high-dimensional biomarkers and clinical outcomes. Novel approaches are urgently needed to tackle these complexities, leading to an improved understanding of non-linear relationships and optimizing patient treatment and care.

In this paper, we have proposed a novel framework that utilizes group concave regularization for feature selection and function estimation in complex modeling, specifically designed for sparse-input neural networks. Unlike the convex penalty LASSO, the concave regularization methods such as MCP and SCAD gradually reduce the penalization rate for large terms, preventing over-shrinkage and improving model selection accuracy. Our optimization algorithm, based on the composite gradient descent, is simple to implement, requiring only an additional thresholding operation after the regular gradient descent step on the smooth component. Furthermore, we incorporate backward path-wise optimization to efficiently navigate the optimization landscape across a fine grid of tuning parameters,
Figure 4: **Comparing feature selection and classification performance by STG, GLASSONet, GMCPNet and GSCADNet. Top left**: the image that takes the average of all the images in the training set and shows relevant pixels in grayscale. **Bottom right**: testing accuracy for clarification of 7s and 8s in the MNIST dataset using training dataset \(d=784\) and \(n=500\). **Other panels**: heatmaps depicting the selection frequencies of each pixel across 100 repetitions for each method. Lighter colors indicate higher selection frequencies, with white highest and darker colors lowest.
generating a smooth solution path from dense to sparse models. This path-wise optimization approach improves stability and computational efficiency, potentially enhancing the applicability of our framework for sparse-input neural networks.

The runtime of our proposed method over a solution path of \(\lambda\)s (with a fixed \(\alpha\)) can be comparable to or even shorter than training a single model with a fixed \(\lambda\), such as the NN method without feature selection (\(\lambda=0\)). To illustrate this, we examine the algorithm complexity of the NN method, which can be approximated as \(\mathcal{O}(ndT)\), where \(T\) denotes the number of epochs for learning the neural network. In contrast, training our proposed method over a solution path of \(m\)\(\lambda\)s has a complexity of \(\mathcal{O}(n\bar{d}T^{\prime}m)\), where \(\bar{d}\) represents the averaged number of inputs along the solution path with dimension pruning, and \(T^{\prime}\) is the number of epochs for each \(\lambda\) in the path. In our simulation with the HD scenario (\(d=1000\)), we set \(T=5000\), \(T^{\prime}=200\), and \(m=50\). Assuming the number of inputs decreases equally along the solution path from the full model to the null model, we have \(\bar{d}=d/2=500\). Thus, \(ndT=n\bar{d}T^{\prime}m\) indicates that solving for an entire path of our proposed method requires a similar computation as training a single model. In real applications, especially in high-dimensional scenarios, the dimensionality usually drops quickly along the solution path. Therefore, \(\bar{d}\) can be much smaller than \(d/2\), and thus solving for a whole solution path can be more computationally efficient. It is worth pointing out that we set \(T^{\prime}\) to be small for the first parameter \(\lambda_{\min}\) as well in the HD setting, to avoid overfitting of an initial dense model.

In our numerical studies, the parameter tunings are limited to \(\lambda\) and \(\alpha\). However, in real-world applications, it may be necessary to tune additional hyperparameters, such as the learning rate, the number of layers, and the number of nodes in each layer. The computation cost associated with tuning these parameters can be reduced by leveraging parallel computing techniques. Furthermore, when the sample size is moderate and the important variables are sparse, we have observed that using a two- or three-layer neural network with a modest number of nodes per layer (e.g., 5 or 10 nodes per layer) is often sufficient for a wide range of datasets.

One limitation of the proposed method arises in ultra-high dimensional scenarios where the number of variables reaches hundreds of millions. Directly applying the proposed sparse-input neural networks in such cases can lead to an exceedingly complex optimization landscape, making it computationally infeasible. To mitigate this limitation, one suggestion is to employ a pre-screening method to reduce the dimensionality to a more manageable size prior to applying the proposed approach.

Another limitation pertains to the proposed group regularized method, which is primarily focused on individual feature selection. This limitation becomes particularly relevant when dealing with covariates exhibiting grouping structures, such as a group of indicator variables representing a multilevel categorical covariate, or scientifically meaningful groups based on prior knowledge. A potential future research direction could involve redefining the groups within the proposed framework. This could be achieved by considering all outgoing connections from a group of input neurons as a single group, enabling group selection and accommodating the presence of grouping structures.

In conclusion, our study exhibits the advantages of employing group concave regularization for sparse-input neural networks. The findings highlight its efficacy in consistently selecting relevant variables and accurately modeling complex non-linear relationships between
covariates and outcomes, across both low and high-dimensional settings. The proposed approach holds the promising potential to enhance modeling strategies and find wide-ranging applications, particularly in diseases characterized by non-linear biomarkers, such as oncology and infectious diseases.

This research was supported in part by the National Institutes of Health Grants R01CA256157, R01CA249279, 1R21CA263950-01A1, the United States Army Medical Research Materiel Command grant Award Number HT9425-23-1-0393, and the Prostate Cancer Foundation Challenge Award.
