Sparse-Input Neural Network using Group Concave Regularization
Authors: Bin Luo
Susan Halabi
Venue: Accepted by TMLR ()

Abstract: Simultaneous feature selection and non-linear function estimation is challenging in modeling, especially in high-dimensional settings where the number of variables exceeds the available sample size. In this article, we investigate the problem of feature selection in neural networks. Although the group least absolute shrinkage and selection operator (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. We provide a rigorous theoretical analysis of the proposed framework, establishing finite-sample guarantees for both variable selection consistency and prediction accuracy. These results are supported by extensive simulation studies and real data applications, which demonstrate the finite-sample performance of the estimator in feature selection and prediction across continuous, binary, and time-to-event outcomes.

Published in Transactions on Machine Learning Research (11/2025)
Sparse-Input Neural Network using Group Concave Regular- ization
Bin Luo
bluo Whoneeenesa
School of Data Science and Analytics
Kennesaw State University
Marietta, GA 30060, USA
Susan Halabi
susan.halabi@dduke. edu
Department of Biostatistics and Bioinformatics
Duke University
Durham, NC 27708, USA
Reviewed on OpenReview: https: openreview. net/ net/ forum? id= LHUee
Abstract
Simultancous feature selection and non-linear function estimation is challenging in modeling, especially in high-dimensional settings where the number of variables exceeds the available sample size, In this article, we investigate the problem of feature selection in neural net- works. Although the group least absolute shrinkage and selection operator (LASSO) has been utilized to select variables for learning with neural networks, tends select unim- portant variables into the model to compensate for its over-shrinkage. To overcome this limitation, propose 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 l2 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 tackle the challenge of complex optimization landscapes. We provide rigorous thcorctical analysis of the proposed framework, establishing finite-sample guarantees for both variable selection consistency and prediction accuracy. These results are supported by extensive simulation studies and real data applications, which demonstrate the finite-sample performance of the estimator in feature selection and prediction across continuous, binary, and time-to-event outcomes.
1 Introduction
Over the past few decades, advancements in molecular, imaging, and other laboratory tests have led growing interest in high-dimensional data analysis (HDDA) (Donoho et al. 2000). This type of data involves large number of observed variables relative to the small sample size, which presents considerable challenge in building accurate and interpretable models. For example, in bioinformatics, hundreds of thousands of RNA expressions, genome-wide association study (GWAS) data, and genomic data are used to understand disease biology and the correlation with clinical outcomes, with only hundreds of patients involved (Visscher et al. 2012. Hertz et al.] 22166) Kim & Halabi 2016 Beltran al. 2017 Sumiyoshi et al., (024).. To address the curse of dimensionality, feature selection is critical step in HDDA modeling. By identilying the most relevant features that capture the relationship with clinical outcomes, feature selection enhances model interpretability and improves generalization.
"Corresponding author.

---

Published in Transactions on Machine Learning Research (11/2025)
There are various methods for feature selection, including filter methods (Koller & Sahami
1996: Guyon &
Flisseeff 2003 Gu et al. 2012 wrapper methods (Kohavi & John 1997 Tixaa et al. 2004 Tang et al. 2014),
and embedded methods (Tibshirani 1996 Zou 2006; Fan & Lij 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 methods can perform simultancous parameter estimation and feature selection by shrinking some of the parameter coefficients to exact zeros. While LASSO has been widely used obtain sparse estimates 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 Zoo 2006), the minimax concave penalty (MCP) (Zhang et al. 2010 and the smoothly clipped absolute devialion (SCAD) Fan & 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-parametrie extensions have been proposed to resolve this problem (Lin & 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. (014) 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 imaging (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 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 genoralization performance for neural networks (Liu et al., 2017).
& Krauus 2018). Their state-of-the-e
Recently, novel approaches have been developed that use 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 imput nodes (Liu et al. 2017 Scardapaae et al. 2017:
Feng & Simon 2017 These methods treat all outgoing connections from a single inpul. neuron as a group
and apply the LASSO penalty to the L2 norm of weight vectors of each group. Other LASSO-regularizee
neural networks in feature selection can be found in the work of Li et al. (2016) and Lemhadri et al. (0221
However, LesSO-regularizedd neural networks tend to over-shrink the weights of relevant variables, leading to the inclusion of many false positives. The adaptive LASSO was employed to alleviate this problem (Dinh & 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 4 regularization by introducing stochastic gates to the input layer of neural networks. They considered li-like regularization based on 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.
Recent thcoretical advancements have provided rigorous guarantecs for sparse deep learning. For instance, Sun al. (2022) proposed method justified under a Bayesian framework that can learn a sparse deep neu- ral network with theoretical guarantees for posterior consistency and variable selection consistency. Their approach first trains dense network with mixture Gaussian prior and then sparsifies its structure using Laplace approximation of marginal posterior inclusion probabilities. Building on a different theoretical perspective, Lederer (024) developed statistical guarantees for sparse deep learning from a high-dimensional statistics perspective, providing oracle inequalities that guarantee near-optimal prediction accuracy for var- ious types of network sparsity, including connection and node sparsity. Despite this progress, however, theory that provides direct frequentist guarantees for variable selection consistency of penalized estimators, cornerstone of classical high-dimensional statistics, remains important area for development in deep learning.
2

---

Published in Transactions on Machine Learning Research (11/2025)
In this paper, we propose novel framework for sparse-input neural networks using group concave regular- ization to overcome the limitations of existing feature selection methods. Although folded concave penalties like MCP and SCAD have been shown to perform well in both theoretical and numerical settings for feature selection and prodiction, 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 polential of concave penalties for feature selection in neural networks by providing a comprehensive approach for simultancous feature selection and function estimation in both low- and high-dimensional settings.
The key contributions of this paper are follows:
A unified framework for simultaneous feature selection and prediction: We introduce structured sparsity in neural networks by applying concave group penalties, treating all outgoing connections from single input neuron as group. An Iz-norm-based concave penalty shrinks entire groups of weights to zero, resulting in a parsimonious neural network that selects only a small subset of input variables.
An effective optimization strategy for stable solution paths: We employ backward path-wise opti- mization, dense-to-sparse approach that traces the solutions of the regularized neural network. This approach improves the stability of model selection and enhances the interpretability of the solution path across the full range of regularization parameters.
We provide rigorous theoretical analysis of our estimator, establishing non-asymplotic bounds on its prediction and estimation error and proving that it possesses the desirable oracle property under standard high-dimensional conditions.
Empirical validation across diverse data types: Through extensive simulations and real-data analysis, we demonstrate that our method outperforms existing approaches in feature selection consistency and prediction accuracy across continuous, binary, and time-to-event oulcomes.
The rest of this article is organized as follows. In Section 2, we formulate the problem of feature selection for generic non 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 establish the theoretical properties of our proposed estimator. Section 5, we conduct extensive simulation studies demonstrate the performance of the proposed method. In Section 6, apply the proposed method to real-world datasets. Lastly, in Section 7, we discuss the results and their implications. The thcoretical proofs, implementation details, and supplementary numerical results are provided in the Appendix.
2 Method
2.1
Problem setup
Let X E Rd be d-dimensional random vector and Y be response variable. We assume the conditional distribution Prix depends on a form of f(Xs) with a function f E F and a subset of variables S C {1,-.. ,d.. We are interested in identifying the true set for important variables and estimating function f so that we can predict Y based on selected variable Xs-
At the population level, we aim to minimize the loss
min Ex,y ((((s),
JEFSS
where é 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 obscrvations (X,, Y;) 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 and aim to approximate f nonparametrically
3

---

Published in Transactions on Machine Learning Research (11/2025)
using neural networks. Thus, our goal is to develop efficient method that simultancously select variable subset S and approximate the solution f for any given class of functions using sparse-input neural network.
2.2 Proposed framework
We consider function estimators based feedforward neural networks. Let Fn be class of feed forward neural networks fw I + IR with parameter w. The architecture of multi-layer perceptron (MLP) can be expressed as a composition of a series of functions
(w( = LD o ID-1 00 ................................ .... [ ER.
where o denotes function composition and (xx is an activation function defined for each component of x. Additionally,
[x(x) W, b,, 0, D....,,
where W, E Rd-1 is weight matrix, is the number of hidden layers, is the width defined the number of neurons of the i-th layer with do d, and b, R4 is the bias vector in the i-th linear transformation I. Note that the vector E IRP is the column-vector concatenation of all parameters in (W,,b; : =,11 0.1...,, D). We define the empirical loss of fw as
n
(w) Σ (x,().))
Ideally, the sparse-input neural network fw should rely only the important variables, meaning that Wa, for j & S, where Wa, denotes the jth column vector of Wo- In order to minimize the empirical loss L.,ww) (w) while inducing sparsity in Wo, propose to train the neural network by minimizing the following group regularized empirical loss
= arg min Cn(w) +
Σ Pa (I aall [a j=1
(1)
wER*
where || ||2 denotes the Euclidean norm of vector.
The objective function in Eq. comprises three components:
(1) Lu (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) PX is concave penalty function parameterized by À 2 0. To simultancously 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 ex is designed to shrink the weight vectors of specific groups to exact zeros, resulting in neural network that utilizes only small subset of the original variables.
(3)
The term a [|| is ridge regularization that plays critical role in ensuring both model stability and the effectiveness of the feature selection mechanism. Note that the group penalty Px acts only the input layer weights Wo. Without the ridge penalty (a = 0), neural network could cireumvent this penalty by shrinking Wa while inflating weights in subsequent layers to compensate, leaving the network output unchanged and rendering feature selection inellective. The ridge penalty prevents this compensatory behavior by regularizing the norms of all weights the network (w), encouraging smaller and well-balancod weights throughout the architecture. Morcover, in Scction 4, the ridge term is expressed as norm constraint on the full parameter vector, which essential for deriving statistical guarantees.

---

Published in Transactions on Machine Learning Research (11/2025)
It should be noted that when the number of hidden layers D 0, the function fw reduces to linear function, and the optimization problem in Eq. (1) becomes the framework of elastic net (Zou & Hastie 2005), SCAD-L2 iZeng & Xie, 2014), and Mnet (Huang et al. 2016), with the choice of ps to be LASSO, SCAD, and MCP, respectively.
2.3 Concave regularization
There are several commonly used penalty functions that encourage sparsily in the solution, such as LASSO Tibshirani 1996p, SCAD Faan & Li 2001,, and MCP (Zhan et al., 2010). When applied to the l-- norm of the coefficients associated with each group of variables, these penalty functions give rise to group regularization methods, inchuding group LASSO (GLASSO) (Yuan & Lin 2006), group SCAD (GSCAD) (uuo al. 2015), and group MCP (GMCP) (Huang et al. 2012). Specifically, LASSO, SCAD, and MCP are defined follows.
LASSO
Px(1) A/1|-
SCAD
A/t| for |t| S X,
Px(t)
3- 2
for |t| S a,, for ||| > ad,
2(a 1)
(a+1)X² 2
where a > fixod.
MCP
(1- (1-000A)A(((0000000000 Xa
2
px(t) sign(t)s
dz,
where a > fixed.
It has been demonnstrated, both theoretically and numerically, that the folded concave regularization mel hods of SCAD and MCP exhibit strong performance in terms of feature selection and prediction (Fan Li, 2001: Zhang 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 / >> aλ. Given the benefits of the concave penalization, we propose using the group concave regularization our framework for simultancous feature selection and function estimation.
3
Implementation
3.1 Composite gradient descent
The optimization in Eq. [ is not a convex optimization problem since both empirical loss function Ln (w)
and the penalty function PA can be non-convex. To obtain the stationary point, we use the composite gradient
descent algorithm [Neteerov] 2013). This algorithm also incorporated in Feng Simon (2017) Lemhadri
ct al. (2021) for sparso-input neural networks basod on the LASSO regularization. The local convergence of the composite gradient descent algorithm for nonconvex regularization was established in Gong al. (2013)
Denote (w) (w) + al|w| the smooth component of the objective function in Eq.
The
composition gradient itcration for epoch t is given by
w ++1 = arg min
- w
w
++
12
Px
(2)
2
w
=
5

---

Published in Transactions on Machine Learning Research (11/2025)
where wh+1 =w' -' is the gradient update only for the smooth component (n,a(w") that be computed using the standard back-propagation algorithm. Here > 0 is the step size for the update and can be set as fixed value or determined by employing the backtracking line search method, as described in Nesterov] (013).. Let Aj represent the index set Wo, within w. We define A as the index set that
includes all weights in the imput layer, given by A = -U Aj}. By solving Eq. (2) obtain the iteration form w A = - and
Ay = h 17, A), for j = 1........ .d.
Here, A refers to the complement of the set A, and the function h represents the thresholding operator, which can be determined by the specific penalty Px By taking Ps 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
(GIAsSO(2175)) = 7A(
GSCAD
Sg(2,
if ||=|l2 & (7 + 1)x, if (7 1) < ||=|| S ad, il |||||| ad.
hGSCAD
w-1- Sg(2,
a-1):
and
GMCP
a a-y Sg(,
if || S at, if |||||| > at,
NGMCP(2;7,N)
where Sg(2;A) is the group soft-thresholding operator defined as
Sy(z;A)= (1
A
z.
|
Therefore, we can efficiently implement the composite gradient descent by integrating an additional thresh- olding operation into the input layer. This operation follows the gradient descent step using the smooth component (w). The caleulation for epoch can be summarized as follows:
(1) Compute the gradient of the smooth component V VL (w") using Jack-propagation
(2) Perform the gradient update for the smooth component to get an intermediate estimate:
H+1 -- w' >VVLn
(3) Apply the thresholding operator to obtain the updated weights w(+1;
++1 w
1+1 do
w ++1 + h W 17,), for j = l....,d. d.
A
The final index set of the selected variables is $ = (j WA, + 0}. Note that we consider y as a scaling factor in the thresholding operator. When 7 = 1 in Step (3), the solutions for GLASSO, GSCAD, and GMCP align with the closed-form results established Wei & Zhu (2012).
3.2 Backward path-wise optimization
We are interested in learning neural networks not only for a specific value of X, but also for a range As where the networks vary by the number of included variables. Specifically, we consider range of from Amin; where the networks include all or an excessively large number of variables, up Amazi where all variables excluded and [W] becomes zero matrix. Since the objective function not convex and has

---

Published in Transactions on Machine Learning Research (11/2025)
multiple local minima, the solution of Eq. ① with random initialization may not vary continuously for À E [Amin: Amaz]: resulting in highly unstable path of solutions that are regularized by X.
To address this issue, we consider a path-wise optimization strategy by varying the regularization parameter along a path. The goal is to generate a stable solution path, which we define as a sequence of solutions where small changes in the regularization parameter result in correspondingly small changes in the model's 's weights and sparsity level, avoiding the large, erratic jumps often seen with random initializations for each (see Figure D. It is important to clarify that the stability of the path is product of this overall strategy, not property of the optimizer (e.g., Adam) used for any single value of X. In this approach, we use the solution of a particular value of À as warm start for the next problem. Regularized linear regression methods (Friedman et al., 2007 2010 Breheny & Huang 2011) typically adopt a forward path-wise optimization, starting from a null model with all variables excluded at Amaz and working forward with decreasing As. However, our numerical studies on sparse-input neural networks revealed that starting with a sparse solution as the initial model does not lead a smooth transition to a dense model. To address this challenge, we implement backward path-wise optimization approach, starting from dense model at. the minimum value of Amin and solving toward sparse models up to Amaz with all variables excluded from the network. This dense-to-sparse warm start approach is also employed in (Lemhadri 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 = f(X) + e, where ex X4) with informative and 16 nuisance variables, and each X, and e follow the standard normal distribution. We present more details of the simulations in Section 5, Figure [ shows the solution paths of GSCAD, GMCP, and GLASSO based different types of optimization. In these plots, each line represents the path of the group weight norm for single input variable; green lines indicate the four true informative variables and gray dashed lines indicate the sixteen muisance variables. It is observed that non-pathwise optimization (top-left) leads to fluctuations or variations in the solution path, whereas forward path-wise optimization tends to remain in the same sparse model (GMCP, top-middle) or experience fluctuation solutions (GLASSO, lop-right), until transitioning to a full model with a sufficiently small A. In contrast, backward path-wise optimization (bottom panels) yields smoother solution paths for informative and mussance variables. Importantly, GLASSO (bottom-right) tends to over-shrink the weight vectors of informative variables and include more variables in the model, while GSCAD (bottom-left) and GMCP (botton-middlee) are designed to prevent over-shrinkage and offer a smooth transition from the full to the null model.
[(X) = ( +X + X2
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 number of gradient descent iterations required for each X. Therefore, the bulk of the computational cost occurs at Amin: and a lower number of iterations for the remaining As 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 imput features, this approach can significantly speed up computation, particularly for higl-dimensiioall data.
3.3 Tuning Parameter Selection
Two tuning parameters are required in our proposed framework: the group penalty coefficient À and the ridge penalty coefficient a. The former controls the number of selected variables and yiclds sparser models for larger values of X, while the latter imposes a penalty on the size of the network weights to prevent overfitting.
7

---

Published in Transactions on Machine Learning Research (11/2025)
Non-pathwise Optimization,
Forward Pathwise Optimization,
Forward Pathwise Optimization, GLASSO
Informativn variale
Information variable
Nuturcee variable
1.5
1.55
11W
0.7
03-
0.30
-25
-15
-05
-lag50
-lag0)
Backward Pathwise Optimization,cm GSCAD
Backward Pathwise Optimization, GMCP
Backward Pathwise Optimization, GLASSO
Infermatve
Naluance variabee
Nalaance veraale
1.50
1.5
0.50
0.2
0.5
0.50
0.0
Figure 1: Solution path of l2 norm of the weight vector associated with each imput node |/oollle Fach line shows the path for single variable's group weight norm; green lines correspond to the four true informative variables and gray dashed lines to the sixteen nuisance variables. The bottom panels illustrate the "smoother" solution paths generated by our backward pathwise strategy, where weight norms change predictably without the erratic jumps seen in the non-pathwise (top left) and forward pathwise (top middle, right) approaches. Note that the x-axis, log(A), differs between plots because the absolute scale of the regularization parameter is not directly comparable across penalty types. The individual panel descriptions are as follows: Top left: Non-pathwise optimization using GMCP. All the neural network weights are initialized by drawing from NO0, 0.1) for each X. Top middle: forward path-wise optimization using GMCP. It starts from the null model and computes the solution with decrcasing X. Random initialization is used before the selection of the first set of variables. Top right:forward path-wise optimization using GLASSO. Bottom left: backward path-wise optimization using GSCAD. Bottom middle: backward path-wise optimization using GMCP. Bottom right: backward path-wisc optimization using GLASSO.

---

Published in Transactions on Machine Learning Research (11/2025)
In all numerical studies presented in this paper, we adopted 20% holdout validation set from the training data. The model was trained using the remaining data, and the optimal values for X and 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.ce wtt.:/ub..cmmm//i/////a..a
4 Theoretical Properties
In this section, we establish the theoretical underpinnings of our proposed estimator, providing non- asymptotic, finite-sample guarantees on its performance. Our primary goal is to show that, under standard high-dimensional regularity conditions, our method achieves strong statistical properties in terms of both prediction accuracy and variable selection consistency.
For the theoretical analysis, we constrain the norm of full vector w instead of using a ridge penalty. Our problem formulation for training the sparse-input neural network becomes:
I Σ
w E arg min WEB
(lw (X),,))) Σ px Poolla)
(4)
where θ = {w& E RP: /|w/22 S K} for a constant K > 0.
As common in the literature for non-convex problems, our analysis focuses on the properties of a suitable local minimizer, leaving the computationally challenging issue of finding a global minimum future work.
4.1 Preliminaries and Notation
Suppose covariates X have support X' C [-X Xaaz] and X contains some open set. We assume the true data generating process depends on function f that uses only subset of important variables S C (l. Let E denote the expectation with respect to the joint distribution of (X,Y). Given i.idd. observations, we denote the empirical expectation E,,
Let the set of optimal neural networks that minimize the expected loss be denoted
FQ = arg min E[e(/((X), Y)]. wee
(5)
Due to the symmetries inherent in neural network architectures (e.g., permutation of hidden nodes), the optimal parameter vector is not unique. Following the framework of Feng & Simon (207)), we assume the parameters in EQ' can be partitioned into a finite number of equivalence classes, denoted by Q 2 1. For any w' EQ", we assume that all weight groups tied the irrelevant features S are zero (i.e., W 0 for all j St). For any parameter vector w, let the closest element in FQ" be defined (()) E arg min ||w w2 We define the excess loss of a neural network with parameters w as
E(w) (X)(((()))))) w'eB = -E)((F(()) -aap(i)--((-aaa.) e(Sw-rom Y),
min [[L(S (X)
(6)
For our theoretical analysis, we partition the input layer's s weight matrix, Wo, column-wise according the true support set S, which contains the indices of the s = |S] important variables. This yields two sub- matrices: W Wo,s, with columns corresponding to variables in S, and Wa.se, with columns for the irrelevant variables in S We define Wupper the vector of all parameters not subject to the group-concave penalty; this includes the input layer's bias vector (bo) and all subsequent weights and biases, which has dimension dupper- This partitioning defines the oracle vector Ws, which combines all parameters of the true model structure (Wo,s and Wupper) for a total dimension d' = (s x di) + dupper: as well as the full vector w with total dimension p = + With slight abuse of notation, we denote the oracle parameter support set of ws by Sparam, where |Sparam| = d' This is defined the union of the index sets corresponding to W a. and Wapper:
9

---

Published in Transactions on Machine Learning Research (11/2025)
4.2 Statistical Guarantees for Estimation and Prediction
We now present our first main result, which provides deterministic upper bounds the estimation and pre- diction performance of any local minimizer whose irrelevant weights lie within the active penalty region. To establish these guarantees, we roly on a set of standard technical conditions widely used in high-dimensional statistics. These assumptions ensure that the loss function is well-behaved around the true parameters, that the model is identifiable, and that our nonconvex penalty can induce sparsity without introducing excessive bias.
Condition 4.1 (Local Strong Convexity). Let the parameter vector w be ordered such that the irrelevant weight groups Wo,s are first. There is a constant Amin > O, which may depend the network architecture, 8, and f*, but not on d, such that for all w EQ",
0 0 0 I
V ((((S(()) Y) w= 2 Mmain
(7)
This condition guarantees that the loss function has sufficient curvature around the true solution with respect to the relevant. parameters. This ensures well-defined minimum, enabling stable and consistent estimation. Condition 4.2 (Compatibility Condition). For all >> there is Xe >> 0 that does not depend such that
-w X S inf E(w) l|ws ws (w 2 and Il
< Σ jES
Woo
wo 122
W wupper
w
Wupper 112
weB
jES=
where ws and Wupper are as defined in Section
4.1
This is a standard requirement for high-dimensional analysis that links parameter error to prediction error. It ensures that a meaningful deviation from the true parameters results a quantifiable increase in the excess risk, which is crucial for model identifiability.
Condition 4.3 (Bounded Third Derivative). The third derivative of the expected loss function is bounded uniformly over θ by some constant G >> 0 that does not depend on d:
+ Jun,,,w,,B. Ow,,
sup max WEB JiJaJs
E(((fw(X)) SG.
(8)
This condition of bounded third derivative ensures the loss function is locally well-approximated by a quadratic. This regularity required for the Taylor series arguments that underpin our proofs.
Remark 4.1. Note that the bounded third-derivalive condition is not striclly satisfied by non-smooth acti- vations like ReLU (Naór & Hinton 2010), which, although computationally efficient, is non- con-differentiable at the origin. Consequently, our theoretical guarantees apply most directly the case of smooth activations, such as Softplus Glooro et al. (0011) or GeLU (Hendrycks & Gimpel 2016). In Section 5, we discuss the practical implemenlation using ReLU, where observe empirical perjormance consislent with the theorelical properties established for smooth activations.
Condition 4.4 (Propertics of the Nonconvex Penalty Function). The penalty is scalar function Px
o0) + [0, oo) satisfies:
(i) (x 0.
(ii) The function t + (t( is non-decreasing.
(iii) The function - ex(t)
non-incrcasing on (0, oo). This implics that px(t) is concave on [0, oo).
(iv) The function + (t) is differentiable on (0, oo).
(v) The derivative satisfies lim (1( A.
10

---

Published in Transactions on Machine Learning Research (11/2025)
(vi) There exists constant 8 > such that dA is the smallest value for which the derivative (1 for all 1 8X.
(vii) The derivative (1) is a concave function on the interval [0, 8A].
These properties formalize the requirements for penalty function that can produce sparse and nearly unbiased estimates. The key features are sharp penalty at the origin to enforce sparsity and vanishing penalty for large coefficients to reduce estimation bias. Common penalties that satisfy these properties include MCP and the SCAD.
Under these conditions, we are ready to state our main deterministic result, which provides explicit bounds on the estimation error, excess risk, and sparsity recovery for local minimizers of the penalized objective.
Theorem 4.1. For any Ā 0 and 2 let
IB - (e)(((w-ww))) (x), (x),
Txr = (1) WEB
sup
S TT}.
\ (I) Wupper
Wapper ||2 2=j=1 |Wo.) wob l12)
Suppose Conditions 4.1 4-4 hold. Let w be a local minimizer of equation 4 such that maxjes= Wooill 8X. Then over the set T5,T for any tuning porameter 2 4TX, such a local minimizer erists and satisfies
(i) Estimation Error: || w'sll S C3 TİJ
(i) Excess Risk: E(W) & (A+ (A+ s)C--
18) Sparsity Recovery: Ejes Wo, (A+TA)(L+0CC)
2-273
where
1 2K
C/
(9)
6O
Xeo
Amin
=
(10)
3 ( + vd
The proof of Theorem 4.1 is included in Section B.1
Note that the statement of Theorem 4.1 is entirely deterministic and the established bounds are expressed in terms of the tuning parameter À and the statistical complexity term 7. By having À vanish at the same rate of X, the estiionion error, excess risk and the norm of irrelevant weights converges the order of (o(As?df/2)) Op( 5522 5//2d C)) d2 22 respectively.
Remark 4.2. Theorem 4.1 establishes statistical consistency only for local minimizers whose irrelevant weights lie within the active penalty region, i.e., maxjes= [|a,ll S 6x. This requirement is not restrictive in praclice: by suilably choosing 5 or X, one can either enlarge the active penally region or drive the irrelevont weights toward zero, so that such a local minimizer may also coincide with a global minimizer. Moreover, it is preferable for the irrelevant weights to fall in this region, as this is the mechanism for inducing sparsity. In general, even if a global minimizer does not satisfy the condition, Theorem 4.1 guarantees the existence of statistically consistent local ménimizer that does.
The bounds in Theorem 4.1 depend on the good-set event Tx and an abstract statistical complexity term
A. The following theorem, adapted from Feng & Simon (2117), provides a high-probability bound on Tx for the common cases of regression and classihiction. This result quantifies the complexity, showing its dependence on the problem dimensions d, the number of obscrvations (n), and properties of the network.
Thcorem 4.2 (Probabilistic Guarantee for Event Tx.)). Consider the following two settings from Feng g
Simon 2017):
Regression setting: Let the loss E(...) be the squared-error loss, and the data be generated as y = (X) +G, where is sub gabsssasnnn random variable with meon zero, independent of X.
11

---

Published in Transactions on Machine Learning Research (11/2025)
2. Classification setting: Let the loss E( be the logistic loss, and the data be generated Y E (0,1) where P(Y = =
Suppose that the set of optimal neural networks EQ" is composed of Q equivalence classes. For each selling, there exists constant > O such that for
i co + dupper)
log Th n
(Vlog +
Xmaz log(ne2) (e((ne2)/oo(eee log((2d)
(11)
C1
we have for any T 2 1,
Prx,y (
where = 1/Vd,, 02 d1 K + Xmaz/c1 and M1 is constant that bounds the gradient of the loss with respect to the upper-layer parameters.
This result provides a concrete convergence rate for Theorem 4.1 The statistical complexity À scales roughly as Vlo(dd)/n, and setting the tuning parameter À /log d) n leads to the estimaaiion error, excess risk and the norm of irrelvant weights converges the order of Op((/og(d)//s$ d G), O, (p(log(d)/ns"// /d/"////////))))))))/)//))))")/ (og(()l)/)ss/////
respectively. These rates demonstrate that the estimator is statistically efficient in high-dimensional, nonparametric settings.
4.3 Oracle Properties
Finally, we verify that our estimator possesses the oracle property, gold standard for variable selection methods. An oracle estimator is an idcalized benchmark that knows the true support set in advance and is defined as the unpenalized minimizer of the empirical loss using only those true variables:
Σ (lw.(())))))
arg min
5
wsE
where Os denotes the subspace of θ restricted to the dimensions corresponding to the true model architec- ture, Og = {w E θ supp(w) C Sparam).
The following theorem shows that, under a minimum signal strength condition, our penalized estimator has stationary point identical to this ideal oracle estimator, when augmented with zeros for the irrelevant weights. Establishing this result demonstrates that our penalty is capable of correctly distinguishing true signals from noise, effectively mimicking the oracle's perfect knowledge.
Theorem 4.3 (Existence of Stationary Point with Oracle Properties). Suppose the conditions of Theorem hold. Assume further that the true nonzero weight groups satisfy the minimum signal strength condition
min
Wa SA + CO JES
where X denotes the statistical complerity of the oracle model, obtained by replacing the full parameter dimension d with the oracle dimension d' in the expression for X from Theorem 4.2
Then, on the event TI. if the regularization parameter satisfies A 2 Cs Csylog(d)/n for some sufficienlly large constant Cs, the augmented oracle estimator
(wS. 0)
is a stationary point of the penalized program equation O(exp(-C log d)) for some constant C 0.
This holds with probability least.
The proof of Theorem 4.3 is included in B.3 The oracle property provides a strong guarantee of variable
selection consistency for method. Specifically, it implies that. when signals sufficiently strong, the
12

---

Published in Transactions on Machine Learning Research (11/2025)
group concave penalty shrinks the coefficients of irrelevant variables exactly to zero, while leaving the coef- ficients of relevant variables essentially unbiased, as if they were unpenalized in an oracle estimator. This minimum signal condition is standard for such oracle results in high-dimensional theory; the required signal strength scales with the order of vs log(d)/n and thus vanishes as the sample size n grows, making the condition weaker for larger datasets. This result, offers a rigorous theoretical justification for the method's 'ss ability to accurately recover the true sparse model.
5 Simulation Studies
We assess the performance of the proposed regularized neural networks in feature selection and prediction through several simulation settings with various types of outcomes. In particular, 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 cstimator GMCPNet and GSCADNet with GLAS- SONet, neural network (NN) without feature selection (A = 0), random (survival) forest (RF), and the STG method proposed in Yamada al. (2020). These competitors were chosen to provide comprehensive eval- uation: GLASSONet serves the most direct alternative using standard convex penalty, NN highlights the challenge of overfitting without sparsity, STG represents a recent state-of-the-ate method with different sparsity mechanism, and RF provides a strong benchmark from the class of tree-based ensembles. We also include the oracle version of NN and RF (Oracle-NN and Oracle-RF) as benchmarks, which are trained with prior knowledge of the true relevant variables, serve as ideal benchmarks for the best achievable performance.
In our implementation, we replace the gradient update in Step (2) with Adam to improve computational efficiency and achieve faster convergence in practico. Specifically, set the scaling factor y = 1 in the thresholding operator and used a base learning rate of LR 0.001 for Adam. For a fair comparison across all neural network methods, we used a ReLU-activated multi-layer perceptron (MLP) with two hidden layers of 10 and 5 units, respectively. As discussed in Section 4.2 ReLU does not strictly satisfy the bounded third-derivative condition of Condition 4.3 since it is non-differentiable at the origin. We used ReLU in our simulations due to its computational efficiency and strong empirical performance. The results obtained with ReLU networks, presented in this section, closely align with the theoretical guarantees, as the non- differentiable points have measure zero for continuous input distributions, and ReLU can be viewed as the limit of smooth activations (e.g., Softplus) for which the condition formally holds.
A sensitivity analysis of the hyperparameter choices, such as network structure and learning rate, is provided in Section 5.3 and additional implementation details can be found in Appendix
5.1 Preliminary Simulation Study
In this section, we consider regression models of XOR-type and hicrarchical signal structures for continuous outcomes. We generale 500 i.i.d. random training samples according to the following models:
XOR-type Signals: Y = X1X2 +
Hierarchical Signals: Y =
In both models, X1 and X2 represent the first two coordinates of the covariate vector X E I2, each taking values in }+1} with equal probability, whilc e denotes standard normal error term. The remaining coor- dinates, X; for = 3, 3.....d, are uninformative variables. We conducted 20 simulations for varying d values ranging from 20 to 500. For each simulation, the performance of the trained model in both prediction and variable selection was evaluated on 500 independently generated random samples. For prediction accuracy, we report the test R² scores. For variable selection, we report the false positive rate (FPR) the percent- age of selected but unimportant covariates, defined FPR |Snga| [SaT
x 100%; and the false negative rate
(FNR) the percentage of important but non-selected covariates, defined as FNR s*ns| x 100%. Recall that S represents the true index sets of important variables and $ {j | ll2 * 0} denote the index
sets of selected variables. Simulation results are shown in Figure
13

---

Published in Transactions on Machine Learning Research (11/2025)
Our proposed methods, GMCPNet and GSCADNet, demonstrate clear advantages over other approaches, achieving relatively high prediction scores with low FPR and FNR across both models. Importantly, neural networks without feature selection tend to overfit as the number of noisy features increases, underscoring the importance of effective feature selection. Additionally, GLASSONet exhibits a tendency to overselect variables, resulting a high FPR due to the bias introduced by the LASSO penalty.
570
0.4
0.5
CMCMfet CSCADNe:
0.6
DMCCPett
0.2
PMR
KDP
0.7
0.7
GMCPlet
OSCADNe
-0.2
-0.2
100
200
400
500
100
400
100
300
400
1.0
STG
STG
GLASSONeT
0.04
ASSION
0.6
CSCAONet
0.02
GBCADNcc
0.4
RSF
S1G
0.3
0.2
GMCPNit
DSCADNet
-0.07
0.0
0.0
0.04
-0.
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
# of features
of festures
of features
Top row: simulation results for the model of XOR-type signals. Bottom row: simulation results for the model of hierarchical signals. The Rt2 scores, false positive rate (FPR), and false negative rate (FNR) are presented in the left, middle, and right columns, respectively. The central lines are the means, while the shaded areas represent standard deviations.
Figure 2:
5.2 High-Dimensional Simulations with Continuous, Binary, and Time-to-Event Outcomes
We evaluate our proposed methods under a more complex noulinear pattern across various outcome types, considering both low- and nigh-dimensionas scenarios. The data are generated through the following function:
f(x) log([ X1| 0.1) + X1X2 X₂ + exp(X3= + X4).
where each component of the covariate vector
(X,.... Xa)" E Rd 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 i.i.d. random samples with continuous outcomes, binary outcomes, and time-to-event outcomes in the following three examples, respectively.
Example 5.1. (Regression Model) The continuous response Y is generated from a standard regression model with additive error follows
y= Y=
where follows a standard normal distribution.
Example 5.2. (Classification Model) The binary response Y {0,1} is generated from Bernoulli distribution with the following conditional probability
P(Y 1(X) 1 exp (-f(X))
14

---

Published in Transactions on Machine Learning Research (11/2025)
Example 5.3. (Proportional Hazards Model) The survival time T follows the proportional hazards
model (PHM) with hazard function
h(t[X) ho(t) exp (f(X)),
(12)
where holl) is the baseline hozard function. Thus, T T= = ( (-log(U) exp [(x)), where U is a uniform random variable in [0, 1], and Ho is the baseline cumuulative hazard function defined Ho(t) - = ho (u)ddu.
We considered Weibull hazard function for Ho, with the scale parameler = 2 and the shape parameter 2. A proportion C of the samples is randomly selected to be censored. The censoring indicator is defined as 6; = for observed events and 8, = O for censored observations. The observed time for the ith individual is
= ==I(1 =) [[]](6,0))
where T, is the event time and C. is the censoring time. For censored individuals, C2 is drawn from an independent uniform distribution (0,T.), T:), ensuring that censoring precedes the event. In our simulation studies, we consider censoring proportions C 0, 0.2, and 0.4.
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 500.
We perform 200 simulations for each scenario. Similar to Section 5.1. the performance of the trained model is evaluated on independently generated n random samples. For prediction accuracy, we report the Ro score, classification accuracy, and C-index for the regression, classification, and survival models, respectively, In addition to FPR and FNR, we also report the model size (MS), which is the average number of selected covariates.
5.2.1
Results
Table D presents a summary of the feature selection performance of the four approaches: STG, NLASSONet, GMCPNet, and GSCADNet, across all simulation scenarios. We exclude the results of the STG method for Example 5.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 FPR and 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, Icading to larger model sizes in both LD and HD settings, which aligns with the inherent nature of the LASSO penalty.
Figure displays the distribution of testing prediction scores for the regression, classification, and PHM with censoring rate of C = 0.2. The complete results of the PHM presented in Appendix 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 scenarios. This strong empirical evidence, in which our proposed estimators perform on par with the Oracle-NN possessing prior knowledge of the true features, aligns with the oracle properties established in our theoretical analysis (Theorem 4.3) STG performs similarly to Oracle-NN the LDD setting of the regression model, but performance deteriorates in the HD setting and other models. Conversely, while GLASSONet outperforms or is comparable to the Oracle-RE 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.
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
15

---

Published in Transactions on Machine Learning Research (11/2025)
capabilities in selecting important variables with low FPR and low FNR, while achieving accurate predictions across various models. These methods show promise for tackling the challenges of feature selection and prediction in higl-dimensional data.
Feature selection results of STG, GLASSONet, GMCPNet, and GSCADNet under the regression, classification, and proportional hazards models. The false positive rate (FPR %), false negative rale (FNR %), and model size (MS) with standard deviation (SD) in hesentheses are displayed.
Table 1:
300, d 20
500, d 20
500, d 1000
Model
Method
FPR, FNR
MS (SD)
FPR, FNR
MS (SD)
FPR, FNR
MS (SD)
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)
Regression
GMCPNet.
2.2, 1.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,
(1.0)
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)
CLASSONel
89.2, 1.0
18.2 (2.8)
94.7, 0.2
19.1 (2.0)
16.3, 21.5
165.4 (92.9)
Classification
GMCPNel
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
(5.9)
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)
Survival (C =
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)
GLASSONet
98.0, 0.1
19.7 (0.8)
99.6,
19.9 (0.3)
16.8, 18.0
170.6 (49.0)
Survival 0.2)
CMCPNel
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)
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)
Survival
GMCPNet
5.8,
4.6 (1.5)
1.2, 0.1
4.2 (0.5)
0.0,
(1.0)
(C 0.4)
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)
5.3 Hyperparameter Sensitivity Analysis
We evaluate the robustness of the proposed method by conducting a sensitivity analysis under high- dimensional regression setting (d 1000, n 500), as described in Example The analysis examines the effect of three key hyperparameters: the thresholding scaling factor (v), the learning rate (LR) used in the Adam optimizer, and the network structure. Similar to previous examples, 200 simulations performed, and the average values of R²,, FPR, and FNR are reported for each configuration. Since GMCPNet exhibits similar performance to GSCADNet, we focus on reporting the results for GSCADNet. The findings are summarized in Figure ! which highlights the sensitivity to each hyperparameter individually while keeping the other parameters fixed. The fixed hyperparameter choices used in this study (2 = =1, LR = 0.001, and network structure [10, 5), i.c., two hidden layers with 10 and nodes, respectively) are markod on the plots for reference.
For 7 the results show that ? and 0.1 yield similar performance across all metrics, achieving relatively high R² while maintaining low FPR and FNR. In contrast, smaller values of (.001,0001) tend to under-select relevant features, as indicated by higher FNR and consequently lower I22 scores. For LR, values of 0.001, 0.01, and 0.1 exhibit stable performance with competitive R², low FPR, and low FNR. Smaller LR values (0.0001) lead slower convergence, resulting lower R² and higher FNR, while larger LR values (0.1) slightly improve R² and reduce FNR but also increase FPR. For network structures, more complex architectures such as [10, 10, and [20, 10, 5] provide modest improvements in R2 though they also
16

---

Published in Transactions on Machine Learning Research (11/2025)
n=300, d=20
n=500, d=20
n=500, d=1000
0.8
0.6
0.4
0.4
0.4
02
0.
0.2
0.0
00
0.0
-0.2
-0.2
-0.2
Cracle-RF
Cracle Cracle-Nig
GLASSONet
GMCtNet
CSCADNet
Cracle-PP Oracie MN
RF
GuAsSotet
GMCreet
CSCADNet
Gracle-RE Oracie dracie-Wlgr
GLASSOmlet
GMCPNet
GSCADNet
d=1000
0.85
0.85
0.85
0.75
0.75
0.75
Acuracy
0.60
0.60
0.60
0.50
Gracle-RF
Cracle NN AF
NN
STG
GLASSONet GMCPNet
GSCADNet
0.50
Cracle-FF
Cracie NN RF
NN
STG
GLASSNet GMCP
GSCADNet
0.50
Cracle-RF
Oracie NN
NN
STG
GLASSOMet
GMOPNet
GSCADNet
n=300, 20 C=
n=500, d=20, C=
n=500, 100000, C=
0.9
0.9
Cindex
Cindex
05
Dracle Oracle Oracle.A0
GLASSONt
Ghoptet
GSCADNet
Gracle Gracle-SFF
Orocle MNN
GLASSONet
GMOPNet
GSCDDNet
Cracle RSF
Gracle
Nii
GLASSONe Glinnet
GSCADNet
Figure 3:
R score of the proposed methods for the regression model outlined in Example Middle row: Accuracy of the proposed methods for the classification model outlined in Example 5.2
Top row:
5.1
Bottom row: C-Index of the proposed methods for the survival model outlined in Example 5.3 The dashed lines represent the median score the Oracle-NN. NN, used as a benchmark for comparison.
17

---

Published in Transactions on Machine Learning Research (11/2025)
result in increased computational costs and potentially higher risk of overfitting. In summary, although we fixed these parameters in our implementation with LR= 0.001, 1, and network structure (10,5)) our analysis reveals that selection accuracy and prediction performance remain stable for between 0.1 and 1, LR. between 0.001 and 0.1, and Network structures ranging from [10,5] to more complex architectures like 120,10,5. These lindings suggest that our method is robust to reasonable range of hyperparameter choices, demonstrating consistent performance across different configurations.
0.8
R²
FPR
0.6
FNR
Metrics
R²
0.4
FPR
FPR
FNR
FNR
0.2
0.0
10-3
10-2
10
10a
10-4
10-3
10-2
10-1
[5]
[10. 5] [10. 10, 5] [20, 10, 5]
Learning (Adam)
Network Structure
7, the scaling factor for the thresh- olding operator. Middle: Learning rate (LR) in the Adam optimizer. Right: Network structure. The network structure [4, 22.4..... 4x] represents the number of nodes in each hidden layer. The fixed choices used in our numerical study (7 1, LR = 0.001, and network structure [10,5]) are marked on the plots with an"x" symbol for each metric (R² FPR, and FNR).
Figure 4: Sensitivity analysis of GSCADNet hyperparameters: Left:
6 Real Data Example
In this section, apply our proposed framework to two real-world datasets to demonstrate its practical utility. The primary goal of these examples is to compare the behavior of different sparse neural network penalties and benchmark our method against strong, standard machine learning approaches, rather than to achieve absolute state-of-the-art performance on these specific tasks. To this end, the CALGB-90401 dataset is used to illustrate the framework's novel application to high-dimensional, time-to-event data, while the MNIST dataset intentionally subsetted to create a high-dimensional (p > n) problem that serves as a visualizable case study for feature selection.
6.1 Survival Analysis on CALGB-90401 dataset
We utilize the data from the CALGB-90401 study, double-blinded phase III clinical trial that compares docetaxel and prednisone with without bevacizumab in with metastatic castration-resistar prostate cancer (mCRPC) to illustrate the performance of our proposed method Kelly et al. (012). The CALGB- 90401 GWAS 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 our interest is studying the DNA damage repair genes, we only consider 625 SNPs
based on an updated literature search (Mateo et al. 2005 Wyatt et al., 2016 Mosquera et al., 2013 Robinson et al.) 2015 Abida et al.] 2019 De Lacre al.] 2017;; Sumiyoshi al.) 2024). We also include the eight clinical variables that have been identified 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 High, hemoglobin (HGB), PSA, and alkaline phosphatase (ALKPHOS). The final dataset contains d = 635 variables, n = 631 patients and a censoring rate C = 6.8%.
We consider the PHM the form of Eq. 12 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
18

---

Published in Transactions on Machine Learning Research (11/2025)
sets (n=526) and testing sets (n=105) using 5:1 allocation ratio. We apply the methods to each of the training sets and calculate the time-dependent area under the receiver operating characteristie curve (1AUC) 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 Our proposed method, Ge, out performs the others survival prediction (left panel). 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 tendency over-select variables and the proposed GMCPNet and GSCADNet select relatively smaller set of variables without compromising prediction performance (middle panel). The right panel of Figure demonstrates that GSCADNet successfully selects most of the key clinical variables and delects some of the important SNPs in predicting overall survival.
0.8
100
HGB
LDH High
0.7-
Veriables
80
ALKPHOS
SAUC
0.6
60-
ECOG
ALB
11011019
40
4964/0494
rs8075760
IS33149290
rs9204011
20
rs4150403
1s1894252
Is12329760
10V
RSF
GLASSONet
GMCPNet
GSCADNet
GLASSONet
GMCPNet
oSCADNet
1.0
Se lelection Proportions of OSCADNet
the number of selected variables for GLASSONet, GMCPNet, and GSCADNet. Right: Variables selected by GSCADNet with selection proportion 10% over 100 random splits.
Figure 5: Left: Boxplots of LAUC from testing set 100 random splits. Middle:
6.2 Classification on High-Dimensional MNIST
We aim to visualize variable selection in a nigh-ddmenssionl binary classification setting using the MNIST dataset. The MNIST dataset is well-known benchmark dataset in computer vision, consisting of grayscale images handwritten digits from to 9. In this study, we focus on distinguishing digits and 8, which share structural similarities that make the classification ion lask nontrivial. While other digit pairs may also exhibit visual similarity, this choice provides a meaningful evaluation of feature selection methods in identifying relevant pixels for classification. 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, resulting in 784 variables. To create a high-dimensional low-sample setting, we construct a training dataset by selecting 250 images of 7s and 8s each, yielding d 784 features and n 500 samples. Importantly, the class labels depend primarily on the central pixels, meaning an effective feature selection method should correctly identify and focus on these relevant regions. To ensure the feature space is not inherently sparse, introduce i.idd. standard Gaussian noise to the images. The trained models are evaluated 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 G. We observe that GLASSONet, GMCPNet, GSCADNet all achieve median accuracies greater than 91%, outperforming the other methods. While the heatmaps of feature selection show that GLASSONet, GMCPNet, GSCADNet consistently select relevant pixels high frequencies, GLASSONet tends over seleet variables and GMCPNet and GSCADNet choose irrelevant pixels in much lower frequencies (indicated by dark red colors).
19

---

Published in Transactions on Machine Learning Research (11/2025)
GLASSONet
0
0
5
5
7
10
15
15
20
20
20
25
25
GMCPNet
GSCADNet
0.91
7
Testing Accuracy
10
0.90
0.89
15
15
20
0.87
0.86
0.85
15
5
20
25
NN
RF
GLASSONet
GMOPNet
GSCADNet
Figure 6: Comparing feature selection and classification performance by STG, GLASSONet, GMCPNet, and GSCADNet. Top left: the image that takes the average of all images in the training set and shows relevant pixels in grayscale. Bottom right: testing accuracy for the classification 8s in a high-dimensional, low-sample MNIST setting with training dataset size = 784 and n = 500. Other panels: heatmaps depicting the selection frequencics of each pixel across 100 repetitions for cach method. Lighter colors indicate higher selection frequencies, with white highest and darker colors lowest.
78 and
20

---

Published in Transactions on Machine Learning Research (11/2025)
7 Discussion
In this paper, we have proposed novel framework that utilizes group concave regularization for feature se- lection 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 gradu- ally 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 opti- mization landscape across fine grid of tuning parameters, 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.
Among numerous feature selection methods, penalized regression has been extensively utilized. However, many of these methods rely on the assumption and application of linear theory, which may not capture the complex relationships ween covariates and the outcome of interest. In biomedical research, for in- stance, rescarchers often normalize data and employ penalized tochniques under a lincar 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-dimenssonal 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. treat- ment and care.
The runtime of our proposed method over solution path of As (with a fixed a) can be comparable to or even shorter than training single model with fixed X, such as the NN method without feature selection (A = 0). To illustrate this, examine the algorithm complexity of the NN method, which can be approximated 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 As has a complexity of O(udT'm), where represents the averaged number of inputs along the solution path with dimension pruning, and T' is the number of epochs for each À in the path. In our simulation with the HD scenario (d = 1000), we set = 5000, T" = 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 d d/2 = 500. Thus, ndT ndT'mm indicates that solving for an entire path of our proposed method requires similar computation as training single model. In real applications, especially in high-dimensional scenarios, the dimensionality usually drops quickly along the solution path. Therefore, 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' to be small for the first parameter Amin as well in the HD setting, to avoid overfitting of an initial dense model.
A central contribution of this work is the rigorous theoretical foundation we establish for the proposed framework. Our analysis provides non-asymptotic, finite-sample guarantees for any suitable local minimizer. In particular, we derive explicit bounds on both estimation error and prediction error (excess risk), showing that the method achieves statistical consistency even in high-dimensional settings. Furthermore, we prove a key oracle property: under regular conditions, the ideal oracle estimator, which is the solution attainable if the true relevant variables were known, is a stationary point of our penalized objective function. Although the non-convexity of neural networks precludes guarantees of convergence to a global minimum, our results demonstrate that the framework is properly formulated to yicld statistically reliable estimates and recover the true sparse model. This is further supported by our empirical studies, which indicate that the algorithm tends to converge to solutions consistent with the predicted theoretical performance and often comparable to the Oracle-NN benchmark.
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 load to an excoedingly complex optimization landscape, making it computationally infcasible. A potential solution is to employ a pre-screening step to reduce dimensionality before applying the proposed
approach (Fan & Lv 2010
21

---

Published in Transactions on Machine Learning Research (11/2025)
Another limitation of the proposed group-regularized method is its focus individual feature selection. This limitation becomes partieularly 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 rescarch dircction 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 single group, enabling group selection and accommodating the presence of grouping structures.
Furthermore, study is limited by its focus standard MLP architecture. This choice intentional, as it provides a clear setting to illustrate the properties of the proposed feature selection framework without the confounding influence of complex architectural priors. We recognize, however, that this choice may not yield state-of sfate-of-the-art predictive performance on tasks where specialized architectures, such as convolutional neural networks for image data, are more effective. A natural direction for future rescarch is to integrate our group concave regularization approach with such architectures, potentially achieving models that combine strong predictive accuracy with the feature-level interprelability offered by our framework.
In conclusion, study demonstrates the advantages of employing group concave regularization for sparse- input neural networks. Our theoretical results establish finite-sample guarantees, demonstrating statistical consistency and oracle-like performance of the method under regular conditions. The findings highlight its effectiveness in consistently selecting relevant variables and accurately modeling complex non-linear rela- tionships between covariates and outcomes across both low- and high-dimensional settings. The proposed approach holds promising potential to enhance modeling strategies and find wide-ranging applications, par- ticularly in diseases characterized by non-linear biomarkers, such oncology and infectious discases.
Acknowledgments
This research was supported in part by the National Institutes of Health grants RO1CA256157, RO1CA249279, R21CA263950, the US Food and Drug Administration grant 1U01DD007557001, the United States Army Medical Research Materiel Command grants HIT9425-23-1-0333 and 2119425-2511-0623, and the Prostate Cancer Foundation Challenge Award.
A Empirical Loss Function
The empirical loss functions Ln (w) for regression, classification, and survival models in Examples 5.1 5.3 are defined as follows:
Mean squared error loss for regression tasks, This loss function measures the average squared difference between the truc values Y, and the predictions (X,)).
L, = Σ fx()))
Cross-entropy loss for classification tasks. It widely used classification problems and quantifies the dissimilarity between the true labels Y; and the predicted proba bibitiiss Y; of class 1. The predicted probability Y; is obtained by applying the sigmoid function to fw (X.):
(w)
1 Σ n
Y log(Y) + (1 Yú) log(1 9:)
Negative log partial likelihood for proportional hazards models. It is derived from survival analysis and aims to maximize the likelihood of observing events while considering censoring information. It incorporates the event indicator 8,, which is 1 if the event interest occurs at time Y, and 0 if the observation right-censorod. The negative log partial likelihood defined as:
(w(w) - Σ
8, 8.fw(Xi) - 6, log Σ exp (f (X,)))
1=1
jen:
22

---

Published in Transactions on Machine Learning Research (11/2025)
Here, R, {i Y, 2 Y represents the risk set just before time Y. The negative log partial likelihood is specifically used in the proportional hazards model.
B Proofs
B.1 Proofs of Theorem 4.1
Lemma B.1 (Quadratic lower bound). Suppose Conditions 4.1
4.2
and
4.3
hold. For some constant
R > 0, suppose w E ④ satisfies
|[w w
w 1/2 S R
and
(13)
Σ jES
WWolll = ||aapper- willed 1/2 + Σ
(14)
jES
Then we have E(w) 2 ||ws w where
R²
C
(15)
60
Xro
Amin
en +V Vd
(16)
Proof. We begin by defining the la distance between the relevant parameters of w and the closest optimal parameters w-(w) as D(w)= = =|lws- Ws 1/1 The objective is to establish quadratic lower bound for the excess loss E(w) in terms of D(w).
Since the gradient of the expected loss vanishes at the optimum, a third-order Taylor expansion of the excess loss around (+)) yields:
E(w)= = 5( -
ww- (w (wwww w w(w)) + Tw
(17)
where Tw is the Lagrange remainder term
From Condition 4.1
the quadratic term is lower- bounded by:
(w *((w))
TI (e( Y)) (w (w w((w) 2 Amin (ww
(18)
From Condition 4.3 the remainder term is bounded as:
G
[rw] ||w (w)
(19)
The subsequent step to bound the 4-norm ||w - w= in terms of D(w). We decompose the norm into components corresponding to the irrelevant weights (S*) and the relevant parameters (ws):
||w w
Σ
||allhh + ||ws w w>(w)])
(20)
JES
noting that W for all E S The two terms of this decomposition are bounded separately.
The first term, corresponding to the irrelevant weights, is bounded using the Cauchy-Schwarz inequality, W s ddlllolllee which leads
Σ JESr
Wo,lhh S ddi Σ W
JeS-
23

---

Published in Transactions on Machine Learning Research (11/2025)
Under 14
assumed in the lemma, the sum of the norms is bounded, yielding:
Σ jEsc
Wosll - dd
||Wupper
w Mapper 1|2 +
Σ
W 0.j
Wol 112
jES
S dg
(|WWuper w Wupperll
112 + 555s| /
Wo.s l12
$ Vd((++533)(()
The second term, for the relevant parameters, is bounded similarly. Let d' be the dimension of ws- By Cauchy-Schwarz:
Ws (w) 11 S Vd - ww 112 w
D(w
Ws
Summing the bounds for these two components provides an upper bound for the total 4 norm:
- "( S 5(5)))) + Vd D(w)
-Vai(1++55 + v D(w)
Substituting this into the remainder incquality equation gives:
3 |nwl S (VdiLL+5V)) + vdt) D²(w)
(21)
By combining the lower bound on the quadratic term 18 with the upper bound for the magnitude of the remainder 21 we establish the following lower bound for the excess risk:
((w) 2 Amin D'(w) - (VT + v D2(w)
Now we use Condition 4.2 and apply Auxiliary Lemma Städler et al. 2010). as it is of the form h(2) 2 D(w). Therefore, applying the lemma yields the desired result:
A Cz3 with
E(w) (www//C3
where
1
R²
C
and
60
Xro
hman
6
(1+ + d3
Lemma B.2 (Lower Bound for a Class of Nonconvex Penaltics). Let the penalty function Px(t) satisfy
Condition 4.4 Then for any & [0, 8A], the penalty satisfies the lower bound:
px(1) 2 S X
Proof. By Condition 4.4- 4.4(iii), the function g(t) = Px(t)tt is non-increasing. Thus, for any E (0, dX), (()) 2 g(x),, which implies (x(11 2 px(2))/8X. The proof reduces to showing Px (x(5))/5X 2 X/2.
By Conditions 4.4 (v), (vi), and (vii), the derivative (*( concave function on [0,6X] with (x( X
and (,(88) O. A concave function's 's graph lies above its secant line. The secant line connecting (0, and ()),0) is T(u)= X(1 u/82). Therefore, p'x( (u) L(u) for E [0,5X].
24

---

Published in Transactions on Machine Learning Research (11/2025)
Integrating both sides from 0 to od gives:
Px(OX) = 1ex pu( > {
Ad L(u) du _oo
Dividing by aA yields px(oA)/XX 2 X/2, which completes the proof.
Now we are ready to prove Theorem 4.1 We first suppose the existence of local minimizer w that satisfies maxjes= WWo,ll I 5λ; we will establish this fact at the end the proof. Then the proof procoeds in three main sleps. First, we establish basic inequality derived from the optimality of w. Second, we use this inequality to verify that the compatibility condition holds. Finally, by combining the basic inequality with the quadratic lower bound from Lemma B.1 we derive the final bounds estimation error, excess risk, and sparsity.
For simplicity, we denote w w'.
Step 1: Basic Inequality. By definition of w as a minimizer and using Condition 4.4 (1) (px = 0):
P (1., v1)
Σ Px(lIWo.llee S n w (0r, Mi) + Σ
Px Cx
n
JES
Rearranging and using the definition of excess loss gives:
E(W) + Σ pa Wo.sll))) s (E,, - [E)(w (w)+ + Σ px (I W1yll)).
JES
Over the event TT. A(w, w*) Wupper w Wüpper lla + Σラ Wo.5 Wü,lle. The proof now proceeds by considering two cases based on the magnitude of A(w, w*).
Case 1: A(w,w") X.
In this scenario, the total estimation error A A(w, w") is already bounded by the statistical complexity term Å. The conclusions of the theorem are thus met, as the claimed upper bounds are of order À or X (or higher), which will be larger than the assumed error for any non-trivial choice of regularization.
To be more explicit, the individual bounds are satisfied as follows:
(i) Estimation Error: The quantity to be bounded is Iws - w3ll2- fiss-wslle
By definition, || - w I hiss-w'll
A(w, w") S. The theorem's bound of + + is of order or Å, and thus the inequality holds.
(ii) Excess Risk: We start from the basic inequality derived from the optimality of w:
Σ PA(WWollll)) S (Ee, E)(Ew- - -ww) + E(sx(IW5, px (lollll)) JES
E(W) +
JES
We bound the terms on the right-hand side. First, the stochastic error term is bounded by the definition of the event Tx
(E, - R) (w) S TT(AVA)
In this case, since < A, this term is bounded by T Next, because the penalty function Ps is Lipschitz with constant (from Conditions 4.4 (iii, vii, the second term bounded by:
Σ (ps((Wollle pa( pa(lWo.lla)) JES
W Woyllz AA.
jES
25

---

Published in Transactions on Machine Learning Research (11/2025)
Again, since A S X, this term is bounded by AX.
Combining these bounds and noting that the penalty on the inactive set 5 non-negative, we can drop it to get:
T
(TX AjA
This derived upper bound is order X² or AX. This is consistent with and smaller than the final bound of (A + (1 + x)CB stated the theorem, which derived from the more critical second case. Thus, the conclusion for the excess risk holds.
(iii) Sparsity Recovery: The quantity to be bounded is Ejes: [ Wo By definition, this is also less than or equal to A(w, w*) S X. The theorem's bound of order À X, so this condition is also satisfied.
Case 2: A A(w, w
We now focus the more challenging case where the estimation error is large. We start again from the basic inequality:
((W) Σ JESc
(..(lWollla) S TTAA w*) [(px(IlWa.la))) - (o)lWWo)))))).
jES
By Lemma B.2 since maxjes- [o S 8A, we have px Wo.lle) 2 / for j S Additionally,
from Conditions 4.4 (ili,)), the penalty Lipschitz with constant X, i.c., (px(t) - px(t2)) s A/1 - -2). This
gives:
E(w)
> Σ
Wo,ill s
TÀ
|I w Wupperlle +
Σ
|Wo,ll
+ +A+
Σ
Wa, - Wöylle
2
JES
JES
jES
Grouping terms yields:
1 1
E(w)
Σ
||Wollla ST upper w Wupperlle 12p (A TT))
Σ
|| o.5
Wo.slla
(22)
2
JES
jES
Step 2: Verifying the Compatibility Condition. We choose À such that A 2 4TX. This ensures
1t - - T}= TA O. From equation 22
Σ JES
Wo,lla
S
W Wupper
upper ll + 1 - TÀ
Wa,
Wo.sll2
TÀ
2
jeS
The first coefficient is < 1 and the second < 5. Thus, Condition 13 is satisfied. In addition, by assumption we have
Wws w3ll2 || + || | ll2 S 2K.
Therefore, the conditions of Lemma B.1 satisfied. Applying the lemma yields
E(w) > [w w
Step 3: Deriving the Final Bound and Proving Existence. Let D(w) = ||ws - w.sl2. From
equation 22 and Cauchy-Schware
E(W) + (
Σ
Woollls [
JESc
26

---

Published in Transactions on Machine Learning Research (11/2025)
Applying Young's inequality the right-hand side gives 1 3(AA+ (X) s) lower bound for E(w):
Diw) 20%
Substituting the
E(w) 1 T
Σ
WWolll - ((++T))3( 5) + [(())
jESr
Multiplying by 2 and rearranging yields the main result:
E(w) + (A- 2TX)
Σ
[Wo, s (A) s)
JES
This holds for any local minimizer satisfying our initial assumption. To show such minimizer exists, we must confirm that the parameter of the penalty function can be chosen to be consistent with this result. From the derived bound, we have:
(++1))(1+
max
Wo &
JES
2Tİ
We require this upper bound to be less than or equal to σλ:
(A+ (1) 6A.
2T
Let cTÀ for a constant e 2 4. After substitution and canceling TÝ on both sides, the condition becomes:
TA(c + () + s)CB S ScT = clc ccc22 2)
S.
c-2 2
Since the left-hand side is a finite constant that depends on the problem parameters (but not on n or d). we can always choose a penalty function with a parameter ó large enough to satisfy this condition. This confirms that local minimizer w satisfying our initial assumption exists for a suitable choice penalty.
B.2 Proof of Theorem 4.2
The proof of this theorem is direct adaptation of the proof of Theorem 2 in Feng & Simon 2017 It relies empirical process theory to bound the stochastic error term
The core logic follows their use entropy bounds for neural network function classes (Lemma 2 in their work) and a peeling argument (Lemma 4 their work). The primary modification for our setting involves adapting their bounds from the sparse group LASSO penalty the group LASSO penalty used here. This corresponds to setting their parameter a in the relevant constants (cı and 62), which leads to the formula for À presented in the theorem statement. We refer the interested reader to the original paper for the complete technical details.
B.3 Proof of Theorem
This section provides complete proof of Theorem stationary point of the penalized problem. The proof follows the Primaa-Dual-Witness (PDW) method, adapted for the non-convex setting of sparse neural networks.
43 which establishes that the oracle estimator is a
B.4 Foundational Lemmas
First, we establish a crucial lemma that provides a tighter estimation error bound for the oracle estimator and shows that it lies in the interior of the feasible set O with high probability. This allows us to use the simpler zero-gradient condition for optimality.
27

---

Published in Transactions on Machine Learning Research (11/2025)
Lemma B.3 (Oracle and Restricted Estimator Properties and Equivalence). Suppose the conditions of Theorem 4.1 are met. Let an onacle estimalor w3 and a restricted penalized estimator ws be defined any solution to their respective programs:
n
W S
arg min w
(fws (X),,
1
n E(fws (X.), + E
Wg E arg min w
n
jES
where Os is the parameter space restricted to the active set S. Then, with probability at least 1 O(L/n):
(a) Estimation Error: The estimators satisfy the following bounds:
(1) [| w'sll2 S CITi-V++S - - & C((A+XX)/++
(8)
(b) Equivalence of Solutions: Suppose in addition that the minimum signal strength condition in Theorem 4.3 holds. Then the set of global minimizers for the oracle program and the restricted penalized progrom are identical.
Proof. The propertics hold on the high-probability event Tx-
Part (a): Estimation Error Bounds
(i) Oracle Estimator Bound: By definition of wg. we have the basic inequality for its excess risk E(WO) S EE,, R)(ww- (wo). On the event T this is bounded by:
TX*
w Wapper
+
> Elloo jES
0.5 Wo.slla
We apply the Cauchy-Schwand inequality to the sum of norms. Letting the vector of norms be v (l) Il....,,Il I1z) and vector of ones be 1, we have v 1 |||||||1||2_ This gives:
I upper
w Wapper lla
Σ Elle, Wö S vI+ JES
/
//we upper Wüpperl1
11 +
Σ JES
we
W
0.j
The term under the square root is preciscly ||w] w's Letting D(w")= = l/ w'sl22, we have E(wO) s TX* &DWO). Applying Young's inequality and the quadratic lower bound E(WO) 2 D2 /CB yields the stated bound.
The restricted penalized program is instance of the general frame- work in (1) with reduced feature set. Therefore, the conclusion of Theorem 1.2(i) applies directly, yielding the stated bound.
(ii) Restricted Estimator Bound:
Part (b): Equivalence of Solutions
We prove that the solution sets are identical by showing that the solutions must have equal loss and penalty values on the high-probability event.
Step 1: Establish Equality of the Penalty Term
We show that for any solutions w0 and Ws, their total penalty values are identical. This requires proving that all their active weight norms exceed the threshold 5A.
28

---

Published in Transactions on Machine Learning Research (11/2025)
For the restricted estimator wg: By the minimum signal strength condition in Theorem 4.3 for each j E S.
Wo,llz 1|2 > 8λ + + T
From Part (i), the restricted estimator satisfies
W g w3ll2 C
Since the error of a single group is bounded by the total error,
| Wo C TX*
Substituting into the signal condition gives
| Ad + Wo. Wo,ll2
By the reverse triangle incquality,
[oolle > Wo - [[o,, Wo > 8X.
For the oracle estimator W From Part (a)(i, the oracle estimator satisfies
w w S C VI+s,
which is strictly smaller than the bound for Ws since À O. Hence, for cach j E S,
[a ll2 > 8A + C3 dd [ Wö.ll2.
Applying the reverse triangle inequality again yields
[Wo > 82.
Since all active component norms for both estimators exceed oλ, they lie in the region where the penalty function Px(t) is constant. This proves the equality of their total penalty values:
E (2) = E (l)llo))ll)).
JES
JES
Step 2:
Establish Equality of the Loss Term and Conclude Equivalence
By the definition of ws as a global minimizer of the penalized objective, its value must be less than or equal that of any other point, including wg:
(ws)) + Σ CP S (w + Σ (
jES
jES
Since we have just shown the penalty terms are equal, leaving:
L. S L., (w
Furthermore, by the definition of the oracle estimator w as global minimizer of the loss, know that Ln(Wg) < (ws). The only way both inequalities can hold is if the loss values are equal:
We have now shown that any restricted solution W achieves the global minimum value of the oracle loss. Therefore, any restricted solution is also an oracle solution. Since their loss and penalty values are identical, their penalized objective values are also identical. As ws is a global minimizer of the penalized program, we must also be one. Therefore, any oracle solution is also a restricted solution.
Since both directions of inclusion hold, the sets of global minimizers for the programs identical. 7
29

---

Published in Transactions on Machine Learning Research (11/2025)
Next establish two crucial lemmas about the properties of the gradient components.
Lemma B.4 (Boundedness and Lipschitz continuity of the gradient). Under Condition 4.3 (bounded third derivalève) and with the compact parameter space θ {wr R |wl]2 I K), the following hold for ony fized (x, and any indices j,k:
(a) Boundedness: The function Sw. M) = Woo. ([w (2), is bounded on O That is, there erists B > 0 such that
I9w.j, (2, v) SB for all w E.
(b) Lipschitz continuity: The map w (2, v) is Lipschitz on θ. There exists Ly > 0 such that for all w, w' E
19w.sdll,, y) 9w Iw' (x, 1)1 s Lg w w'll.
Proof. Preliminaries. Let.
G(w) == (w(((f(z).)) E RP and H(w) :- V E R
Condition 4.3 states that the third derivatives of to fw with respect to uniformly bounded on θ In particular, the Hessian H(w) is contimuuusly differentiable on θ and hence continuous. Moreover, bounded third derivatives imply that H-) is Lipschitz on θ; i.e. there exists M < xo such that for all w. O
[[(w) H(v)llop M//w -aww-vlll v
where || llop denotes the operator norm. The compactness of then ensures sup [H())]lop < co.
Part (a) Boundedness. The coordinate function Iw.j y) a component of the vector-valued function G(w)) From the preliminaries, G(w) is continnous on θ. Since θ is compact, the Extreme Value Theorem implies that each coordinate of G is bounded O. Therefore there exists B 0 such that 19w.j.& (z, y)1 < B for all w E O, proving part (a).
Part (b) Lipschitz continuity. Fix w, w' E O. Using the fundamental theorem calculus (integral form) for vector-valued functions,
G(w) G(w') [[ H(w + +(w -w)) (w w') dt.
Taking the (j, coorddiatte and absolute values gives
19w.5.6(3,,)) v) Iw (. v)1 s 1.5/ H(w + +(w w')) dt 1|w w'lz.
Because H-) continuous on the compact set θ, its operator norm attains finite supremum θ. Define
Lg == sup [H])]lop < oc. wEE
Then
1| 1 H(w + ((w-w)) de < L [[ + 1(w w dt S Lg
and hence
19w.j.k (2. Iw (., y)1 < L wwwwwlla
This proves that Iw.g.u (z, 8) is Lipschitz on with constant Lg, completing the proof of part (b) and of the lemma.
30

---

Published in Transactions on Machine Learning Research (11/2025)
Lemma B.5 (Gradient concentration over fixed support). Let Y.) be i.i.d. observations and the function class Fs be the set of scalar components of the gradients, restricted to the oracle porameter subspace. Os = (w E θ supp(w) C Sparam):
Fs {9..
WV. (z)(w)) v)lx w E Os j = =.......d d, k 1= d1]
Assume there exist B, and a compact set TRP in a ball of radius such that:
1. (Uniform boundedness) for all g Fs, Io(a, v)1 s B
(Uniform Lipschitz) for all j, ls, (2. y), I[V (( (2), y)4a WW.,, V)lel -|[ w/l2.
(Uniform variance bound) supgeTs Varlg(X, Y) S o².
Then there erist absolute constants C₁ C2 C3 such that for
d' + log(d d1)
Acc
n
the following holds:
PP
sup [(P, gEFs
E)gl s C1 BAco) 2 1- C2Cexxp (-
Proof. The proof procoods three steps.
Step (Covering). We cover the set of functions Fs. These functions are generated by parameter vectors w from the set Os. While the vectors in Os are elements of R2, they are constrained to a fixed d' df-dimensional subspace defined by the support Sparam- The geometric complexity of this set is therefore equivalent to that of d d"=diieeniionll ball of radius K. The covering number for this set is given by:
3K
N(e, Os, Os, (12) 5
By the Lipschitz mapping property, the metric entropy (E, (2)) is bounded by:
(() < log(d d1) + log (E/ Os, ll (2)
3KL
s log(d di) + dr log
Step 2 (Symmetrization and Dudley's Bound). We use Dudley's entropy integral to bound the Rademacher complexity:
C J
Vn. /" (Hs du.
ER, (Fs)
Using our entropy bound, the integral is bounded by:
+ log(d di) + d' log(3KL d log(d dl).
This yields the expectation bound:
B
E sup (IP, E]gl
S
d' + log(d d1) BAGc.
Vn
9
Step 3 (Concentration to High Probability). We use Bousquet's incquality, which states that for any
> O, with probability at least -t,
202t 2
sup (P, g
E)gl < E,
n
3n
31

---

Published in Transactions on Machine Learning Research (11/2025)
where En is the expectation bound from Step 2. We choose terms on the right-hand side become:
enAGe for suitable constant > 0. The
En E BAGe,
V S SoveAoe,
20
2Bt 3n
X BcX
For n large enough, the Aoo term is dominated by the Aee terms. Thus, there exist constants C1 C2 C3 > 0
such that:
P sup (TP, - E)gl I C B A00) 2 1 C2exp(-Cyndac)) gEFs
which completes the proof.
B.5 Main Proof
Proof of Theorem 4-3 We verify that the augmented oracle estimator w (W3,0) 0) is stationary point
of the penalized program equation by checking the probability.
(KKT) conditions with high
(a) Primal feasibility: w θ is satisfied by construction.
(b) Stationarity the active set (S): First, we establish that w is an interior point of its program. Assuming the true parameters not on the boundary of the feasible set (i.e., l|wšlla S K/2), the interior point property follows from the triangle inequality and the estimation error bound in Lemma
B.3 (i):
K + +TX VI+
1 ll2 S Iw's 12 l/ w'll2
2
For sufficiently large sample size n, the second term is bounded by K/2, ensuring that [w < K. As an interior-point solution, its umpenalized loss gradient must be zero: [VC,(w))) 0.
Furthermore, Lemma B.3 b) establishes the equivalence of the oracle and the restricted penalized programs. This means w is also global minimizer of the restricted penalized program and must satisfy its KKT conditions. The stationarity condition for any j E S is Vwa, (w + 2j = 0, where Zj is a vector in the subgradient the penalty. Since the loss gradient is zero, the penalty subgradient zj must also be zero. Therefore, the stationarity condition for the full penalized program equation is satisfied for all weights in the active set S.
(c) Strict dual feasibility on the inactive set (S"): For any j S we must show that V Wo. Fn (w0) A. We decompose the gradient into its stochastic and deterministic parts:
V., (VL, V((W))) + (VL(W)) VL((w))
(T) Stochastic Error
Bias
Bounding (I): Stochastic Error. The conditions required to apply Lemma B.5
uniform bound-
edness, uniform Lipschitz continuity, and uniform variance bound- follow directly from Lemma Consequently, Lemma B.5 applies and guarantees that, with high probability, the stochastic error is uniformly controlled:
B.4
max JES
VWa,
& di sup (PP. - E]gl S Vd geFs
for some constant Ca, where the statistical rate Ace of order O( ((/log(d)/n).
Bounding (II): Bias. The population gradient is Lipschitz continuous (Lemma B.4 [VL[(O) vL(w*)ll < L w Since both wo and w* are zero the inactive set S", this distance simplifies to l w3ll2. Using the sharp estimation error bound for the oracle estimator from Lemma B.3 (i), we can bound the bias:
VV((>) VL(w")ll & L. COT''
32

---

Published in Transactions on Machine Learning Research (11/2025)
Crucially, the rate of the restricted complexity term X depends the oracle dimension d' (which scales with 8), not the ambient feature dimension d.
Combining the Bounds and Concluding the Proof. By combining these results, the total gradient norm for any inactive feature j E S is bounded by the sum of the two parts:
Cn Cn <
Vaicadoc + +CCTAVVI+
Stochastic
Bias
In the high-dimensional setting where d > the stochastic error term dominates because its rate Age scales with vlog d, while the bias term's 's rate X scales with the much smaller oracle dimension vlog d The entire gradient norm is therefore bounded by an expression of order (,((/og(d)/m) As established in Lemma B.5 the gradient concentration holds on an event with probability at least (cxp([CCloo log d)) for some constant C > Thus, by selecting a tuning paramcter À that satisfies Cs (lolog())/ for a large enough constant Cs, ensure that with high probability, the gradient norm is smaller than the penalty threshold:
S.n fn (w A.
This establishes strict dual feasibility for all j E 8 completing the proof.
Since all KKT conditions satisfied with high probability, the augmented oracle estimator w is sta- tionary point of the penalized program.
C Complete Results for Survival Model
Figure shows that larger variations in C-index are associated with larger censoring rates overall. GMCPNet and GSCADNet achieve comparable results to Oracle-NN while surpassing all other methods, including Oracle-SSF.
D Simulation with Correlated Variables
The simulation study in Scction focuses on independent covariates. However, in real-world applications, particularly in high-dimennsooal settings, the presence of correlations among covariates is common and presents a challenge for feature selection. In this section, we assess the effectiveness of the proposed method using simulated data that incorporates correlated variables.
To be more specifie, extend the high-diiensional scenario described in Section by generating correlated covariate vector, denoted as X~ (0, Σ). The correlation structure is defined using a power decay pattern, where Σij = 0-5 This modification allows us to examine the performance of our method in the presence of correlation among the covariates. Comparing the results of feature selection for independent covariates in Table to the outcomes presented in Table 2 becomes evident that STG and GL.SSONet exhibit larger variations in selected model sizes, along with FNR and FPR in the regression model. This behavior can be attributed to the presence of correlated features. In contrast, the proposed GMCPNet and GSCADNet methods effectively identify relevant variables while maintaining relatively low false positive and negative rates across all models. Furthermore, Figure demonstrates that both GMCPNet and GSCADNet perform comparably to the Oracle-NN method in the regression and survival models, while outperforming other non- oracle approaches the classification model. These findings indicate that the proposed methods exhibit robustness against correlations among covariates in terms feature selection and model prediction.
E
Understanding the Impact of LASSO Bias on Feature Selection
Our numerical study demonstrates that GLASSONet tends to select many noisy variables due to the bias introduced by the LASSO penalty. In contrast, the group-concave regularization used in our proposed framework (e.g,, GMCPNet and CSCADNet) reduces this bias and improves feature selection accuracy.
33

---

Published in Transactions on Machine Learning Research (11/2025)
n=300, d=20,
n=500, d=20,c=0
n=500, d=1000,
0.90
0.90
+
1
0.75
0.75
0.75
Cirdas
0.70
0.70
0.70
0.65
0.65
0.65
0.60
0.45
0.45
0.45
Oracle ASF
Oracle NN RSF
NN
GLASSONet GMCPNet
GSCADNet
Oracle-RS RSF
Oracle NNN RSF
NN
GLASSONet
GMC
GSCADNet
Cracle-RSF
OracleNNi RSF
NN
GLASSONet
GMOPNet
GSCADNet
d=20,
d=20,
d=1000,
0.90
0.90
0.90
1
Ciidex
0.70
0.70
0.70
0.65
0.65
0.65
0.60
0.60
0.60
0.55
0.55
0.55
Oracle-RSF
Oracle NN ASF
NN
GLASSONet GMOPNet
GSCADNet
Oracle-RS RSF
Oracle-NNN RSF
NN
GLASSONet GMct
ESCADNet
Oracle-RS5
OracleNN RSF
NN
GLASSONet
GMOPNet
GSCADNet
n=300, 2= C=0.4
n=500,
2= C=0.4
n=500, d=1000, c=0.4
0.90
0.90
0.90
0.85
0.85
0.85
0.80
Cliddx
0.65
0.65
0.65
0.60
0.60
0.60
055
055
0.5
Cracle-RSP RSF
Oracle NN
RSF
GLASSNnt GMOPNet
GSCADNet
Cracle R5F Oracle 104 RSF
GAASONet
GMOPNet
GSCADNet
Oracle R5F
Oracle NN RSF
GLASONNe GMOtNet
GSCADNet
Figure C-Index of the proposed methods for the survival model outlined in Example 5.3
The dashed line represents the median C-Index of the Oracle-NN, used benchmark for comparison.
34

---

Published in Transactions on Machine Learning Research (11/2025)
Table 2: Feature selection results of STG, GLASSONet, GMCPNet, and GSCADNet using correlated features in high-dimensional scenario (n 500, d = 1000). The False positive rate (FPR %), False negative rate (FNR %), and model size (MS) with standard deviation (SD) in parentheses are displayed.
Regression
Classification
Survival (C 0.2)
Method
FPR, FNR
MS (SD)
FPR, FNR
MS (SD)
FPR, FNR
MS (SD)
STG
8.4, 16.6
86.8(132.6)
1.5, 21.0
18.6(121.1)
-(-)
7
GLASSONet
28.8, 26.6
290.0(144.6)
19.3, 22.4
195.7(116.4)
16.0, 1.9
163.1
GMCPNet
0.1, 13.4
4.0(1.4)
0.2, 13.9
5.5(4.5)
0.0, 0.0
4.1(0.4)
GSCADNet
0.1, 13.2
4.0(1.2)
0.1, 11.8
4.8(2.9)
0.0, 0.0
4.1(0.6)
10
085
0.9
0.8
0.80
0.0
0.6
Acuracy
0.70
Cindex
0.65
0.6
0.60
05
-0.2
0.50
0.4
Oracle
Crocle NN
NN
STG
GLAssoeet
GMCP
GSCDDNet
Cracle RF Oracle- NN RF
NW
STG
GLAsSNett GMCP
GSCADNet
Cracle
Gracle NMN RSF
NN
GLASSONet
GMOPNet
GSCADNet
Figure Prediction scores of the proposed methods for the regression, classification, and survival models (C = 0.2) using correlated features in high-dimensional scenario (n 500, 1000). The dashed lines represent the median of the Oracle-NN, used benchmark for comparison.
35

---

Published in Transactions on Machine Learning Research (11/2025)
To further investigate the impact of LASSO's bias feature selection, propose modified version of GLASSONet, applying a relaxed LASSO approach and terming it relaxed-GTAASSNNe..
The GeLSxeet method follows a two-stage procedure: for each group LASSO parameter X, we first select features using GLASSONet. Then, we refit a standard neural network with only the selected features by setting À O, thereby reducing bias during model fitting. The final model selected based on its predictive performance on a validation set. Our goal to explore whether the Geaxed-GLASSOO mitigate the feature overselection observed in the LASSO-regularizer approach by removing the bias, and ultimately enhance prediction performance.
can
We apply the model, repeating the simulation described in Section 5.1 We compare relaxed-GLASSONNe with GLAS- SONet, as well as our proposed GMCPNet and GSCADNet methods, with results presented Figure g Our results indicate that counts, thereby improving prediction accuracy. Importantly, Gelaxed-GLASSONee performs comparably to GMCPNet and GLASSONet in low-dimensional settings, but its performance declines dimensionality increases. These findings confirm that the LASSO penalty tends to over-select features to compensate for its inherent bias. This overselection can be mitigated by reducing bias through model refitting with only the selected features, leading to more accurate feature selection.
method to synthetic data generated from the XOR-type signal regression
seleets significantly fewer false positives across varying feature
0.6
GLASSONEt
GLASSONEC
GMCPNet
GMCPNet
0.5
GSCADNet
OSCADNet
0.8
Pelaxed GLASSONet
0.6
Relaxed -LAAS---LASSONt
0.4
0.6
0.4
1
FNR
0.4
GLASSONet
GMCPNet
GSCADNet
-0.2
100
300
of features
of features
of features
Figure 9: Simulation results of relaxed-GLASSONe for XOR-Type signal model. The Rt scores, false positive rate (FPR), and false negative rate (FNR) are presented in the left, middle, and right panels, respectively. The contral lines are the mcans while the shaded areas represent standard deviations.
F Implementation Details
F.1 Simulation studies
To ensure fair comparison among all the neural-net-based methods, we adopted a ReLU-activated Multi- Layer Perceptron (MLP) with two hidden layers consisting 10 and 5 units, respectively. The network weights were initialized by sampling from a Gaussian distribution with mean and standard deviation 0.1, while the bias terms were sot following the Xavier initialization technique (Glorot & Bengio 2010). The optimization of the neural networks was performed using the Adam optimizer with base learning rate (LR) of 0.001.
For all the methods falling within the framework of Equation (1) in the paper, we selected the optimal values of X and a from a two-dimensional grid, with À and ranging over 50 and 10 evenly spaced values on a logarithmic scale, respectively. The selection was based on their performance on the validation set, which consisted of 20% of the training set. The parameter scarch ranges are displayed in Table 3 We set À = 0 for NN and Oraclo-NN to deactivate feature selection. For GLASSONet, GMCPNet, and GSCADNet, the number of epochs at Amia set to 2000 for the LD and 200 for the HD scenarios. For all other values of
36

---

Published in Transactions on Machine Learning Research (11/2025)
X, the number of epochs was set to 200 for both LD and HD settings. The number of epochs for NN consistently fixed at 5000.
We employed RF with 1000 decision trees for the model fitting process. We implemented the STG method as described in Yamada et al. (2020) that the LR and regularization parameter were optimized via Optuna with 500 trials, using 10% of the training set as a validation set. The number of epochs was 2000 for each trial.
Table 3: List of the search range for the tuning parameters used in our simulation.
Search range
Param
LD
HD
X
[1e-3, 0.5]
[1-2, 0.5
le-3, 0.1]
[l-2, 0.1
LR (STG)
1c-4,0.1]
[1c44 0.1
(STG)
1e-3, 10
1e-2, 100]
À (LASSONet)
5e-4, 2e-3
5e-4, 2e-3]
F.2
Real Data Example
In the analysis of real data examples, the implementation details remain the same as the HD scenario in the simulation studies, with the following modifications:
For the survival analysis on the CALGB-90401 dataset, we utilized the MLP with two hidden layers, cach consisting of 10 nodes. In hyperparameter tuning, we explored 100 values of ranging from 0.01 to 0.1 for GMCPNet and GSCADNet. Additionally, we increase the number of candidates for to 50.
In the classification task on the MNIST dataset, we adjust the search range of a to [1e-3, 0.1].
The data from CALGB 90401 is available from the NCTN Data Archive at https://neta-data-rreed nci pei.nih.ge. niw.goo//////o. The MNIST dataset was downloaded using the built-in tordavision.datastt...N. interface in PyTorrh, which automatically retrieves the original dataset from Yann LeCun's repository.
References
Wassim Abida, Joanna Cyrta, Glenn Heller, Davide Prandi, Joshua Armenia, Ilsa Coleman, Marcin Cieslik, Mattoo Benelli, Dan Robinson, Eliezer M Van Allen, al. Genomic correlates of clinical outcome in advanced prostale cancer. Proceedings of the National Academyy of Sciences, 116(23):11428 11436, 2019.
Himisha Beltran, Alexander W Wyatt, Edmund C Chedgy, Adam Donoghue, Matti Annala, Evan W Warner, Kevin Beja, Michael Sigouros, Fan Mo, Ladan Fauli, et al. Impact of therapy on genomics and transcrip- tomics in high-risk prostate cancer treated with neoadjuvant docetaxel and androgen deprivation therapy- molecular analysis high-risk pea after neoadjuvant therapy. Clinical Cancer Research, 23(22)::822-6811 2017.
Patrick Breheny and Jian Huang. Coordinate descent algorithms for nonconvex penalized regression, with applications biological feature selection. The annals of applied statistics, 5(1):232, 2011.
William Chan, Navdeep Jaitly, Quoe Le, and Oriol Vinyals. Listen, attend and spell: A neural network for large vocabulary conversational speech recognition. In 2016 IEEE international conference on acoustics, speech and signal processing (ICASSP), pp. 4960 4964. IEEE, 2016.
Bram De Laere, Pieter-Jan van Dam, Tom Whitington, Markus Mayrhofer, Emanuela Henao Diaz, Gert
Van den Fynden, Jean Vandebroek, Jurgen Del-Favero, Steven Van Laere, Lue Dirix, et al. Comprehensive
37

---

Published in Transactions on Machine Learning Research (11/2025)
profiling of the androgen receptor in liquid biopsies from castration-resistant prostate cancer reveals novel intra-ar structural variation and splice variant expression patterns. European wrology, 72(2):192 2000 2017.
Jacob Devlin, Ming- Wei Chang, Kenton Lee, and Kristina Toutanova. Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXio:1810.04800 040000 2018.
Vu Dinh and Lam Si Tung Ho. Consistent feature selection for neural networks via adaptive group lasso.
arXiv preprint ar Xic:2000 00344, 2020.
David L Donoho et al. High-dimensional data analysis: The curses and blessings of dimonsionality, AMS math challenges lecture, (2000):32, 2000.
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties.
Journal of the Americon statistical Association, 96(456):1348 1360, 2001.
Jianqing Fan and Jinchi Lv. A selective overview of variable selection in high dimensional feature space. Statistica Sinica, 20(1):101, 2010.
Jean Feng and Noah Simon. Sparse-input neural networks for high-dimensional nonparametric regression and classification. ar preprint orXf::111.00759) 2017.
Thomas Fischer and Christopher Kranss. Deep learning with long short-term memory networks for financial market predictions. European journal of operational research, 270(2):654-00 669, 2018.
Jerome Friedman, Trevor Hastie, Holger Höfling, and Robert Tibshirani. Pathwise coordinate optimization. The annals of applied statistics, 122):002-33 -32, 2007.
Jerome Friedman, Trevor Hastic, and Rob Tibshirani. Regularization paths for generalized linear models via coordinate descent. Journal of statislical sofbware, 33(1):1, 2010.
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep foedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249- 256. JMLR Workshop and Conference Proceedings, 2010.
Xavier Glorot, Antoine Bordes, and Yoshua Bengio. Deep sparse rectifier neural networks. In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pp. 315-323. JMLR Workshop and Conference Proceedings, 2011.
Pinghua Gong, Changshui Zhang, Zhaosong Lu, Jianhua Huang, and Jieping Ye. A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems. In international conference on machine learning, pp. 37-45. PMLR, 2013.
Alex Graves, Abdel-rahman Mohamed, and Gcoffrey Hinton. Speech recognition with deep recurrent neural networks. In 2013 IEEE international conference on acoustics, speech and signal processing, pp. 6645 6649..
leee, 2013.
Quanquan Gu, Zhenhui Li, and Jiawei Han. Generalized fisher score for feature selection. ar preprint a7X 2012.
Xiao Guo, Hai Zhang, Yao Wang, and Jiang-Lun Wu. Model selection and estimation in high dimensional regression models with group scad. Statistics & Probability Letters, 103:86-92 92, 2015.
Isabelle Guyon and André Elisseeff. An introduction to variable and feature selection. Journal of machine learning research, 3(Mar):1157-1112) 2003.
Susan Halabi, Chen Yen Lin, W Kevin Kelly, Karim Fizazi, Judd W Moul, Ellen B Kaplan, Michael J Morris, and Eric J Small. Updated prognostic model for predicting overall survival first-line chemother- apy for patients with metastatic rastration-resistann prostate cancer. Journal of Clinical Oncology, 32(7):
671, 2014.
38

---

Published in Transactions on Machine Learning Research (11/2025)
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognilion, pp. 770-778, 2016.
Dan Hendrycks and Kevin Gimpel. Gaussian error linear units (gelus). arXiv preprint arXiv:166..00 08406.0841, 2016.
Daniel Hertz, Kouros Owzar, Sherrie Lessans; Claudia Wing, Chen Jiang, William Kevin Kelly, Jai Patel, Susan Halabi, Yoichi Furukawa, Heather E Wheeler, al. Pharmacogenetic discovery in calgb (alliance) 90401 and mechanistic validation of vac14 polymorphism that increases risk of docetaxel- induced neuropathyvac snp predicts docetaxel-induced neuropathy. Clinical Cancer Research, 22(19): 2016.
Jian Huang, Patrick Breheny, and Shuangge Ma. A selective review of group selection in high-dimensional models. Statistical science: a review journal of the Institute of Mathematical Statistics, 27(4), 2012.
Jian Huang, Patrick Breheny, Sangin Lee, Shuangge Ma, and Cun-Hui Zhang. The mnet method for variable selection. Statistica Sinica, pp. 903-923, 2016.
Inaki Iuza, Pedro Larranaga, Rosa Blanco, and Antonio J Cerrolaza. Filter versus wrapper gene selection approaches in dna microarray domains. Artificial intelligence in medicine, 31(2):91-103, 2004.
William Kevin Kelly, Susan Halabi, Michael Carducci, Daniel George, John F Mahoney, Walter M Stadler, Michael Morris, Philip Kantoff, J Paul Monk, Ellen Kaplan, et al. Randomized, double-blind, placebo- controlled phase ii trial comparing docetaxel and prednisone with or without bevacizumab in men with metastatic castration-resistanl prostate cancer: CALGB 90401. Journal of Clivical Oncology, 30(13): 1534-1540, 2012.
Sangjin Kim and Susan Halabi. High dimensional variable selection with error control. BioMed research international, 2016(1):8209453 2016.
Ron Kohavi and George H John. Wrappers for feature subset selection. Arlificial intelligence, 97(1-2): 773-324, 1997.
Daphne Koller and Mchran Sahami. Toward optimal feature selection. Technical report, Stanford InfoLab, 1996.
Alex Krizhevsky, Ilya Sutskever, and Geoffrey Hinton. Imagenet classification with deep convolutional neural networks. Communications of the ACM, 60(6):84-90) 90,,8049--00000 2017.
Johannes Lederer. Statistical guarantees for sparse deep learning. AStA Advances in Statistical Analysis, 108(2):231-258, 2024.
Ismael Lemhadri, Feng Ruan, Louis Abraham, and Robert Tibshirani. Lassonet: A neural network with feature sparsity. The Journal of Machine Learning Research, 22(1):56333-66) 56611 2021.
Yifeng Li, Chih-Yu Chen, and Wyeth W Wasserman. Deep feature selection: theory and application to identify enhancers and promoters. Journal of Computational Biology, 23(5):322 336, 2016.
Yi Lin and Hao Helen Zhang. Component selection and smoothing in multivariate nonparametric regression.
The Annals of Statistics, 32(5):2272 2297, 2006. doi: 10. 10.1144/00003300000007222 URL https:
7/doi org/10. 10. 1214 1214/0990536600000772
Bo Liu, Ying Wei, Yu Zhang, and Qiang Yang. Deep neural networks for high dimension, low sample size data. In ICAI, pp. 2287-2293, 2017.
Joaquin Matco, Suzanne Carreira, Shahneen Sandhu, Susana Miranda, Helen Mossop, Raquel Porez-Lopez, Daniel Nava Rodrigues, Dan Robinson, Aurelius Omlin, Nina Tunariu, et al. Dna-repair defects and olaparib in metastatic prostate cancer. New England Journal of Medicine, 373(18):1697-1777 2015.
39

---

Published in Transactions on Machine Learning Research (11/2025)
Lukas Meier, Sara Van de Geer, and Peter Btihlmann. High-dimensional additive modeling. The Annals of Stalistics, 37(6B):3779 -821, 2009.
Juan Miguel Mosquera, Himisha Beltran, Kyung Park, Theresa Y MacDonald, Brian D Robinson, Scott T Tagawa, Sven Perner, Tarek A Bismar, Andreas Erbersdobler, Rajiv Dhir, et al. Concurrent aurka and myen gene amplifications are harbingers of lethal treatmentrelated neuroendocrine prostate cancer. Neoplasia, 15(1):1 IN4, 2013.
Vinod Nair and Geoffrey E Hinton. Rectified linear units improve restricted boltzmann machines. Proceedings of the 27th international conference machine learning (ICMI-10), pp. 807-814, 2010.
Yu Nesterov. Gradient methods for minimizing composite functions. Mathemalical Programmming 140(1): 125-161, 2013.
Pradeep Ravikumar, John Lafferty, Han Liu, and Larry Wasserman. Sparse additive models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71 (55):1000-11300 2009.
Dan Robinson, Eliezer M Van Allen, Yi-Mi Wu, Nikolaus Schultz, Robert Lonigro, Juan-Miguu Mosquera, Bruce Montgomery, Mary-Ellen Taplin, Colin C Pritchard, Gerhardt Attard, et al. Integrative clinical genomics of advanced prostate cancer. Cell, 161(5):1215-122) 2015.
Simone Scardapane, Danilo Comminiello, Amir Hussain, and Aurelio Uncini. Group sparse regularization for deep neural networks. Neurocomputing, 241:81-80 2017.
Nicolas Städler, Peter Bühhlaan, and Sara Van De Geer. Ll-penalization for mixture regression models. Test, 19(2):209-25) 256, 2010.
Takayuki Sumiyoshi, Xiaofei Wang, Evan W Warner, Andrea Sboner, Matti Annala, Michael Sigouros, Kevin Beja, Kei Miauno, Shengyu Ku, Ladan Fazli, James Eastham, Mary-Ellen Taplin, Jeffrey Simko, Susan Halabi, Michael J Morris, Martin E Gleave, Alexander W Wyatt, and Himisha Beltran. Molecular features of prostate cancer after neoadjuvant therapy in the phase 3 calgb 90203 trial. JNCI: Journal of the National Cancer Institute, 116(1):115 1166))111-12,,,, 2024.
Yan Sun, Qifan Song, and Faming Liang. Consistent sparse deep learning: Theory and computation. Jowrnal of the American Statistical Association, 117(540):198111995, 2022.
Jiliang Tang, Salem Alelyani, and Huan Liu. Feature selection for classification: A review. Data classification: Algorithms and applications, pp. 37, 2014.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methonlological) pp. 267 288, 1996.
Hajime Uno, Tianxi Cai, Lu Tian, and Lee-Jen Wei. Evaluating prediction rules for t-year survivors with censored regression models. Journal of the American Statistical Association, 102(478):527-533 2007.
Peter M Visscher, Matthew A Brown, Mark McCarthy, and Jian Yang. Five years of gwas discovery. The American Journal of Human Genetics, 90(1):7-24 2012.
Fengrong Wei and Hongxiao Zhu. Group coordinate descent algorithms for nonconvex penalized regression. Computational statistics & data analysis, 56(2):316-32) 2012.
Alexander W Wyatt, Arun A Azad, Stanislav V Volik, Matti Annala, Kevin Beja, Brian McConeghy, Anne Hacgert, Evan W Warner, Fan Mo, Sonal Brahmbhatt, et al. Genomic alterations in cell-free dna and enzalutamide resistance in castration-resistant prostate cancer. JAMA oncology, 2(12):1598-1606 2016.
Makoto Yamada, Wittawat Jitkrittum, Leonid Sigal, Eric P Xing, and Masashi Sugiyama. High-dimensional feature selection by feature-wise kernelized lasso. Neural computation, 26(1):185 207, 2014.
Yutaro Yamada, Ofr Lindenbaum, Sahand Negahban, and Yuval Kluger. Feature selection using stochastic gates. In International Conference Machine Learning, pp. 10648-10651 PMT.R, 2020.
40

---

Published in Transactions on Machine Learning Research (11/2025)
Tom Young, Devamanyu Hazarika, Soujanya Poria, and Erik Cambria. Recent trends in deep learning based natural language processing. seee Compulational intelligenCe magazine, 13(3):55 75, 2018.
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 68(1):49 67, 2006.
Lingmin Zeng and Jun Xie. Group variable selection via scad-1 2. Statistics, (5(1):49-000) 66, 2014.
Cun-Hui Zhang et al. Nearly unbiased variable selection under minimax concave penalty. The Annals of statistics, 38(2):894-942, 2010.
Hui Zou. The adaptive lasso and oracle properties. Journal of the American statistical association, (476):1418 1429, 2006.
Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. Journal of the royal statistical society: series (statistical methodology), 67(2):301 320, 2005.
41


Reviews:
[{"review_id": "uQo0uVxhlc", "reviewer": "Reviewer_tuMf", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength :\n- The article is well written and easy to follow. \n- The way to specify a group for sparse feature selection is a novelty compared to existing works, such as Feng & Simon 2017 on regularized neural networks.\n- Empirical validation across diverse data types. \n\nWeakness :\n- Although the authors claim that a stable optimization algorithm is proposed in the introduction, it is not so clear what this means. In particular, there is a lack of convergence analysis for the optimization algorithm used in the numerical results.\n\n- The baseline methods chosen in real data example in Section 5 seem to be too weak. It would be better to discuss whether there are other baseline methods to consider.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This article proposes a sparse input feature selection method to improve the performance of neural networks. \nIn the litterature of feature selection, many methods are developed for linear models\nand they may not able to capture complex input-output relations. To model complex relationships, \nregularized neural networks are studied, but existing works are limited either to regression-only problems\nor to non-sparse feature selection. This article proposes a way to encounter these limitations by using group concave regularization \nsimilar to group LASSO on both synthetic and real datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Fig 1. Please explain what smoother solution paths mean. It is not clear as the x-axis uses a different choice of log (lambda). It would be better to use the same x-axis to compare. Also what is the difference between the curves with the same color in Fig 1? \n- According to Section 4, Adam method is used in numerical results as the optimization algorithm. Could you clarify in what sense it is stable and why?\n- The comparison with SOTA in Section 5 seems to be of limited impact. In the MNIST case, the test accuracy is not so high compared to convolutional neural networks. At the current test accuracy level, linear SVM could also learn something like Fig 6. Why is this method not compared? What is the SOTA AUC performance in the survival dataset of Section 5.1 ? I think that these SOTA performances should be mentioned as a too large gap suggests that the considered neural networks in the article are not able to capture complex input-output relationships (which could be future work).", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "VB7xviXv1c", "reviewer": "Reviewer_t1FU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, this work seems well motivated and technically sound. It is also clear to read.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work provides a new framework for feature selection and prediction in high-dimensional settings, for several tasks and losses (regression, classification, survival analysis). The authors propose and compare several losses for the selection of features, and illustrate their properties in several simulated and real data situations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "as the authors point out themselves, there is a strong lack of theoretical grounds for the proposed method, and this should be compensated for with a much stronger empirical analysis. Here are two main suggestions\n- include baselines, both from the framework that predicts and selects in a single algorithm, but also that does variable selection in a post-hoc way (variable importance, knockoffs, there is a huge litterature). The authors may look into the International Seminar on Selective Inference for references https://sites.google.com/view/selective-inference-seminar/home?authuser=0\n- For all methods and experiments, the authors should include more metrics, in particular, proper metrics. Accuracy or the c-index are known to provide a very limited view of the problem (https://arxiv.org/abs/2506.02075 and https://arxiv.org/pdf/2412.10288 are good introductions to this question). \n\nThe objective would be to show empirically that the proposed approaches are SOTA both in prediction and selection; otherwise, it is not clear what the point is of using this combined approach compared to existing methods. A study on the computing cost can also be done, which could be a strong point of this method.", "broader_impact_concerns": "Without strong claims of performance compared to the state of the art, nor theory, the impact of this work is limited.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5JgXP2XK35", "reviewer": "Reviewer_Vnaj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: The numerical experiments are detailed. \n\nWeakness: No theoretial guarantees are provided for the proposed method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes conducting input variable selection using a combined penalty, which is a combination of an L2-penalty and a group concave penalty. While empirical results are presented, no theoretical results are provided. The overall contribution of the paper appears to be incremental, considering the sparse input work [1]  cited and other existing sparse deep learning methods, see e.g., [2]. \n\n[1] Feng, J. and Simon, N. (2017). Sparse input neural networks for high-dimensional nonparametric regression and classification. arXiv:1711.07592.\n\n[2] Sun, Y., Song, Q., & Liang, F. (2021). Consistent Sparse Deep Learning: Theory and Computation, arXiv:2102.13229.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "1, In recent years, theoretical guarantees have been developed for sparse deep neural networks, such as in references [2] and [3]. I believe these works merit discussion in the paper.\n\n[2] Sun, Y., Song, X., & Liang, F. (2021). Consistent Sparse Deep Learning: Theory and Computation, arXiv:2102.13229.\n[3] Lederer, J. (2022). Statistical guarantees for sparse deep learning, arXiv:2212.05427.\n\n2. The impacts of α and 𝜆 on the performance of the proposed method need to be assessed. For example, if α=0, the method may fail to select any meaningful input variables, given the universal approximation ability of deep neural networks. I understand that the authors determine their values using a validation method. Are there any insights regarding the choice of their values?\n\n3. The comparsions with more baseline methods might be helpful to strengrgen the paper due to the lack of theoretical guarantees.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)