End-to-end Feature Selection Approach for Learning Skinny Trees

Shibal Ibrahim, Kayhan Behdin, Rahul Mazumder Massachusetts Institute of
Technology, Cambridge, MA, USA

###### Abstract

Joint feature selection and tree ensemble learning is a challenging task.
Popular tree ensemble toolkits e.g., Gradient Boosted Trees and Random Forests
support feature selection post-training based on feature importances, which
are known to be misleading, and can significantly hurt performance. We propose
Skinny Trees: a toolkit for feature selection in tree ensembles, such that
feature selection and tree ensemble learning occurs simultaneously. It is
based on an end-to-end optimization approach that considers feature selection
in differentiable trees with Group $\ell_{0}-\ell_{2}$ regularization. We
optimize with a first-order proximal method and present convergence guarantees
for a non-convex and non-smooth objective. Interestingly, dense-to-sparse
regularization scheduling can lead to more expressive and sparser tree
ensembles than vanilla proximal method. On 15 synthetic and real-world
datasets, Skinny Trees can achieve $1.5\\!\times\\!-~{}620~{}\\!\times\\!$
feature compression rates, leading up to $10\times$ faster inference over
dense trees, without any loss in performance. Skinny Trees lead to superior
feature selection than many existing toolkits e.g., in terms of AUC
performance for 25% feature budget, Skinny Trees outperforms LightGBM by
$10.2\%$ (up to $37.7\%$), and Random Forests by $3\%$ (up to $12.5\%$).

##  1 Introduction

Decision trees have been popular in various machine learning applications [20,
16] for their competitive performance, interpretability, robustness to
outliers, and ease of tuning [26]. Many scalable toolkits for learning tree
ensembles have been developed [11, 16, 39, 54]. While these toolkits are good
candidates for building tree ensembles, they do not enable feature selection
during training.

Feature selection is a fundamental problem in machine learning and has seen
widespread usage across various real-world tasks [68, 12, 44]. The popular
tree ensemble toolkits only allow selecting informative features post-training
based on feature importances, which produce misleading feature rankings for
feature selection [65, 9, 75]. Recently, there has been some work on
optimization-based sparse feature selection in trees. For example, the local-
search-based TAO algorithm [74] imposes $\ell_{1}$ penalty in oblique decision
trees, such that each node is a hyperplane involving only few features.
However, this only achieves node-level feature selection, and is sub-optimal
for tree-level or ensemble-level feature selection as we show in Sec. 7.1.
[45] proposed ControlBurn that considers a group-lasso type penalty for
feature selection on a pre-trained forest. This can be viewed as a two-stage
procedure (unlike an end-to-end training procedure we propose here), where
feature selection occurs after a tree ensemble is formed using all original
features. While these methods serve as promising candidates for feature
selection, these works highlight that identifying relevant features while
learning compact trees remains an open challenge—an avenue we address in this
work.

Many problems assume that the data is already stored in datasets and available
without charge. However, in practice, data are not free in real-world
applications. It costs both money and time to collect feature values [49, 71].
Under this context, one would like to collect a compact set of features to
reduce experimental costs. Additionally, selecting a compact set of relevant
features can lead to enhanced interpretability [59], faster inference,
decreased memory footprint, and even improved model generalization on unseen
data [14].

In this paper, we propose an end-to-end learning framework for global feature
selection in trees such that feature selection is performed jointly with the
tree ensemble optimization in _one_ stage. Our framework builds upon
differentiable (a.k.a. soft) tree ensembles [38, 41, 27, 35]. These works do
not address feature selection in trees. In this paper, we use a combinatorial
penalty function for feature selection in tree ensembles. We use a group
$\ell_{0}-\ell_{2}$ penalty function. While group $\ell_{0}-\ell_{2}$ penalty
has been found to be useful in recent work on high-dimensional linear [29] and
additive models [29, 36], training for an ensemble of trees presents a
different set of challenges. To obtain models with superior generalization and
sparsity properties, we need to pay special attention to dense-to-sparse
training, which is contrary to what is employed in linear/additive models
above. We demonstrate that our end-to-end learning approach leads to superior
feature selection in tree ensembles than that achieved by popular toolkits.

Contributions. As discussed earlier, joint feature selection and tree ensemble
learning is a challenging task. Existing tree toolkits perform feature
selection post-training based on feature importances, which are misleading
measures for feature selection. They can also hurt performance in (a) settings
with $N<p$ and (b) settings with correlated features (see Sec. 6). In this
context, our contributions can be summarized as:

  * •

We propose a joint optimization approach, where feature selection is performed
simultaneously with tree ensemble learning. This joint training makes it
different from post-training feature selection in trees. We propose an end-to-
end feature selection approach for learning sparse (differentiable) tree
ensembles with a group $\ell_{0}$-based regularizer.

  * •

Our algorithmic workhorse is based on proximal mini-batch gradient descent
(GD). We also discuss the convergence properties of our approach in the
context of a nonconvex and nonsmooth objective. Interestingly, augmenting
first-order optimization with _dense-to-sparse_ regularization scheduling
leads to sparser tree ensembles with improved test AUC.

  * •

We introduce a new toolkit: Skinny Trees. We consider 15 synthetic and real-
world datasets, showing that our toolkit can lead to superior feature
selection and test AUC compared to popular toolkits. In particular, for 25%
feature budget, Skinny Trees outperforms LightGBM by $10.2\%$ (up to
$37.7\%$), XGBoost by $3.1\%$ (up to $17.4\%$), and Random Forests by $3\%$
(up to $12.5\%$) in test AUC. Skinny Trees will be open-sourced.

##  2 Related Work

Trees. Learning tree ensembles is a hard combinatorial optimization problem
[32]. Historically, tree ensembles have been trained with bagging [11] or
boosting methods [26]. These have led to various popular toolkits, e.g.,
Random Forests [11], Gradient Boosted Trees [16, 39, 54]. Another line of
approach, most related to the approach we pursue here, considers
differentiable relaxations at the split nodes [38] and performs end-to-end
learning of trees with first-order methods [41, 27]. Joint training of soft
trees results in compact representations (i.e., fewer trees) compared to
boosting-based procedures [27]. Despite the success of both lines of
approaches, interesting challenges remain: earlier approaches do not perform
feature selection while training in an end-to-end fashion. Next, we summarize
popular feature selection methods.

Feature selection. Feature selection methods are classified into three major
categories: (i) Filter methods remove irrelevant features before learning a
model. These methods rely on statistical measures to filter features using a
per-feature relevance score [5, 52, 21, 63, 62, 15]. These methods look at the
marginal effect of features without considering the joint effect of feature
interactions. (ii) Wrapper methods [40, 64, 76, 57, 3, 60] use the outcome of
a model to determine the relevance of each feature. These methods require
recomputing the model for each subset of features and is therefore
computationally expensive. Feature-importance based selection from pre-trained
forests/ensembles [11, 16, 39, 54] also fall under this category. Similarly,
[45] proposed ControlBurn, which formulates an optimization problem with
group-lasso regularizer to perform feature selection on a pre-trained forest.
They rely on a commercial solver to optimize their formulation. (iii) Embedded
methods simultaneously learn the model and the relevant subset of relevant
features. Notable among these methods include $\ell_{0}$-based methods [28,
29, 36, 33] and their relaxations via lasso-type methods [67, 56, 72] in the
linear and additive model settings. Embedded nonlinear feature selection
methods have also been explored for neural networks. [17] consider a
cardinality-constrained feature selection approach using an active-set style
method. Others consider grouped lasso type methods [61, 23, 18, 43] or
consider reparameterizations of $\ell_{0}$ with stochastic gates, which act as
“soft” masking of features during training [46, 69].

We propose an embedded approach that simultaneously performs feature selection
and tree ensemble learning. This can have important implications for
compression, efficient inference, and/or generalization.

Organization. The rest of the paper is organized as follows. Section 3
summarizes relevant preliminaries. Section 4 presents a formulation for
feature selection in soft tree ensembles. Section 5 discusses our optimization
algorithm and its convergence properties. Later, we discuss a scheduling
approach that can result in more compact and higher quality feature selection.
Sections 6 and 7 perform experiments on a combination of 15 synthetic and
real-world datasets to highlight the usefulness of our proposals.

##  3 Preliminaries

We learn a mapping $\bm{f}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{c}$, from
input space $\mathcal{X}\subseteq\mathbb{R}^{p}$ to output space
$\mathcal{Y}\subseteq\mathbb{R}^{c}$, where we parameterize function $\bm{f}$
with a soft tree ensemble. In a regression setting $c=1$, while in multi-class
classification setting $c=C$, where $C$ is the number of classes. Let $m$ be
the number of trees in the ensemble and let $\bm{f}^{j}$ be the $j$th tree in
the ensemble. We learn an additive model with the output being sum over
outputs of all the trees: $\bm{f}(\bm{x})=\sum_{j=1}^{m}\bm{f}^{j}(\bm{x})$
for an input feature-vector $\bm{x}\in\mathbb{R}^{p}$. A summary of the
notation can be found in Table S1 in Supplement.

Contrary to classical trees, soft trees allow significantly more flexibility
in catering to more application-focused loss functions [35], sparse routing in
Mixture of Experts [34] etc. A soft tree is a variant of a classical decision
tree that is differentiable, so learning can be done using gradient-based
methods. It was proposed by [38], and further developed by [41, 27, 35] for
end-to-end optimization. Soft trees typically perform soft routing, i.e., a
sample is fractionally routed to all leaves. However, they can also perform
hard routing [27, 34]. Soft trees are based on hyperplane splits, where the
routing decisions rely on a linear combination of the features. Particularly,
each internal node is associated with a trainable weight vector that defines
the node’s hyperplane split. Formulations [35] for training soft tree
ensembles more efficiently have been proposed, which exploit tensor structure
of a tree ensemble. A summary of soft trees and tensor formulation is in
Supplement Sec. S2.

##  4 Problem Formulation

Feature selection plays a ubiquitous role in modern statistical regression,
especially when the number of predictors is large relative to the number of
observations. We describe the problem of global feature selection. We assume a
data-generating model $p({\mathbf{x}};{\mathbf{y}})$ over a $p$-dimensional
space, where ${\mathbf{x}}\in\mathbb{R}^{p}$ is the covariate and
${\mathbf{y}}$ is the response. The goal is to find the best function
$\bm{f}({\mathbf{x}})$ for predicting ${\mathbf{y}}$ by minimizing:

$\displaystyle{\min}_{\bm{f}\in\mathcal{F},\mathcal{Q}}~{}~{}\mathbb{E}[L({\mathbf{y}},\bm{f}({\mathbf{x}}_{\mathcal{Q}}))]$
(1)

where $\mathcal{Q}\subseteq\\{1,2,\cdots,p\\}$ is an unknown (learnable)
subset of features of size at most $K$, $\bm{f}$ is a learnable non-parametric
function from the function class $\mathcal{F}$, and
$L:\mathbb{R}^{p}\times\mathbb{R}^{c}\rightarrow\mathbb{R}$ is a loss
function. The principal difficulty in solving (1) lies in the joint
optimization of $(f,\mathcal{Q})$—the number of subsets $\mathcal{Q}$ grows
exponentially with $p$. In addition, the function class $\mathcal{F}$ needs to
exhibit strong expressive power (here, $\mathcal{F}$ is the class of soft tree
ensembles with fixed ensemble size $m$ and depth $d$).

Figure 1: Illustration of Skinny Trees. Each horizontal slice
$\bm{\mathcal{W}}_{k,:,:}$ depicts a single feature. White slices indicate
features filtered out by the ensemble while training. Each vertical slice
(along the depth of the page), $\bm{\mathcal{W}}_{:,:,j}=\bm{W}_{j}$
corresponds to parameters in $j$-th (splitting) supernode (blue circles) in
the ensemble, eventually producing the routing decisions. The red squares
depict leaf nodes. $S(\cdot)$ denotes an activation function, which can be
Sigmoid [38] or Smooth-Step [27].

In the context of tree ensembles, the goal is to select a small subset of
features across all the trees in the ensemble for global feature selection.
More specifically, we consider the framework with response $\bm{y}$ and
prediction $\bm{f}(\bm{x};\bm{\mathcal{W}},\bm{\mathcal{O}})$, where the
function $\bm{f}$ is parameterized by learnable hyperplane parameters
$\bm{\mathcal{W}}\in\mathbb{R}^{p,m,|\mathcal{I}|}$ across all the split nodes
and learnable leaf parameters $\bm{\mathcal{O}}$ across all trees. Note that
$m$ is the number of trees in the ensemble and $|\bm{\mathcal{I}}|$ represents
the number of (split) nodes in each tree. This parameterization of soft trees
points to a group structure in $\bm{\mathcal{W}}$, where the whole slice
$\bm{\mathcal{W}}_{k,:,:}$ in the tensor formulation has to be zero to
maintain feature sparsity both across all split nodes in each tree and across
the trees in the ensemble — see Fig. 1. This is a natural feature-wise group
structure in soft trees and allows adaptation of the grouped selection problem
in linear models [48, 7] to soft trees.

Mixed Integer Problem (MIP) Formulation. Let us consider the tensor
$\bm{\mathcal{W}}$. Feature selection is described via binary variables. This
leads to the regularized loss function:

$\displaystyle\min_{\bm{\mathcal{W}},\bm{\mathcal{O}},{\bm{z}}}$
$\displaystyle~{}\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})]+\lambda_{0}\sum_{k\in[p]}z_{k},$
(2) s.t.
$\displaystyle~{}||\bm{\mathcal{W}}_{k,:,:}||(1-z_{k})=0,z_{k}\in\\{0,1\\}~{}~{}k\in[p],$

where, the binary variable $z_{k}$ controls whether the $k$-th feature is on
or off via the constraint $||\bm{\mathcal{W}}_{k,:,:}||(1-z_{k})=0$.
$\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))]:=(1/N)\sum_{n\in[N]}L(\bm{y}_{n},\bm{f}(\bm{x}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}}))$
is the empirical loss; and $\lambda_{0}$ is regularization strength. Note MIP
formulations [6] can also be setup with classical trees under feature
selection, but these are intractable from optimization standpoint.

Unconstrained Reformulation of Problem 2. For training purposes, we consider
an unconstrained version of (2) involving variables
($\bm{\mathcal{W}},\bm{\mathcal{O}}$), via the grouped $\ell_{0}$ (pseudo)
norm so that it is amenable to end-to-end training via first order methods —
more details in Sec. 5. [29] has highlighted overfitting with vanilla version
of $\ell_{0}$ (and group $\ell_{0}$) penalization in the context of
nonparametric additive models. One potential solution to ameliorate this
problem is to include an additional ridge regularization for shrinkage [47,
29, 36]. We consider the following group $\ell_{0}-\ell_{2}$ regularized
problem:

$\displaystyle{\min}_{\bm{\mathcal{W}},\bm{\mathcal{O}}}$
$\displaystyle~{}~{}\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))]$
(3)
$\displaystyle+\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]+(\lambda_{2}/m|\mathcal{I}|)||\bm{\mathcal{W}}||_{2}^{2}$

where $1[\cdot]$ is the indicator function, $\lambda_{0}\geq 0$ controls the
number of features selected, and $\lambda_{2}$ is a non-negative
regularization parameter that controls the amount of shrinkage of each group.
We normalize $\lambda_{2}$ by the product $m|\bm{\mathcal{I}}|$ for
convenience in hyperparameter tuning.

##  5 End-to-end learning for sparse trees

We propose a fast approximate algorithm to obtain high quality local solutions
to Problem (3). As we discuss here, this algorithm admits a simple and
efficient implementation. We use a proximal (mini-batch) gradient based
algorithm. Training sparse soft trees involves two operations. First, a
vanilla mini-batch GD step is applied to all model parameters. Then, a
proximal operator is applied to the hyperplane parameters $\bm{\mathcal{W}}$,
which minimizes a local approximation of the loss with the $\ell_{0}$ penalty
unchanged. This sequential operation on top of backpropagation makes the
procedure simple to implement in popular ML frameworks e.g. Tensorflow [1].

###  5.1 Proximal mini-batch gradient descent

We first present the proximal mini-batch GD algorithm for solving Problem (3)
in Algorithm 1. We also discuss computation of the Prox operator in line 7 of
Algorithm 1. Prox finds the global minimum of the optimization problem:

$\displaystyle\bm{\mathcal{W}}^{(t)}=\text{argmin}_{\bm{\mathcal{W}}}$
$\displaystyle~{}~{}(1/2\eta)||\bm{\mathcal{W}}-\bm{\mathcal{Z}}^{(t)}||_{2}^{2}$
$\displaystyle+\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq 0]$
(4)

where
$\bm{\mathcal{Z}}^{(t)}=\bm{\mathcal{W}}^{(t-1)}-\eta\nabla_{\bm{\mathcal{W}}}h$.
As it turns out, the proximal operator (Prox) is decomposable across features.
This leads to a feature-wise hard-thresholding operator given by:

$\displaystyle
H_{\eta\lambda_{0}}(\bm{\mathcal{Z}}^{(t)}_{k,:,:})=\bm{\mathcal{Z}}^{(t)}_{k,:,:}\odot
1\left[||\bm{\mathcal{Z}}^{(t)}_{k,:,:}||\geq\sqrt{2\eta\lambda_{0}}\right].$
(5)

1:Data: $X,Y$; Hyperparameters: $\lambda_{0},\lambda_{2}$, epochs, batch-size,
learning rate ($\eta$);

2:Initialize ensemble with $m$ trees of depth $d$ ($|\mathcal{I}|=2^{d}-1$):
$\bm{\mathcal{W}},\bm{\mathcal{O}}$

3:for $\text{epoch}=1,2,\ldots,\text{epochs}$ do

4: for $batch=1,\ldots,N/\text{batch-size}$ do

5: Randomly sample a batch: $\bm{X}_{batch},\bm{Y}_{batch}$

6: Compute gradient of loss $h$ w.r.t. $\bm{\mathcal{O}}$, W.

7: Update leaves:
$\bm{\mathcal{O}}\leftarrow\bm{\mathcal{O}}-\eta\nabla_{\bm{\mathcal{O}}}h$

8: Update hyperplanes:
$\bm{\mathcal{W}}\leftarrow\textit{Prox}(\bm{\mathcal{W}}-\eta\nabla_{\bm{\mathcal{W}}}h,\eta,\lambda_{0}))$

9: end for

10:end for

where
$h=\hat{\mathbb{E}}[L(\bm{Y}_{batch},\bm{F}_{batch}]+(\lambda_{2}/m|\mathcal{I}|)||{\bm{\mathcal{W}}}||_{2}^{2}$

Algorithm 1 Proximal Mini-batch Gradient Descent for Optimizing (3).

We use the feature-wise separability for faster computation. The cost of Prox
is of the order $\mathcal{O}(v)$, where $v$ is the total hyperplane parameters
being updated (i.e. $v=pm|\mathcal{I}|$). This overhead is negligible compared
to the computation of the gradients with respect to the same parameters. Also,
implementing the optimizer is straightforward in standard deep learning APIs.

To our knowledge, we are the first to study iterative-hard-thresholding based
approach for group $\ell_{0}$-based nonlinear feature selection in soft tree
ensembles. Note that [17] considers group $\ell_{0}$ based cardinality
constrained formulation for feature selection in neural networks. However,
their active-set style approach is very different from our iterative-hard-
thresholding based approach. Their toolkit also appears to be up to
$900\times$ slower than our toolkit. Also, to our knowledge, in a related
context in neural network pruning, modifications of iterative-hard-
thresholding based approaches [37, 53] have appeared for individual weight
pruning in neural networks, but not for a grouped problem as the one we
pursue.

###  5.2 Convergence Analysis of Algorithm 1

In this section, we analyze the convergence properties of Algorithm 1. For
simplicity, we assume that the outcomes are scalar, i.e. $c=1$ and consider
the least squares loss (see Assumption (A2) below). We also analyze the full-
batch algorithm, i.e., we assume $\text{batch-size}=N$. Extending our result
to vector outputs and a mini-batch algorithm is straightforward and therefore
omitted here. Before stating our formal results, we discuss our assumptions on
the model.

  1. (A1)

(Activation function, $S(\cdot)$). $S(\cdot)$ is a differentiable piece-wise
polynomial function,

$S(x)=\begin{cases}1&\mbox{if}~{}~{}~{}x>\theta\\\
p(x)&\mbox{if}~{}~{}~{}-\theta\leq x\leq\theta\\\
0&\mbox{if}~{}~{}~{}x<-\theta\\\ \end{cases}$

for some $\theta>0$ and a polynomial function $p(x)$. We assume the derivative
of $S(x)$ is continuous.

  2. (A2)

(Loss function, $L$). We use the least squares loss $L(x,y)=(x-y)^{2}/2$.

  3. (A3)

(Solutions). There exists a numerical constant $B>0$ only depending on the
data, such that $\|\bm{\mathcal{O}}\|_{2}\leq B$ for all iterations in
Algorithm 1.

Assumption (A1) encompasses a general class of activation functions.
Particularly, the smooth-step function [27] used in soft trees falls under
this class. The least squares loss in Assumption (A2) is a standard loss for
regression problems. Finally, we argue the boundedness of leaf weights in
Assumption (A3) is a weak assumption as the data is bounded. Alternatively,
one can add additional projection steps for leaf weights in Algorithm 1,
ensuring the boundedness, and remove Assumption (A3).

Our main result in this section is stated in Theorem 1.

######  Theorem 1.

Let $\lambda_{2}>0$ and suppose Assumptions (A1), (A2) and (A3) hold. Then:

  1. 1.

There exists a sufficiently small $\eta>0$ for which Algorithm 1 (full-batch)
is a descent algorithm with non-increasing objective.

  2. 2.

The sequence of $\bm{\mathcal{W}}$ generated from Algorithm 1 is bounded.

  3. 3.

The sequence of $\bm{\mathcal{W}},\bm{\mathcal{O}}$ generated from Algorithm 1
converges if $\eta>0$ is chosen as in Part 1.

Theorem 1 shows that Algorithm 1 is indeed a convergent descent method for a
properly selected learning rate. The proof of Theorem 1 is presented in
Supplement S3. Here, we note that Problem (3) is non-convex and non-smooth.
Moreover, the activation function and therefore
$\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))]$
are not twice differentiable everywhere. These lead to technical challenges in
proof, as discussed in the Supplement.

###  5.3 Dense-to-Sparse Learning (DSL)

Prior work in feature selection [43] suggest a multi-stage approach: Train a
dense model completely and then learn a progessively sparser model. At each
sparsity level, the model is trained to convergence. This approach appears to
effectively leverage favorable generalization properties of the dense solution
and preserves them after drifting into sparse local minima [43]. For example,
this approach appears to be superior than sparsely learning from scratch.
However, this approach can be expensive as it requires learning a dense model
completely before the sparser training even begins. For more compact feature
selection (higher sparsity settings), this approach is also time-consuming.

Lung | Madelon
---|---
|
|
Figure 2: Trajectory of validation loss and feature sparsity during training with dense-to-sparse learning. Random Forests | XGBoost |  Skinny Trees (ours)
---|---|---
N=100 | N=100 | N=100
|  |
N=200 | N=200 | N=200
|  |
N=1000 | N=1000 | N=1000
|  |
Figure 3: Features selected by Random Forests, XGBoost and Skinny Trees for
different sample sizes

To reduce the computational overhead of the above approach, we propose dense-
to-sparse learning in a single stage. To perform dense-to-sparse learning
(DSL) for feature selection, we annealed the sparsity-inducing penalty
$\lambda_{0}$ from small to large ($0\rightarrow\gamma$) during the course of
training. Our approach requires much lesser computational budget to reach the
desired sparsity level as that required by [43]. We resort to an exponential
annealing scheduler of the form: $\lambda_{0}=\gamma(1-\text{exp}(-s*t))$,
where $\gamma$ is the largest regularization penalty that leads to a fully
sparse soft tree, $s$ controls the rate of increase of the regularization
penalty and $t$ denotes the iteration step.

We show the validation loss trajectory and the number of features selected
during the course of training with dense-to-sparse learning in Figure 2. We
empirically observed this scheduler to result in better out-of-sample
performance and more compact trees in terms of feature selection (see Figure 4
in Sec. 7.5).

##  6 Simulation

We first evaluate our proposed method using data with correlated features. In
real-world, high-dimensional datasets, many features are correlated. Such
correlations introduce challenges for feature selection. For instance, if
there is set of highly correlated features, then it is not clear which to
select. Existing tree ensemble toolkits, e.g., XGBoost and Random Forests,
based on feature importances, may produce misleading results — any of the
correlated features can work as a splitter, the feature importance scores can
get distributed among the correlated features. Below, we consider a correlated
feature setting and demonstrate the strong performance of Skinny Trees in
terms of true support recovery because of end-to-end learning.

Table 1: Test MSE, feature sparsity and support recovery metrics (F1-score)
for a linear setting with correlated design matrix. Skinny Trees outperforms
feature-importance based methods across all metrics.

$\sigma$ | $p$ | $N$ | Model | Test MSE | #features | F1-score
---|---|---|---|---|---|---
0.7 | 512 | 100 | Random Forests | 6.49 $\pm$ 0.19 | 79 $\pm$ 17 | 0.21 $\pm$ 0.02
XGBoost | 8.65 $\pm$ 0.27 | 32 $\pm$ 9 | 0.18 $\pm$ 0.03
Skinny Trees |  0.65 $\pm$ 0.12 |  12 $\pm$ 1 |  0.86 $\pm$ 0.04
200 | Random Forests | 4.90 $\pm$ 0.15 | 40 $\pm$ 8 | 0.35 $\pm$ 0.03
XGBoost | 5.97 $\pm$ 0.12 | 110 $\pm$ 21 | 0.18 $\pm$ 0.02
Skinny Trees |  0.34 $\pm$ 0.00 |  11 $\pm$ 1 |  0.89 $\pm$ 0.03
1000 | Random Forests | 2.97 $\pm$ 0.03 | 11 $\pm$ 1 | 0.84 $\pm$ 0.02
XGBoost | 1.81 $\pm$ 0.02 | 24 $\pm$ 1 | 0.50 $\pm$ 0.02
Skinny Trees |  0.26 $\pm$ 0.00 |  8 $\pm$ 0 |  1.00 $\pm$ 0.00
0.5 | 256 | 100 | Random Forests | 6.24 $\pm$ 0.13 | 42 $\pm$ 11 | 0.35 $\pm$ 0.03
XGBoost | 7.93 $\pm$ 0.27 | 35 $\pm$ 10 | 0.25 $\pm$ 0.02
Skinny Trees |  0.45 $\pm$ 0.06 |  10 $\pm$ 1 |  0.89 $\pm$ 0.03
200 | Random Forests | 4.40 $\pm$ 0.13 | 18 $\pm$ 3 | 0.61 $\pm$ 0.04
XGBoost | 5.61 $\pm$ 0.13 | 67 $\pm$ 12 | 0.26 $\pm$ 0.02
Skinny Trees |  0.31 $\pm$ 0.00 |  12 $\pm$ 1 |  0.87 $\pm$ 0.04
1000 | Random Forests | 2.90 $\pm$ 0.02 | 9 $\pm$ 0 | 0.94 $\pm$ 0.01
XGBoost | 1.45 $\pm$ 0.01 | 10 $\pm$ 0 | 0.91 $\pm$ 0.01
Skinny Trees |  0.26 $\pm$ 0.00 |  8 $\pm$ 0 |  1.00 $\pm$ 0.00

We evaluate our proposed method in a setting where the underlying data comes
from a sparse linear model with equi-spaced nonzeros with correlated design
matrix. Details are in Supplement Sec. S4. We experiment with a range of
training set sizes $N\in\\{100,200,1000\\}$, correlation strengths
$\sigma\in\\{0.5,0.7\\}$, and number of total features $p\in\\{256,512\\}$. We
evaluate the final performance averaged across 25 runs in terms of (i) test
MSE, (ii) number of features selected and (iii) support recovery via f1-score
between the true support and the recovered feature set. Skinny Trees
significantly outperforms both Random Forests and XGBoost in all three
measures across a range of settings. With Skinny Trees, we observe a 5-15 fold
improvement in mse performance and $9\%-65\%$ improvement in support recovery
metric (f1-score). The table shows that even if the features are correlated,
Skinny Trees successfully recovers the true support with high probability. We
also visualize this in Figure 3. Indices corresponding to those in the true
support are depicted in red. This confirms the usefulness of our end-to-end
feature selection approach.

##  7 Experiments

We study the performance of Skinny Trees on real-world datasets and compare
against popular baselines. We make the following comparisons: (i) Single
Skinny Tree versus sparse tree approaches. (ii) Skinny Trees versus dense soft
trees. (iii) Skinny Trees vs wrapped-based feature selection tree toolkits.
(iv) Skinny Trees against embedded feature selection toolkit from neural
networks. (v) Ablation study for dense-to-sparse learning for feature
selection.

Datasets. We use 14 open-source classification datasets (binary and
multiclass) from various domains with a range of number of features:
$20-100000$. Dataset details are in Table S3 in Supplement.

Tuning, Toolkits, and Details. For all the experiments, we tune the
hyperparameters using Optuna [2] with random search. The number of selected
features affects the AUC. Therefore, to treat all the methods in a fair
manner, we tune the hyperparameter that controls the sparsity level using
Optuna which optimizes the AUC across different $K$’s (feature selection
budget) e.g., $0.25p$ or $0.50p$ on a held-out validation set. Details are in
the Supplement.

###  7.1 Studying a single tree

We first study feature selection for a single tree on 4 classification tasks.
We study the performance of Skinny Tree (a single soft tree with group
$\ell_{0}-\ell_{2}$).

Competing Methods. We compare against:

  1. 1.

Classical Tree with TAO [13] with hyperplane splits that uses L1
regularization for node-level feature selection.

  2. 2.

Soft Tree with Group Lasso. Group Lasso is popularly used in high-dimensional
statistics and machine learning literature [70, 61]. We add
$\frac{\lambda_{1}}{\sqrt{m|\mathcal{I}|}}{\sum}_{k\in[p]}\left\lVert\bm{\mathcal{W}}_{k,:,:}\right\rVert_{2}$
as a regularization.

Table 2: Test AUC for TAO with $\ell_{1}$ regularization, single soft tree
with Group Lasso and Skinny Tree (a single soft tree with Group
$\ell_{0}-\ell_{2}$).

| Classical Tree | Soft Tree w/ | Skinny Tree
---|---|---|---
| TAO | Group Lasso |
Churn | 58.36 | 76.23 |  89.35$\pm$0.15
Satimage | 58.53 | 83.89 |  88.66$\pm$0.05
Texture | 58.90 | 93.83 |  98.42$\pm$0.01
Mice-protein | 57.13 | 87.88 |  99.19$\pm$0.00

Results. The numbers for classical-tree based TAO with $\ell_{1}$
regularization and soft tree with Group Lasso regularization are shown in
Table 2 for 50% sparsity budget. Results for Skinny Tree are also shown. We
see a huge gain in test AUC performance across all 4 datasets with Skinny Tree
in comparison with TAO and group lasso variant of a soft tree. This confirms
that in the context of feature selection at the ensemble level, a node-level
$\ell_{1}$ penalty is sub-optimal. Similarly, it also suggests that the joint
selection and shrinkage in Group Lasso can be less useful than Group
$\ell_{0}-\ell_{2}$, which decouples the selection and shrinkage.

###  7.2 Skinny Trees vs Dense Soft Trees

In this section, we compare our sparse trees with dense soft trees. For dense
soft trees, we use FASTEL [35] (an efficient state-of-the-art toolkit for
training soft tree ensembles). We present test AUC performance in Table 3.
Skinny Trees matches or outperforms dense soft trees in 10 datasets. Notably,
we observe a $30\%$ gain in test AUC on Madelon dataset with Skinny Trees. We
also observe $11\%$ improvements in test AUC on Cll and Gli datasets.
Additionally, sparse trees achieve $1.3\\!\times\\!-620\times$ feature
compression on 10 datasets. Note that in soft trees, feature compression has a
direct impact on model compression for reduced storage requirements and faster
inference. We observed up to $10\times$ faster inference times for Skinny
Trees compared to dense soft trees for compression rates of
$1.5\\!\times\\!-620\times$.

Table 3: Test AUC for Skinny Trees vs _dense_ Soft Trees. We also report
feature compression.

| Dense Trees | Skinny Trees | Compression
---|---|---|---
Churn | 91.15$\pm$0.09 |  93.20$\pm$0.08 | 1.8$\times$
Gisette |  99.81$\pm$0.003 |  99.81$\pm$0.002 | 1.5$\times$
Arcene | 89.57$\pm$0.11 |  90.80$\pm$0.30 | 2$\times$
Dorothea | 90.67$\pm$0.03 |  92.15$\pm$0.25 | 2.7$\times$
Madelon | 65.32$\pm$0.15 |  95.44$\pm$0.05 | 26$\times$
Smk |  84.10$\pm$0.16 | 79.29$\pm$0.22 | 253$\times$
Cll | 81.70$\pm$0.82 |  92.86$\pm$0.31 | 189$\times$
Gli | 88.65$\pm$0.90 |  99.80$\pm$0.07 | 619$\times$
Lung | 99.40$\pm$0.09 |  99.80$\pm$0.03 | 253$\times$
Tox | 99.19$\pm$0.04 |  99.74$\pm$0.02 | 189$\times$

Table 4: Test AUC (%) performance of Skinny Trees and feature-importance-based
toolkits for trees for $25\%$ feature budget ($K=0.25p$). Bold and italics
indicates best and runner-up models respectively.

Case |  Dataset |  Random Forests |  XGBoost |  LightGBM |  CatBoost |  Skinny Trees
---|---|---|---|---|---|---
$N<p$ |  Lung |  93.80$\pm$0.28 |  86.38$\pm$0.48 |  80.83$\pm$1.87 |  94.72$\pm$0.56 |  99.80$\pm$0.03
Tox |  94.52$\pm$0.14 |  97.10$\pm$0.09 |  95.94$\pm$0.54 |  95.95$\pm$0.14 |  99.74$\pm$0.02
Arcene |  74.80$\pm$0.36 |  76.36$\pm$0.16 |  76.92$\pm$0.36 |  76.64$\pm$0.22 |  80.80$\pm$0.30
Cll |  94.08$\pm$0.27 |  94.21$\pm$0.18 |  55.17$\pm$1.14 |  94.41$\pm$0.26 |  92.86$\pm$0.31
Smk |  77.78$\pm$0.20 |  76.88$\pm$0.40 |  67.29$\pm$0.91 |  78.44$\pm$0.41 |  79.29$\pm$0.22
Gli |  87.35$\pm$1.08 |  82.37$\pm$1.47 |  71.28$\pm$2.05 |  91.31$\pm$0.73 |  99.80$\pm$0.07
Dorothea |  89.71$\pm$0.12 |  89.09$\pm$0.09 |  88.14$\pm$0.18 |  88.50$\pm$0.27 |  90.87$\pm$0.02
$N>p$ |  Churn |  83.79$\pm$0.24 |  88.68$\pm$0.06 |  86.33$\pm$0.08 |  83.73$\pm$0.06 |  91.38$\pm$0.08
Satimage |  97.62$\pm$0.005 |  98.23$\pm$0.01 |  94.00$\pm$0.05 |  95.11$\pm$0.05 |  98.05$\pm$0.01
Texture |  99.60$\pm$0.003 |  99.94$\pm$0.001 |  96.14$\pm$0.03 |  94.90$\pm$0.07 |  99.97$\pm$0.002
Mice-protein |  99.30$\pm$0.01 |  99.77$\pm$0.01 |  89.59$\pm$0.22 |  95.03$\pm$0.07 |  99.59$\pm$0.02
Isolet |  99.17$\pm$0.002 |  99.86$\pm$0.002 |  97.62$\pm$0.003 |  99.89$\pm$0.001 |  99.94$\pm$0.01
Madelon |  94.11$\pm$0.02 |  94.65$\pm$0.01 |  86.46$\pm$0.08 |  96.41$\pm$0.01 |  94.14$\pm$0.09
Gisette |  98.99$\pm$0.004 |  99.64$\pm$0.004 |  98.09$\pm$0.50 |  99.57$\pm$0.01 |  99.81$\pm$0.002
|  Average |  91.75 |  91.65 |  84.56 |  91.76 |  94.72

###  7.3 Skinny Trees vs Classical Trees

We compare Skinny Trees against wrapper methods for feature selection from
classical trees (Random Forests, XGBoost, LightGBM, and CatBoost) on real-
world datasets. For Skinny Trees, we use the combined dense-to-sparse
scheduler. The tuning protocol and hyperparameters for all methods are
reported in the Supplement Sec. S5.3. The results are in Table 4. Skinny Trees
leads on 10 datasets. In contrast, others leads on 2 datasets. In terms of
test AUC, Skinny Trees outperforms LightGBM by $10.2\%$ (up to $37.7\%$),
XGBoost by $3.1\%$ (upto $17.4\%$), Random Forests by $3\%$ (up to $12.5\%$)
and CatBoost by $3\%$ (up to $8.5\%$). Overall, Skinny Trees provides a strong
alternative to existing wrapper-based methods.

Additional comparison with ControlBurn [45] is included in Supplement Sec. S7.
Skinny Trees also outperforms ControlBurn, achieving $2\%$ (up to $6\%$)
improvement in AUC.

###  7.4 Skinny Trees vs Neural Networks

In this paper, we pursue embedded feature selection methods for _tree
ensembles_. However, for completeness, we compare Skinny Trees against some
state-of-the-art embedded feature selection methods from neural networks,
namely LassoNet [43], AlgNet [18] and DFS [17]. Details are in Supplement Sec.
S5.4.

Results. We report AUC performance for $25\%$ feature budget in Table 5.
Skinny Trees leads across many datasets with statistical significance. In
terms of test AUC, Skinny Trees outperforms LassoNet by $4.3\%$ (up to
$24\%$), AlgNet by $25.9\%$ (up to $49\%$), and DFS by $2.6\%$ (up to
$11.4\%$).

Table 5: Test AUC (%) performance of Skinny Trees and embedded feature
selection methods from NN (LassoNet, AlgNet, DFS) for $25\%$ feature budget.

Dataset | LassoNet | AlgNet | DFS | Skinny Trees
---|---|---|---|---
Lung | 99.56$\pm$0.02 | 56.72$\pm$1.34 |  98.05$\pm$0.36 |  99.80$\pm$0.03
Tox | 99.63$\pm$0.03 | 51.01$\pm$0.70 |  99.13$\pm$0.24 |  99.74$\pm$0.02
Arcene | 66.26$\pm$0.29 | 51.00$\pm$1.77 | 69.37$\pm$0.59 |  80.80$\pm$0.30
Cll |  95.04$\pm$0.24 | 64.96$\pm$2.54 | 92.85$\pm$0.32 |  92.86$\pm$0.31
Smk |  85.56$\pm$0.31 | 53.92$\pm$2.16 |  79.62$\pm$0.29 | 79.29$\pm$0.22
Gli |  97.78$\pm$0.56 | 61.37$\pm$4.27 | 92.25$\pm$0.64 |  99.80$\pm$0.07
Dorothea | out of mem. | 81.74$\pm$0.91 | 85.18 |  90.87$\pm$0.02
Churn | 67.34$\pm$1.58 | 70.10$\pm$0.91 |  85.70$\pm$0.52 |  91.38$\pm$0.08
Satimage | 94.73$\pm$0.19 | 95.30$\pm$0.20 |  97.39$\pm$0.04 |  98.05$\pm$0.01
Texture | 98.02$\pm$0.40 | 76.24$\pm$1.94 |  99.63$\pm$0.04 |  99.97$\pm$0.002
Mice-protein | 94.90$\pm$0.26 | 89.07$\pm$0.59 |  99.04$\pm$0.03 |  99.59$\pm$0.02
Isolet | 99.64$\pm$0.01 | 70.21$\pm$2.92 |  99.92$\pm$0.00 |  99.94$\pm$0.01
Madelon | 81.15$\pm$2.53 | 68.55$\pm$1.42 |  92.73$\pm$0.45 |  94.14$\pm$0.09
Gisette | 99.81$\pm$0.002 | 73.49$\pm$1.54 | $\textit{99.72}^{*}$ |  99.81$\pm$0.002
Average | 90.43** | 68.83 | 92.18 | 94.72
∗DFS is very time-consuming to run, we report the test AUC for best
trial (based on validation AUC) during tuning on Gisette and Dorothea.
∗∗Adjusted Average: $\frac{90.72}{(94.72*14-90.87)/13}*94.72=90.43$.

###  7.5 Dense-to-Sparse Learning

Smk | Dorothea
---|---
|
Figure 4: Performance without/with Dense-to-sparse learning for different
feature selection budgets.

We perform an ablation study in which we compare the predictive performance
achieved with dense-to-sparse learning (DSL) over a range of feature selection
budgets. Tuning details are in the Supplement Sec. S5.5. The results are
reported in Figure 4. Interestingly, we significantly improve in test AUC
across a range of feature selection budgets with dense-to-sparse learning over
fixed regularization tuning.

###  7.6 Discussion on training times.

Skinny Trees is very competitive in terms of training times in comparison to
existing toolkits. We compared timings on a single Tesla V100 GPU. For
example, on dorothea dataset, Skinny Trees trained in under 3 minutes for
optimal hyperparameter setting. XGBoost took 10 minutes. In contrast, DFS took
45 hours.

##  8 Conclusion

We introduce an end-to-end optimization approach for joint feature selection
and tree ensemble learning. Our approach is based on differentiable trees with
group $\ell_{0}-\ell_{2}$ regularization. We use a simple but effective
proximal mini-batch gradient descent algorithm and present convergence
guarantees. Interestingly, a dense-to-sparse regularization scheduling
approach can lead to more expressive and sparser tree ensembles. We
demonstrate on various datasets that our toolkit Skinny Trees can
significantly improve feature selection over several state-of-the-art wrapper-
based feature selection methods in trees and embedded feature selection
methods in neural networks.

##  9 Acknowledgments

This research was supported in part, by grants from the Office of Naval
Research (N000142112841), and Liberty Mutual Insurance. The authors
acknowledge the MIT SuperCloud [58] for providing HPC resources that have
contributed to the research reported within this paper.

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, et al.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.
  * [2] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, et al.  Optuna: A next-generation hyperparameter optimization framework.  In Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
  * [3] Genevera I. Allen.  Automatic feature selection via weighted kernels and regularization.  Journal of Computational and Graphical Statistics, 22(2):284–299, 04 2013.
  * [4] Hedy Attouch, Jérôme Bolte, and Benar Fux Svaiter.  Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods.  Mathematical Programming, 137(1-2):91–129, 2013.
  * [5] R. Battiti.  Using mutual information for selecting features in supervised neural net learning.  IEEE Transactions on Neural Networks, 5(4):537–550, 1994.
  * [6] D. Bertsimas and J. Dunn.  Machine Learning Under a Modern Optimization Lens.  Dynamic Ideas LLC, 2019.
  * [7] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  The Annals of Statistics, 44(2):813–852, April 2016.
  * [8] Jérôme Bolte, Aris Daniilidis, and Adrian Lewis.  The łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems.  SIAM Journal on Optimization, 17(4):1205–1223, 2007.
  * [9] A.-L. Boulesteix, A. Bender, J. Lorenzo Bermejo, and C. Strobl.  Random forest gini importance favours SNPs with large minor allele frequency: impact, sources and recommendations.  Briefings in Bioinformatics, 13(3):292–304, September 2011.
  * [10] L. Breiman, J. Friedman, C.J. Stone, and R.A. Olshen.  Classification and Regression Trees.  Taylor & Francis, 1984.
  * [11] Leo Breiman.  Random forests.  Machine Learning, 45(1):5–32, 2001.
  * [12] Jie Cai, Jiawei Luo, Shulin Wang, and Sheng Yang.  Feature selection in machine learning: A new perspective.  Neurocomputing, 300:70–79, 2018.
  * [13] Miguel A. Carreira-Perpinan and Pooya Tavallali.  Alternating optimization of decision trees, with application to learning sparse oblique trees.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
  * [14] Girish Chandrashekar and Ferat Sahin.  A survey on feature selection methods.  Computers & Electrical Engineering, 40(1):16–28, 2014.  40th-year commemorative issue.
  * [15] Jianbo Chen, Mitchell Stern, Martin J Wainwright, and Michael I Jordan.  Kernel feature selection via conditional covariance minimization.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [16] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 785–794, New York, NY, USA, 2016. Association for Computing Machinery.
  * [17] Yao Chen, Qingyi Gao, Faming Liang, et al.  Nonlinear variable selection via deep neural networks.  Journal of Computational and Graphical Statistics, 30(2):484–492, 2021.
  * [18] V. Dinh and L. S. T. Ho.  Consistent feature selection for analytic deep neural networks.  In NIPS’20, NIPS’20, Red Hook, NY, USA, 2020. Curran Associates Inc.
  * [19] L. P. D. van den Dries.  SEMIALGEBRAIC SETS, page 31–42.  London Mathematical Society Lecture Note Series. Cambridge University Press, 1998.
  * [20] Chandra Erdman and Nancy Bates.  The low response score (LRS).  Public Opinion Quarterly, 81(1):144–156, December 2016.
  * [21] Pablo A. Estevez, Michel Tesmer, Claudio A. Perez, and Jacek M. Zurada.  Normalized mutual information feature selection.  IEEE Transactions on Neural Networks, 20(2):189–201, 2009.
  * [22] Rong-En Fan, Kai-Wei Chang, Cho-Jui Hsieh, Xiang-Rui Wang, and Chih-Jen Lin.  Liblinear: A library for large linear classification.  J. Mach. Learn. Res., 9:1871–1874, jun 2008.
  * [23] Jean Feng and Noah Simon.  Sparse-input neural networks for high-dimensional nonparametric regression and classification.  arXiv: Methodology, 2017.
  * [24] Nicholas Frosst and Geoffrey Hinton.  Distilling a neural network into a soft decision tree.  ArXiv, abs/1711.09784, 2017.
  * [25] Isabelle Guyon, Steve Gunn, Asa Ben-Hur, and Gideon Dror.  Result analysis of the nips 2003 feature selection challenge.  In L. Saul, Y. Weiss, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 17. MIT Press, 2004.
  * [26] T. J. Hastie, R. J. Tibshirani, and J. Friedman.  The Elements of Statistical Learning.  Springer, 2 edition, 2009.
  * [27] H. Hazimeh, N. Ponomareva, P. Mol, et al.  The tree ensemble layer: Differentiability meets conditional computation.  In Hal Daumé III and Aarti Singh, editors, ICML’20, volume 119 of Proceedings of Machine Learning Research, pages 4138–4148. PMLR, 13–18 Jul 2020.
  * [28] Hussein Hazimeh and Rahul Mazumder.  Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms.  Oper. Res., 68(5):1517–1537, September 2020.
  * [29] Hussein Hazimeh, Rahul Mazumder, and Peter Radchenko.  Grouped variable selection with discrete optimization: Computational and statistical perspectives.  arXiv preprint arXiv:2104.07084 (Annals of Statistics, to appear), 2022.
  * [30] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [31] Thomas M. Hehn, Julian F. P. Kooij, and Fred A. Hamprecht.  End-to-end learning of decision trees and forests.  International Journal of Computer Vision, 128:997–1011, 2019.
  * [32] Laurent Hyafil and Ronald L. Rivest.  Constructing optimal binary decision trees is np-complete.  Information Processing Letters, 5(1):15–17, 1976.
  * [33] Shibal Ibrahim, Gabriel Afriat, Kayhan Behdin, and Rahul Mazumder.  GRAND-SLAMIN’ interpretable additive modeling with structural constraints.  In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
  * [34] Shibal Ibrahim, Wenyu Chen, Hussein Hazimeh, Natalia Ponomareva, Zhe Zhao, and Rahul Mazumder.  Comet: Learning cardinality constrained mixture of experts with trees and local search.  KDD ’23, page 832–844, New York, NY, USA, 2023. Association for Computing Machinery.
  * [35] Shibal Ibrahim, Hussein Hazimeh, and Rahul Mazumder.  Flexible modeling and multitask learning using differentiable tree ensembles.  In KDD’22, KDD ’22, page 666–675, New York, NY, USA, 2022. Association for Computing Machinery.
  * [36] Shibal Ibrahim, Rahul Mazumder, Peter Radchenko, and Emanuel Ben-David.  Predicting census survey response rates via interpretable nonparametric additive models with structured interactions, 2021.
  * [37] Xiaojie Jin, Xiaotong Yuan, Jiashi Feng, and Shuicheng Yan.  Training skinny deep neural networks with iterative hard thresholding methods, 2016.
  * [38] Michael I. Jordan and Robert A. Jacobs.  Hierarchical mixtures of experts and the em algorithm.  Neural Comput., 6(2):181–214, mar 1994.
  * [39] Guolin Ke, Qi Meng, Thomas Finley, et al.  Lightgbm: A highly efficient gradient boosting decision tree.  In I. Guyon, U. V. Luxburg, S. Bengio, et al., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [40] Ron Kohavi and George H. John.  Wrappers for feature subset selection.  Artificial Intelligence, 97(1):273–324, 1997.
  * [41] Peter Kontschieder, Madalina Fiterau, Antonio Criminisi, et al.  Deep neural decision forests.  In 2015 IEEE International Conference on Computer Vision (ICCV), pages 1467–1475, 2015.
  * [42] Krzysztof Kurdyka.  On gradients of functions definable in o-minimal structures.  Annales de l’institut Fourier, 48(3):769–783, 1998.
  * [43] Ismael Lemhadri, Feng Ruan, and Rob Tibshirani.  Lassonet: Neural networks with feature sparsity.  In Arindam Banerjee and Kenji Fukumizu, editors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 10–18. PMLR, 13–15 Apr 2021.
  * [44] Jundong Li, Kewei Cheng, Suhang Wang, Fred Morstatter, Robert P. Trevino, Jiliang Tang, and Huan Liu.  Feature selection: A data perspective.  ACM Comput. Surv., 50(6), dec 2017.
  * [45] Brian Liu, Miaolan Xie, and Madeleine Udell.  Controlburn: Feature selection by sparse forests.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD ’21, page 1045–1054, New York, NY, USA, 2021\. Association for Computing Machinery.
  * [46] Christos Louizos, Max Welling, and Diederik P. Kingma.  Learning sparse neural networks through l0 regularization.  In International Conference on Learning Representations, 2018.
  * [47] Rahul Mazumder, Peter Radchenko, and Antoine Dedieu.  Subset selection with shrinkage: Sparse linear modeling when the snr is low.  arXiv preprint arXiv:1708.03288 (Operations Research, to appear), 2017.
  * [48] Alan J. Miller.  Subset selection in regression.  Chapman & Hall/CRC, 1990.
  * [49] Fan Min, Qinghua Hu, and William Zhu.  Feature selection with test cost constraint.  International Journal of Approximate Reasoning, 55(1, Part 2):167–179, 2014.  Special issue on Decision-Theoretic Rough Sets.
  * [50] Sreerama K. Murthy, Simon Kasif, and Steven Salzberg.  A system for induction of oblique decision trees.  J. Artif. Int. Res., 2(1):1–32, aug 1994.
  * [51] Randal S. Olson, William La Cava, Patryk Orzechowski, Ryan J. Urbanowicz, and Jason H. Moore.  Pmlb: a large benchmark suite for machine learning evaluation and comparison.  BioData Mining, 10(1):36, 2017.
  * [52] Hanchuan Peng, Fuhui Long, and C. Ding.  Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 27(8):1226–1238, 2005.
  * [53] Alexandra Peste, Eugenia Iofinova, Adrian Vladu, and Dan Alistarh.  AC/DC: Alternating compressed/decompressed training of deep neural networks.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [54] Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, et al.  Catboost: Unbiased boosting with categorical features.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, page 6639–6649, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * [55] J. Ross Quinlan.  C4.5: Programs for Machine Learning.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1993.
  * [56] P. Ravikumar, J. Lafferty, H. Liu, et al.  Sparse additive models.  Journal of the Royal Statistical Society, B., 71:1009–1030, 2009\.
  * [57] Juha Reunanen.  Overfitting in making comparisons between variable selection methods.  J. Mach. Learn. Res., 3(null):1371–1382, mar 2003.
  * [58] Albert Reuther, Jeremy Kepner, Chansup Byun, Siddharth Samsi, William Arcand, David Bestor, Bill Bergeron, Vijay Gadepally, Michael Houle, Matthew Hubbell, Michael Jones, Anna Klein, Lauren Milechin, Julia Mullen, Andrew Prout, Antonio Rosa, Charles Yee, and Peter Michaleas.  Interactive supercomputing on 40,000 cores for machine learning and data analysis.  In 2018 IEEE High Performance extreme Computing Conference (HPEC), pages 1–6. IEEE, 2018.
  * [59] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  ”why should i trust you?”: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 1135–1144, New York, NY, USA, 2016. Association for Computing Machinery.
  * [60] Debaditya Roy, K. Sri Rama Murty, and C. Krishna Mohan.  Feature selection using deep neural networks.  In 2015 International Joint Conference on Neural Networks (IJCNN), pages 1–6, 2015.
  * [61] Simone Scardapane, Danilo Comminiello, Amir Hussain, and Aurelio Uncini.  Group sparse regularization for deep neural networks.  Neurocomputing, 241:81–89, 2017.
  * [62] Le Song, Alex Smola, Arthur Gretton, Justin Bedo, and Karsten Borgwardt.  Feature selection via dependence maximization, 2012.
  * [63] Le Song, Alex Smola, Arthur Gretton, Karsten M. Borgwardt, and Justin Bedo.  Supervised feature selection via dependence estimation.  In Proceedings of the 24th International Conference on Machine Learning, ICML ’07, page 823–830, New York, NY, USA, 2007. Association for Computing Machinery.
  * [64] Gary Stein, Bing Chen, Annie S. Wu, and Kien A. Hua.  Decision tree classifier for network intrusion detection with ga-based feature selection.  In Proceedings of the 43rd Annual Southeast Regional Conference \- Volume 2, ACM-SE 43, page 136–141, New York, NY, USA, 2005. Association for Computing Machinery.
  * [65] Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis, and Torsten Hothorn.  Bias in random forest variable importance measures: Illustrations, sources and a solution.  BMC Bioinformatics, 8(1):25, 2007.
  * [66] Ryutaro Tanno, Kai Arulkumaran, Daniel C. Alexander, Antonio Criminisi, and Aditya V. Nori.  Adaptive neural trees.  ArXiv, abs/1807.06699, 2019.
  * [67] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 2022/05/03 1996.
  * [68] Julia D. Wulfkuhle, Lance A. Liotta, and Emanuel F. Petricoin.  Proteomic applications for the early detection of cancer.  Nature Reviews Cancer, 3(4):267–275, 2003.
  * [69] Yutaro Yamada, Ofir Lindenbaum, Sahand Negahban, and Yuval Kluger.  Feature selection using stochastic gates.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 10648–10659. PMLR, 13–18 Jul 2020.
  * [70] M. Yuan and Y. Lin.  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society, Series B, 68:49–67, 2006\.
  * [71] Shichao Zhang.  Cost-sensitive classification with respect to waiting cost.  Know.-Based Syst., 23(5):369–378, jul 2010.
  * [72] Tuo Zhao and Han Liu.  Sparse additive machine.  In Neil D. Lawrence and Mark Girolami, editors, Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, volume 22 of Proceedings of Machine Learning Research, pages 1435–1443, La Palma, Canary Islands, 21–23 Apr 2012. PMLR.
  * [73] Arman Zharmagambetov, Suryabhan Singh Hada, and Miguel Á. Carreira-Perpiñán.  An experimental comparison of old and new decision tree algorithms.  CoRR, abs/1911.03054, 2019.
  * [74] Arman Serikuly Zharmagambetov and Miguel Á. Carreira-Perpiñán.  Smaller, more accurate regression forests using tree alternating optimization.  In ICML, 2020.
  * [75] Zhengze Zhou and Giles Hooker.  Unbiased measurement of feature importance in tree-based methods.  ACM Transactions on Knowledge Discovery from Data (TKDD), 15:1 – 21, 2021.
  * [76] Zexuan Zhu, Yew-Soon Ong, and Manoranjan Dash.  Wrapper–filter feature selection algorithm using a memetic framework.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 37(1):70–76, 2007.

## Supplementary Material for “End-to-end Feature Selection Approach for
Learning Skinny Trees”

##  Appendix S1 Notations and Acronyms

Notation Table S1 lists the notation used throughout the paper.

Table S1: List of notation used.

Notation | Space or Type |  Explanation
---|---|---
$[n]$ | Set |  The set of integers $\\{1,2,....,n\\}$.
$\bm{1}_{m}$ | $\mathbb{R}^{m}$ |  Vector with all coordinates equal to 1.
$\bm{U}$ | $\mathbb{R}^{m,n}$ |  Matrix with elements $((U_{ij}))$
$\bm{u}\cdot\bm{v}$ | $\mathbb{R}$ |  A dot product between two vectors $\bm{u},\bm{v}$.
$\bm{U}\cdot\bm{v}$ | $\mathbb{R}^{n}$ |  A dot product between a matrix $\bm{U}\in R^{m,n}$ and a vector $\bm{v}\in\mathbb{R}^{m}$ is denoted as $\bm{U}\cdot\bm{v}=\bm{U}^{T}v\in\mathbb{R}^{n}$.
$\mathcal{X}$ | $\mathbb{R}^{p}$ |  Input feature space.
$\mathcal{Y}$ | $\mathbb{R}^{c}$ |  Output (label) space.
$m$ | $\mathbb{Z}_{>0}$ |  Number of trees in Skinny Trees.
$\bm{f}(\bm{x})$ | Function |  The output of Skinny Trees, a function that takes an input sample and returns a logit which corresponds to the sum of all the trees in the ensemble. Formally, $\bm{f}:\mathcal{X}\to\mathbb{R}^{c}.$
$\bm{f}^{j}(\bm{x})$ | Function |  A single perfect binary tree which takes an input sample and returns a logit, i.e., $\bm{f}^{j}:\mathcal{X}\to\mathbb{R}^{c}$.
$d$ | $\mathbb{Z}_{>0}$ |  The depth of tree $\bm{f}^{j}$.
$\mathcal{I}^{j}$ | Set |  The set of internal (split) nodes in $\bm{f}^{j}$.
$\mathcal{I}$ | Set |  The set of internal (split) supernodes in $\bm{f}$.
$\mathcal{L}^{j}$ | Set |  The set of leaf nodes in $\bm{f}^{j}$.
$\mathcal{A}(i)$ | Set |  The set of ancestors of node $i$.
$\\{x\to i\\}$ | Event |  The event that sample $x\in\mathbb{R}^{p}$ reaches node $i$.
$\bm{w}_{i}$ | $\mathbb{R}^{p}$ |  Weight vector of internal node $i$ (trainable). Defines the hyperplane split used in sample routing.
$\bm{W}_{i}$ | $\mathbb{R}^{p,m}$ |  Matrix of all weights in the internal supernode $i$ of the ensemble in the tensor-formulation.
$\bm{\mathcal{W}}$ | $\mathbb{R}^{p,m,|\mathcal{I}|}$ |  Tensor of all weights across all internal supernodes in the ensemble.
$\bm{\mathcal{W}}_{k,:,:}$ | $\mathbb{R}^{m,|\mathcal{I}|}$ |  Matrix of all weights for $k$-th feature/covariate across all internal supernodes in the ensemble.
$S$ | Function |  Activation function $\mathbb{R}\to[0,1]$
$S(\bm{w}_{i}\cdot\bm{x})$ | $[0,1]$ |  Probability (proportion) that internal node $i$ routes $x$ to the left.
$S^{\prime}(v)$ | Function |  The derivative of $S(v)$
$[l\swarrow i]$ | Event |  The event that leaf $l$ belongs to the left subtree of node $i\in\mathcal{I}$.
$[l\searrow i]$ | Event |  The event that leaf $l$ belongs to the right subtree of node $i\in\mathcal{I}$.
$\bm{o}_{l}$ | $\mathbb{R}^{c}$ |  Leaf $l$’s weight vector (trainable).
$\bm{O}_{l}$ | $\mathbb{R}^{m,c}$ |  Matrix of weights in superleaf $l$.
$\bm{\mathcal{O}}$ | $\mathbb{R}^{m,c,|\mathcal{L}|}$ |  Tensor of weights across the superleaves in the ensemble.
$L$ | Function |  Loss function for training (e.g., cross-entropy).
$\mathcal{Q}$ | Set |  Unknown (learnable) subset of features of size at most $K$.
$z_{k}$ | $\\{0,1\\}$ |  Binary variable controls whether $k$-th feature is on or off in Problem (2).
$\bm{z}$ | $\\{0,1\\}^{p}$ |  Binary vector controlling which features are on or off in Problem (2).
$\lambda_{0}$ | ${\mathbb{R}}_{\geq 0}$ |  Non-negative $\ell_{0}$ regularization parameter controlling the number of features selected in Problem (3)
$\lambda_{2}$ | ${\mathbb{R}}_{\geq 0}$ |  Non-negative $\ell_{2}$ regularization parameter controlling the shrinkage in Problem (3)
$\lambda_{1}$ | ${\mathbb{R}}_{\geq 0}$ |  Non-negative $\ell_{0}$ regularization parameter controlling the number of features selected and shrinkage in Problem (S19)
$\gamma$ | ${\mathbb{R}}_{\geq 0}$ |  Non-negative scaling parameter for the exponential ramp-up of $\ell_{0}$-penalty in dense-to-sparse learning.
$s$ | ${\mathbb{R}}_{\geq 0}$ |  Non-negative temperature parameter for controlling the ramping rate of the $\ell_{0}$-penalty in dense-to-sparse learning.
$\eta$ | ${\mathbb{R}}_{\geq 0}$ |  Learning rate parameter for proximal mini-batch gradient descent.
$b$ | ${\mathbb{Z}}_{>0}$ |  Batch size scaling parameter for large-batch training.

Acronyms Table S2 lists the acronyms used throughout the paper.

Table S2: List of Acronyms used. Terms | Acronyms
---|---
Gradient Descent | GD
Dense-to-sparse learning | DSL

##  Appendix S2 Background: Differentiable (a.k.a. soft) decision tree

A soft tree is a variant of a classical decision tree that performs soft
routing, i.e., a sample is fractionally routed to all leaves. It was proposed
by [38], and further developed in [41, 27] for end-to-end optimization. Soft
routing makes soft trees differentiable, so learning can be done using
gradient-based methods.

Let us fix some $j\in[m]$ and consider a single tree $\bm{f}^{j}$, which takes
an input sample and returns an output vector (logit), i.e.,
$\bm{f}^{j}:X\in\mathbb{R}^{p}\rightarrow\mathbb{R}^{c}$. Moreover, we assume
that $\bm{f}^{j}$ is a perfect binary tree with depth $d$. Let
$\mathcal{I}^{j}$ and $\mathcal{L}^{j}$ denote the sets of internal (split)
nodes and the leaves of the tree, respectively. For any node
$i\in\mathcal{I}^{j}\cup\mathcal{L}^{j}$, we define $A^{j}(i)$ as its set of
ancestors and use the notation $\bm{x}\rightarrow i$ for the event that a
sample $\bm{x}\in\mathbb{R}^{p}$ reaches $i$. See Table S1 for detailed
notation summary.

Routing Following existing work [41, 31, 27], we present routing in soft trees
with a probabilistic model. Although the sample routing is formulated with a
probabilistic model, the final prediction of the tree $\bm{f}$ is a
deterministic function as it assumes an expectation over the leaf predictions.
According to this probabilistic model, internal (split) nodes in a soft tree
perform soft routing, where a sample is routed left and right with different
probabilities. Classical decision trees are modeled with either axis-aligned
splits [10, 55] or hyperplane (a.k.a. oblique) splits [50]. Soft trees are
based on hyperplane splits, where the routing decisions rely on a linear
combination of the features. Particularly, each internal node
$i\in\mathcal{I}^{j}$ is associated with a trainable weight vector
$\bm{w}_{i}^{j}\in\mathbb{R}^{p}$ that defines the node’s hyperplane split.

Given a sample $\bm{x}\in\mathbb{R}^{p}$, the probability that internal node
$i$ routes $\bm{x}$ to the left is defined by $S(\bm{w}_{i}^{j}\cdot\bm{x})$.
Now we discuss how to model the probability that $\bm{x}$ reaches a certain
leaf $l$. Let $[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]$ (resp.
[$i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l$]) denote the event that leaf
$l$ belongs to the left (resp. right) subtree of node $i\in\mathcal{I}^{j}$.
The probability that $\bm{x}$ reaches $l$ is given by $P^{j}(\\{x\rightarrow
l\\})=\prod_{i\in A(l)}r_{i,l}^{j}(\bm{x})$, where $r_{i,l}^{j}(\bm{x})$ is
the probability of node $i$ routing $\bm{x}$ towards the subtree containing
leaf $l$, i.e.,
$r_{i,l}^{j}(x):=S(\bm{w}_{i}^{j}\cdot\bm{x})^{1[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]}\odot(1-S(\bm{w}_{i}^{j}\cdot\bm{x}))^{1[i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l]}$. Let $S:R\rightarrow[0,1]$ be
an activation function. Popular choices for $S$ include logistic function [38,
41, 24, 66, 31] and Smooth-step function (for hard routing) [27]. Next, we
define how the root-to-leaf probabilities can be used to make the final
prediction of the tree.

Prediction As with classical decision trees, we assume that each leaf stores a
learnable weight vector $\bm{o}_{l}^{j}\in R^{c}$. For a sample
$\bm{x}\in\mathbb{R}^{p}$, prediction of the tree is defined as an expectation
over the leaf outputs, i.e.,
$\bm{f}^{j}(\bm{x})=\sum_{l\in\mathcal{L}^{j}}P^{j}(\\{\bm{x}\rightarrow
l\\})\bm{o}_{l}^{j}$.

###  S2.1 Tree Ensemble Tensor Formulation

[35] proposed a tensor formulation for modeling tree ensembles more
efficiently, which can lead to faster training times than classical
formulations [41, 27]: $\sim 10\times$ on CPUs and $\sim 20\times$ on GPUs. We
use a similar tensor formulation, as we discuss below. The internal nodes in
the trees across the ensemble are jointly modeled as a “supernode”. In
particular, an internal node $i\in\mathcal{I}^{j}$ at depth $d$ in all trees
can be condensed together into a supernode $i\in\mathcal{I}$. Let
$\bm{W}_{i}\in\mathbb{R}^{p,m}$ be learnable weight matrix, where each $j$-th
column of the weight matrix contains the learnable weight vector
$\bm{w}_{i}^{j}$ of the original $j$-th tree in the ensemble. Similarly, the
leaves in the trees across the ensemble are jointly modeled as a superleaf.
Let $\bm{O}_{l}\in\mathbb{R}^{m,c}$ be the learnable weight matrix to store
the leaf nodes, where each $j$-th row contains the learnable weight vector
$\bm{o}_{l}^{j}$ in the original $j$-th tree in the ensemble. The prediction
of the tree ensemble is
$\bm{f}(\bm{x})=(\sum_{l\in\mathcal{L}}\bm{O}_{l}\odot\prod_{i\in
A(l)}\bm{R}_{i,l})\cdot\bm{1}_{m},$ where $\odot$ denotes the element-wise
product,
$\bm{R}_{i,l}=S(\bm{W}_{i}\cdot\bm{x})^{1[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]}\odot(1-S(\bm{W}_{i}\cdot\bm{x}))^{1[i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l]}\in\mathbb{R}^{m,1}$ and the
activation function $S$ is applied element-wise. $\bm{1}_{m}\in\mathbb{R}^{m}$
is a vector of ones that combines the predictions of trees in the ensemble. We
denote all the hyperplane parameters across all the supernodes of the tree
ensemble as a tensor $\bm{\mathcal{W}}\in R^{p,m,|\mathcal{I}|}$ and the
parameters across all the superleaves as a tensor $\bm{\mathcal{O}}\in
R^{m,c,|\mathcal{L}|}$.

##  Appendix S3 Proof of Theorem 1

Overview and Preliminaries. Let us denote the training loss corresponding to
the sample $n$ as

$\Phi_{n}(\bm{\mathcal{O}},\bm{\mathcal{W}})=L(y_{n},\bm{f}({\mathbf{x}}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}})).$
(S1)

We also use the notation ${\mathbb{X}}$ to denote the set of all decision
variables in the model, ${\mathbb{X}}:=(\bm{\mathcal{W}},\bm{\mathcal{O}})$.
For $j\in[m]$ and $i\in\mathcal{I}^{j}$ and $t\in[p]$, we let $w^{j}_{i,t}$ be
the $t$-th coordinate of $\bm{w}_{i}^{j}$.

In this proof, we follow the general steps outlined below:

  1. 1.

First, we show that $\Phi_{n}$ is $\mathcal{M}$-smooth for some
$\mathcal{M}>0$ only depending on the data and the constants appearing in
Assumptions (A1), (A2) and (A3). That is, there exists $\mathcal{M}>0$ such
that for any two
${\mathbb{X}}_{1}=(\bm{\mathcal{W}}_{1},\bm{\mathcal{O}}_{1}),{\mathbb{X}}_{2}=(\bm{\mathcal{W}}_{2},\bm{\mathcal{O}}_{2})$
with $\|\bm{\mathcal{O}}_{1}\|_{2},\|\bm{\mathcal{O}}_{2}\|_{2}\leq B$, we
have

$\|\nabla\Phi_{n}(\bm{\mathcal{O}}_{1},\bm{\mathcal{W}}_{1})-\nabla\Phi_{n}(\bm{\mathcal{O}}_{2},\bm{\mathcal{W}}_{2})\|_{2}\leq\mathcal{M}\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2}.$
(S2)

This will prove the descent property of the algorithm.

  2. 2.

Next, we show that as long as $\lambda_{2}>0$, the sequence of solutions
generated by the algorithm is bounded.

  3. 3.

Finally, we show that $\Phi_{n}$ is semi-algebraic [19, 4, Chapter 2] and
therefore satisfies the Kurdyka-Łojasiewicz (KL) property [8, 42]. This will
complete the proof of convergence.

Before continuing with the proof, we derive some results that will be useful.
For notational convenience, we drop the sample index $n$ as our results will
be true for all samples.

First, for $j\in[m]$, $i\in\mathcal{I}^{j}$ and $l\in\mathcal{L}^{j}$

$\displaystyle\frac{\partial
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})}{\partial\bm{w}_{i}^{j}}$
$\displaystyle=\sum_{l\in\mathcal{L}^{j}}o_{l}^{j}\frac{\partial
P^{j}(\\{{\mathbf{x}}\rightarrow l\\})}{\partial\bm{w}_{i}^{j}}$ (S3)
$\displaystyle\frac{\partial
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})}{\partial o_{l}^{j}}$
$\displaystyle=P^{j}(\\{{\mathbf{x}}\rightarrow l\\}).$

Thus,

$\displaystyle\frac{\partial\Phi}{\partial\bm{w}_{i}^{j}}$
$\displaystyle=L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))\frac{\partial
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})}{\partial\bm{w}_{i}^{j}}$
$\displaystyle=L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))\sum_{l\in\mathcal{L}^{j}}o_{l}^{j}\frac{\partial
P^{j}(\\{{\mathbf{x}}\rightarrow l\\})}{\partial\bm{w}_{i}^{j}}.$ (S4)

If $i\notin\mathcal{A}(l)$, then $\frac{\partial
P^{j}(\\{{\mathbf{x}}\rightarrow l\\})}{\partial\bm{w}_{i}^{j}}=\bm{0}$.
Otherwise,

$\displaystyle\frac{\partial P^{j}(\\{{\mathbf{x}}\rightarrow
l\\})}{\partial\bm{w}_{i}^{j}}=\prod_{\begin{subarray}{c}k\in\mathcal{A}(l)\\\
k\neq i\end{subarray}}r_{k,l}^{j}({\mathbf{x}})\frac{\partial
r_{i,l}^{j}({\mathbf{x}})}{\partial\bm{w}_{i}^{j}}.$ (S5)

Moreover, by the definition of $r_{i,l}^{j}$, we have

$\frac{\partial
r_{i,l}^{j}({\mathbf{x}})}{\partial\bm{w}_{i}^{j}}=\begin{cases}S^{\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i\\\
-S^{\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i.\end{cases}$ (S6)

In addition, for $j\in[m]$ and $l\in\mathcal{L}^{j}$,

$\displaystyle\frac{\partial\Phi}{\partial o_{l}^{j}}$
$\displaystyle=L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))\frac{\partial
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})}{\partial o_{l}^{j}}$
$\displaystyle=L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))P^{j}(\\{{\mathbf{x}}\rightarrow
l\\}).$ (S7)

We define

$\displaystyle\phi_{i,t}^{j}({\mathbb{X}})$
$\displaystyle=\frac{\partial\Phi}{\partial{w}_{i,t}^{j}},$ (S8)
$\displaystyle\phi_{l}^{j}({\mathbb{X}})$
$\displaystyle=\frac{\partial\Phi}{\partial{o}_{l}^{j}}.$

Next, we state a few technical lemma that will be useful in our proof.

######  Lemma 1.

Define
$E_{i,+}^{j},E_{i,-}^{j}\subseteq\mathbb{R}^{p,m,|\mathcal{I}|}\times\mathbb{R}^{m,|\mathcal{L}|}$
for $i\in\mathcal{I}^{j},j\in[m]$:

$\displaystyle
E_{i,+}^{j}=\\{\bm{w}_{i}^{j}\cdot{\mathbf{x}}=\theta\\},~{}~{}E_{i,-}^{j}=\\{\bm{w}_{i}^{j}\cdot{\mathbf{x}}=-\theta\\}$
(S9)

and

$\mathcal{D}=\left(\left[\mathbb{R}^{p,m,|\mathcal{I}|}\times\mathbb{R}^{m,|\mathcal{L}|}\right]\bigcap\\{\|\bm{\mathcal{O}}\|_{2}\leq
B\\}\right)\setminus\bigcup_{\begin{subarray}{c}j\in[m]\\\
i\in\mathcal{I}^{j}\end{subarray}}\left\\{E_{i,+}^{j}\cup
E_{i,-}^{j}\right\\}$

where $B$ is defined in Assumption (A3). If ${\mathbb{X}}\in\mathcal{D}$, then
$\Phi({\mathbb{X}})$ is infinitely differentiable. Moreover, $\mathcal{D}$ can
be partitioned into finitely many subsets.

###### Proof.

Note that $\Phi({\mathbb{X}})$ may not have infinitely many derivatives only
if $S(\bm{w}_{i}^{j}\cdot{\mathbf{x}})$ is not smooth for some
$j\in[m],i\in\mathcal{I}^{j}$. By the construction of $S(\cdot)$ from
Assumption (A1), the activation $S(x)$ function is not infinitely
differentiable only for $x=\pm\theta$. As a result, if
${\mathbb{X}}\in\mathcal{D}$, all activation functions are smooth and
therefore $\Phi({\mathbb{X}})$ is infinitely differentiable.
Moreover, note that each $E_{i,\pm}^{j}$ is an affine space with codimension
1. Therefore, each $E_{i,\pm}^{j}$ partitions the space (excluding
$E_{i,\pm}^{j}$) into two subsets
$\\{\bm{w}_{i}^{j}\cdot{\mathbf{x}}>\pm\theta\\},\\{\bm{w}_{i}^{j}\cdot{\mathbf{x}}<\pm\theta\\}$.
Therefore, all $E_{i,\pm}^{j}$ can partition $\mathcal{D}$ into finitely many
subsets, as there are finitely many of sets $E_{i,\pm}^{j}$. ∎

######  Lemma 2.

Suppose ${\mathbb{X}}\in\mathcal{D}$. Under the assumption of Theorem 1, there
exists a numerical constant $C>0$, only depending on the data and the
constants appearing in the assumptions of the theorem, such that the following
functions and their gradients are bounded by $C$:

$\displaystyle
r^{j}_{i,l}({\mathbf{x}}),~{}j\in[m],l\in\mathcal{L}^{j},i\in\mathcal{A}(l)$
(S10) $\displaystyle\frac{\partial r^{j}_{i,j}({\mathbf{x}})}{\partial
w_{i,t}^{j}},~{}j\in[m],l\in\mathcal{L}^{j},i\in\mathcal{A}(l),t\in[p]$
$\displaystyle P^{j}(\\{{\mathbf{x}}\rightarrow
l\\}),~{}j\in[m],l\in\mathcal{L}^{j}.$

###### Proof.

First, note that $r_{i,l}^{j}({\mathbf{x}})\in[0,1]$ and therefore
$r_{i,l}^{j}({\mathbf{x}})$ is uniformly bounded. Moreover, from equation S6

$\frac{\partial
r_{i,l}^{j}({\mathbf{x}})}{\partial\bm{w}_{i}^{j}}=\begin{cases}S^{\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i\\\
-S^{\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i\end{cases}$ (S11)

with other partial derivatives of $r_{i,l}^{j}({\mathbf{x}})$ being zero. As a
result, $r_{i,l}^{j}({\mathbf{x}})$ has uniformly bounded derivative, where
the bound only depends on Assumption (A1) and the data. This is true as
$S^{\prime}$ is bounded by the assumption. This also shows that
$\frac{\partial r^{j}_{i,l}({\mathbf{x}})}{\partial w_{i,t}^{j}}$ is uniformly
bounded. Next,

$\frac{\partial^{2}r_{i,l}^{j}({\mathbf{x}})}{\partial(\bm{w}_{i}^{j})^{2}}=\begin{cases}S^{\prime\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}{\mathbf{x}}^{T}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i\\\
-S^{\prime\prime}(\bm{w}_{i}^{j}\cdot{\mathbf{x}}){\mathbf{x}}{\mathbf{x}}&\mbox{if
}~{}~{}~{}l\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i\end{cases}$ (S12)

which is similarly uniformly bounded by Assumption (A1). Therefore,
$\frac{\partial r^{j}_{i,j}({\mathbf{x}})}{\partial w_{i,t}^{j}}$ has a
uniformly bounded gradient.
Finally, $P^{j}(\\{{\mathbf{x}}\rightarrow l\\})\in[0,1]$ and by equation S5,
the gradient of $P^{j}(\\{{\mathbf{x}}\rightarrow l\\})$ is the product of
bounded functions, and therefore bounded. ∎

######  Lemma 3.

Suppose ${\mathbb{X}}\in\mathcal{D}$. Under the assumption of Theorem 1, there
exists a numerical constant $C>0$, only depending on the data and the
constants appearing in the assumptions of the theorem, such that the following
functions and their gradients are bounded by $C$:

$\displaystyle f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$ (S13)
$\displaystyle
L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))$

###### Proof.

First,

$\displaystyle|f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})|$
$\displaystyle=\left|\sum_{j=1}^{m}\sum_{l\in\mathcal{L}^{j}}o_{l}^{j}P^{j}(\\{{\mathbf{x}}\rightarrow
l\\})\right|$
$\displaystyle\leq\sum_{j=1}^{m}\sum_{l\in\mathcal{L}^{j}}|o_{l}^{j}|$
$\displaystyle\leq m|\mathcal{L}^{1}|B.$ (S14)

Moreover, from equation S3 and Lemma 2, each coordinate of the gradient of
$f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$ is finite summation and
product of uniformly bounded functions, which is bounded.
Next, as $f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$ is bounded,
$L^{\prime}({y},f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))$ is bounded
by Assumption (A2). Moreover,

$\nabla
L^{\prime}({y},f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})=L^{\prime\prime}({y},f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))\nabla
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$

which is bounded as $L^{\prime\prime}$ is bounded by Assumption (A2) and
$\nabla f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$ is bounded as we
showed above. ∎

######  Lemma 4.

There exists a numerical constant $M>0$, only depending on the data and
constants introduced in Assumptions (A1), (A2) and (A3), such that if for
${\mathbb{X}}_{1},{\mathbb{X}}_{2}$ and $\alpha\in(0,1)$,
$\alpha{\mathbb{X}}_{1}+(1-\alpha){\mathbb{X}}_{2}\in\mathcal{D}$, then one
has

$|\phi_{i,t}^{j}({\mathbb{X}}_{1})-\phi_{i,t}^{j}({\mathbb{X}}_{2})|\leq
M\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2},~{}t\in[p],j\in[m],i\in\mathcal{I}^{j}$

where $\phi_{i,t}^{j}$ is defined in equation S8.

###### Proof.

Note that by the definition of $\phi_{i,t}^{j}$,

$\displaystyle\phi_{i,t}^{j}({\mathbb{X}})$
$\displaystyle=L^{\prime}(y,f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))\frac{\partial
f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})}{\partial w_{i,t}^{j}}$

which is the product of two bounded functions with bounded derivatives by
Lemma 3. As a result, there exists a constant $M>0$ such that

$\left\|\frac{\partial\phi^{j}_{i,t}({\mathbb{X}})}{\partial{\mathbb{X}}}\right\|_{2}\leq
M.$ (S15)

For $\alpha\in[0,1]$, let
$\hat{\phi}_{i,t}^{j}(\alpha)=\phi_{i,t}^{j}((1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2})$.
Then, $\hat{\phi}_{i,t}^{j}(\alpha)$ is differentiable for $\alpha\in(0,1)$,
$\hat{\phi}_{i,t}^{j}(0)=\phi_{i,t}^{j}({\mathbb{X}}_{1})$ and
$\hat{\phi}_{i,t}^{j}(1)=\phi_{i,t}^{j}({\mathbb{X}}_{2})$. Moreover, by the
chain rule,

$\displaystyle\frac{d\hat{\phi}^{j}_{i,t}(\alpha)}{d\alpha}$
$\displaystyle=\bigg{\langle}\frac{\partial\hat{\phi}^{j}_{i,t}(\alpha)}{\partial\left[(1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2}\right]},\frac{\partial\left[(1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2}\right]}{\partial\alpha}\bigg{\rangle}$
$\displaystyle=\bigg{\langle}({\mathbb{X}}_{2}-{\mathbb{X}}_{1}),\frac{\partial\phi^{j}_{i,t}(\alpha)}{\partial\left[(1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2}\right]}\bigg{\rangle}$

where $\langle\cdot,\cdot\rangle$ denotes the inner product. As a result, by
the fundamental theorem of calculus,

$\displaystyle\left|\hat{\phi}_{i,j}^{j}(\alpha)-\hat{\phi}_{i,j}^{j}(0)\right|$
$\displaystyle=\left|\int_{0}^{\alpha}\frac{d\hat{\phi}^{j}_{i,t}(u)}{du}du\right|$
$\displaystyle\leq\|{\mathbb{X}}_{2}-{\mathbb{X}}_{1}\|_{2}\int_{0}^{\alpha}\underbrace{\left\|\frac{\partial\hat{\phi}^{j}_{i,t}(u)}{\partial\left[(1-u){\mathbb{X}}_{1}+u{\mathbb{X}}_{2}\right]}\right\|_{2}}_{\leq
M~{}\text{by~{}equation~{}\ref{app-grad-m}}}du$ $\displaystyle\leq\alpha
M\|{\mathbb{X}}_{2}-{\mathbb{X}}_{1}\|_{2}.$

In particular, by setting $\alpha=1$ the proof is complete. ∎

######  Lemma 5.

Take any two ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$ such that
$\|\bm{\mathcal{O}}_{1}\|_{2},\|\bm{\mathcal{O}}_{2}\|_{2}\leq B$. Suppose for
all $j\in[m],i\in\mathcal{I}^{j}$,
$\\{{\mathbb{X}}_{1},{\mathbb{X}}_{2}\\}\not\subseteq E_{i,+}^{j}$ or
$\\{{\mathbb{X}}_{1},{\mathbb{X}}_{2}\\}\not\subseteq E_{i,-}^{j}$. Then

$|\phi_{i,t}^{j}({\mathbb{X}}_{1})-\phi_{i,t}^{j}({\mathbb{X}}_{2})|\leq
M\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2}~{}t\in[p],j\in[m],i\in\mathcal{I}^{j}$

where $\phi_{i,t}^{j}$ is defined in equation S8 and $M$ is defined in Lemma
4.

###### Proof.

By Lemma 1, the sets $E_{i,\pm}^{j}$ are affine and therefore if both
${\mathbb{X}}_{1},{\mathbb{X}}_{2}$ do not belong to one of these sets, the
segment connecting ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$ will intersect each set
$E_{i,\pm}^{j}$ at most once. Thus, there exist $K\geq 1$ and
$\tilde{{\mathbb{X}}}_{0},\cdots,\tilde{{\mathbb{X}}}_{K}$ such that
$\tilde{{\mathbb{X}}}_{0}={\mathbb{X}}_{1},\tilde{{\mathbb{X}}}_{K}={\mathbb{X}}_{2}$,
and $\tilde{{\mathbb{X}}}_{0},\cdots,\tilde{{\mathbb{X}}}_{K}$ lie in the
segment connecting ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$. Moreover, we have that

$\\{(1-\alpha)\tilde{{\mathbb{X}}}_{k}+\alpha\tilde{{\mathbb{X}}}_{k+1}:\alpha\in(0,1)\\}\subseteq\mathcal{D}~{}~{}\forall
k\in[K-1].$

By triangle inequality, we have

$\displaystyle|\phi_{i,t}^{j}({\mathbb{X}}_{1})-\phi_{i,t}^{j}({\mathbb{X}}_{2})|$
$\displaystyle=\left|\sum_{k=0}^{K-1}\\{\phi_{i,t}^{j}(\tilde{{\mathbb{X}}}_{k})-\phi_{i,t}^{j}(\tilde{{\mathbb{X}}}_{k+1})\\}\right|$
$\displaystyle\leq\sum_{k=0}^{K-1}\left|\phi_{i,t}^{j}(\tilde{{\mathbb{X}}}_{k})-\phi_{i,t}^{j}(\tilde{{\mathbb{X}}}_{k+1})\right|$
$\displaystyle\leq
M\sum_{k=0}^{K-1}\|\tilde{{\mathbb{X}}}_{k}-\tilde{{\mathbb{X}}}_{k+1}\|_{2}$
$\displaystyle=M\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2}$

where the last inequality is by Lemma 4 and the last equality is by the fact
that $\tilde{{\mathbb{X}}}_{0},\cdots,\tilde{{\mathbb{X}}}_{K}$ lie in the
segment connecting ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$. ∎

######  Lemma 6.

There exists a numerical constant $\mathcal{M}>0$, only depending on the data
and constants introduced in Assumptions (A1), (A2) and (A3), such that
$\Phi(\bm{\mathcal{O}},\bm{\mathcal{W}})$ is $\mathcal{M}$-smooth for
$\|\bm{O}\|_{2}\leq B$.

###### Proof.

First, note that

$\|\nabla\Phi(\bm{\mathcal{O}}_{1},\bm{\mathcal{W}}_{1})-\nabla\Phi(\bm{\mathcal{O}}_{2},\bm{\mathcal{W}}_{2})\|_{2}\leq\sum_{j=1}^{m}\sum_{l\in\mathcal{L}^{j}}|\phi^{j}_{l}({\mathbb{X}}_{1})-\phi^{j}_{l}({\mathbb{X}}_{2})|+\sum_{j=1}^{m}\sum_{i\in\mathcal{I}^{j}}\sum_{t=1}^{p}|\phi^{j}_{i,t}({\mathbb{X}}_{1})-\phi^{j}_{i,t}({\mathbb{X}}_{2})|.$
(S16)

Take ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$ with
$\|\bm{O}_{1}\|_{2},\|\bm{O}_{2}\|_{2}\leq B$. If these two solutions
simultaneously do not belong to some $E_{i,\pm}^{j}$, by Lemma 5 we have

$|\phi^{j}_{i,t}({\mathbb{X}}_{1})-\phi^{j}_{i,t}({\mathbb{X}}_{2})|\leq
M\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2}.$

A similar result follows for $\phi^{j}_{l}$, hence by equation S16 we achieve

$\|\nabla\Phi(\bm{\mathcal{O}}_{1},\bm{\mathcal{W}}_{1})-\nabla\Phi(\bm{\mathcal{O}}_{2},\bm{\mathcal{W}}_{2})\|_{2}\leq
m(|\mathcal{L}^{1}|+p|\mathcal{I}^{1}||)M\|{\mathbb{X}}_{1}-{\mathbb{X}}_{2}\|_{2}$

which completes the proof. Suppose there exists $i_{0},j_{0}$ such that
${\mathbb{X}}_{1},{\mathbb{X}}_{2}\in E^{j_{0}}_{i_{0},+}$. Then,

$\left((1-\alpha)(\bm{w}_{1})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}+\left(\alpha(\bm{w}_{2})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}=\theta$

therefore
$S\left(\left((1-\alpha)(\bm{w}_{1})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}+\left(\alpha(\bm{w}_{2})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}\right)=1$
by Assumption (A1). As a result,
$S(\bm{w}_{i_{0}}^{j_{0}}\cdot{\mathbf{x}})=1$ for the whole segment
connecting ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$. Let
$\tilde{{\mathbb{X}}}_{1},\tilde{{\mathbb{X}}}_{2}$ be such that all their
weights are the same as ${\mathbb{X}}_{1},{\mathbb{X}}_{2}$, except that

$(\tilde{\bm{w}}_{1})_{i_{0}}^{j_{0}}=2({\bm{w}}_{1})_{i_{0}}^{j_{0}},(\tilde{\bm{w}}_{2})_{i_{0}}^{j_{0}}=2({\bm{w}}_{2})_{i_{0}}^{j_{0}}.$

As a result,

$\left((1-\alpha)(\tilde{\bm{w}}_{1})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}+\left(\alpha(\tilde{\bm{w}}_{2})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}=2\theta$

so
$S\left(\left((1-\alpha)(\tilde{\bm{w}}_{1})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}+\left(\alpha(\tilde{\bm{w}}_{2})_{i_{0}}^{j_{0}}\right)\cdot{\mathbf{x}}\right)=1$.
Therefore, one has

$\displaystyle{f}({\mathbf{x}};(1-\alpha)\tilde{{\mathbb{X}}}_{1}+\alpha\tilde{{\mathbb{X}}}_{2})$
$\displaystyle={f}({\mathbf{x}};(1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2}),~{}\alpha\in[0,1],$
$\displaystyle\frac{\partial
f({\mathbf{x}};(1-\alpha)\tilde{{\mathbb{X}}}_{1}+\alpha\tilde{{\mathbb{X}}}_{2})}{\partial\bm{w}_{i}^{j}}$
$\displaystyle=\frac{\partial
f({\mathbf{x}};(1-\alpha){\mathbb{X}}_{1}+\alpha{\mathbb{X}}_{2})}{\partial\bm{w}_{i}^{j}},~{}\alpha\in[0,1],j\in[m],i\in\mathcal{I}^{j}$

as for all $i,j$,

$S\left(\left((1-\alpha)(\tilde{\bm{w}}_{1})_{i}^{j}\right)\cdot{\mathbf{x}}+\left(\alpha(\tilde{\bm{w}}_{2})_{i}^{j}\right)\cdot{\mathbf{x}}\right)=S\left(\left((1-\alpha)({\bm{w}}_{1})_{i}^{j}\right)\cdot{\mathbf{x}}+\left(\alpha({\bm{w}}_{2})_{i}^{j}\right)\cdot{\mathbf{x}}\right).$

As a result, for all $i,j,t$,

$\phi^{j}_{i,t}(\tilde{{\mathbb{X}}}_{1})=\phi^{j}_{i,t}({{\mathbb{X}}}_{1}),\phi^{j}_{i,t}(\tilde{{\mathbb{X}}}_{2})=\phi^{j}_{i,t}({{\mathbb{X}}}_{2})$

However, $\tilde{{\mathbb{X}}}_{1},\tilde{{\mathbb{X}}}_{2}\notin
E_{i_{0},+}^{j_{0}}$. In words, the new points effectively replace the
problematic coefficient $\bm{w}_{i_{0}}^{j_{0}}$ in the model coefficients.
Repeat this process until all such coefficients are removed and therefore
$\tilde{{\mathbb{X}}}_{1},\tilde{{\mathbb{X}}}_{2}$ fit into the assumptions
of Lemma 5. This completes the proof. ∎

######  Lemma 7.

Under Assumption (A1), the activation function $S(x)$ is semi-algebraic.

###### Proof.

Consider the epigraph of the activation function

$\mathcal{G}(S)=\\{(x,y):y\geq S(x)\\}.$

Let

$\displaystyle A_{1}$ $\displaystyle=\\{y\geq 0,x\leq-\theta\\}$
$\displaystyle A_{2}$ $\displaystyle=\\{y\geq 1,x\geq\theta\\}$ $\displaystyle
A_{3}$ $\displaystyle=\\{y\geq p(x),-\theta\leq x\leq\theta\\}.$

Then,

$\mathcal{G}(S)=A_{1}\cup A_{2}\cup A_{3}$

showing $\mathcal{G}(S)$ and consequently, $S$ are semi-algebraic. ∎

######  Lemma 8.

The function $1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]$ is semi-algebraic.

###### Proof.

Consider the epigraph:

$G=\mathcal{G}(1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}])=\\{(y,\bm{\mathcal{W}}):y\geq
1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]\\}.$

Let $A_{1}=\\{y\geq 1\\}$ and
$A_{2}=\\{\|\bm{\mathcal{W}}_{k,:,:}\|_{2}^{2}=0,y\geq 0\\}.$ Then,
$G=A_{1}\cup A_{2}$ showing $1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]$ is semi-
algebraic. ∎

###### Proof of Theorem 1..

Part 1)  Note that by Assumption (A3), Algorithm 1 can be equivalently run on
the problem

${\min}_{\bm{\mathcal{W}},\bm{\mathcal{O}}}~{}~{}\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))]+\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]+(\lambda_{2}/m|\mathcal{I}|)\|\bm{\mathcal{W}}\|_{2}^{2}~{}~{}~{}~{}~{}~{}~{}\text{s.t.}~{}~{}~{}~{}~{}~{}~{}\|\bm{\mathcal{O}}\|_{2}\leq
B$ (S17)

and result in the same sequence of solutions. Let

$\displaystyle h(\bm{\mathcal{W}},\bm{\mathcal{O}})$
$\displaystyle=\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))]+(\lambda_{2}/m|\mathcal{I}|)\|\bm{\mathcal{W}}\|_{2}^{2}$
(S18) $\displaystyle g(\bm{\mathcal{W}},\bm{\mathcal{O}})$
$\displaystyle=\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]+\chi(\|\bm{\mathcal{O}}\|_{2}\leq
B)$

where $\chi(\|\bm{\mathcal{O}}\|_{2}\leq B)=0$ if
$\|\bm{\mathcal{O}}\|_{2}\leq B$ and $\chi(\|\bm{\mathcal{O}}\|_{2}\leq
B)=\infty$ otherwise. Then, Problem S17 can be written as minimizing $g+h$. By
Lemma 6, the function $h$ is
$(\mathcal{M}+2\lambda_{2}/m|\mathcal{I}|)$-smooth. Hence, if
$\eta<1/(\mathcal{M}+2\lambda_{2}/m|\mathcal{I}|)$, by Lemma 3.1 of [4] (and
calculations leading to (52) of [4]) the descent property is proved.

Part 2)  Let $c_{0}$ be the objective value for the initial solution to
Algorithm 1. That is,

$\frac{1}{N}\sum_{n=1}^{N}L(y_{n},\bm{f}({\mathbf{x}}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}}))]+\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]\\\
+(\lambda_{2}/m|\mathcal{I}|)\|\bm{\mathcal{W}}\|_{2}^{2}=c_{0}.$

As the sequence of the objectives is non-increasing by the first part of the
theorem, at each iteration we have

$\displaystyle(\lambda_{2}/m|\mathcal{I}|)\|\bm{\mathcal{W}}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}L(y_{n},\bm{f}({\mathbf{x}}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}}))]+\lambda_{0}{\sum}_{k\in[p]}1[\bm{\mathcal{W}}_{k,:,:}\neq\bm{0}]+(\lambda_{2}/m|\mathcal{I}|)\|\bm{\mathcal{W}}\|_{2}^{2}$
$\displaystyle\leq$ $\displaystyle c_{0}$

completing the proof as $\lambda_{2}>0$.

Part 3)  Note that the function
$f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})$ is a composition, finite
sum and finite product of semi-algebraic functions by Lemma 7, therefore it is
semi-algebraic. As $L$ is polynomial,
$L(f({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}}))$ is semi-algebraic. As a
result, by Lemma 8 the functions $g,h$ defined in equation S18 are semi-
algebraic and hence, possess the KL property. Therefore, $g,h$ satisfy the
required conditions of Theorem 5.1 of [4], proving the convergence. ∎

##  Appendix S4 Appendix for Sec. 6

Data Design We generate the data matrix, $\bm{X}\in\mathbb{R}^{N\times p}$
with samples randomly drawn from a multivariate normal distribution
$\mathcal{N}(0,\bm{\Sigma})$ with a correlated design matrix $\bm{\Sigma}$
whose values are defined by $\Sigma_{ij}=\sigma^{|i-j|}$. We construct the
response variable $\bm{y}=\bm{X}\bm{\beta}^{*}+\epsilon$ where the values of
the noise $\epsilon_{i},i=[N]$ are drawn independently from
$\mathcal{N}(0,0.5)$. We use a known sparsity $\bm{\beta}^{*}$ with equi-
spaced nonzeros, where $\left\lVert\bm{\beta}^{*}\right\rVert_{0}=8$.

Simulation Procedure We experiment with a range of training set sizes
$N\in\\{100,200,1000\\}$, correlation strengths $\sigma\in\\{0.5,0.7\\}$, and
number of total features $p\in\\{256,512\\}$. For each setting, we run 25
simulations with randomly generated samples. We evaluate on $10,000$ test
samples drawn from the data generating model described above. Out of $N$
samples, we allocate $80\%$ for training and $20\%$ for validation for model
selection. For each simulation, we perform 500 tuning trials (hyperparameters
are given below) and select the model with the smallest validation mean
squared error (mse). We evaluate the final performance in terms of (i) test
MSE, (ii) number of features selected and (iii) support recovery via f1-score
between the true support and the recovered feature set. We compute averages
and standard errors across the 25 simulations to report final results.

Skinny Trees with DSL

  * •

Number of trees: Discrete uniform with range $[1,50]$,

  * •

Depths: Discrete uniform with range $[1,5]$,

  * •

Batch sizes: $16b$ with b uniform over the range $\\{1,2,\ldots,8\\}$,

  * •

Number of Epochs: Discrete uniform with range $[5,500]$,

  * •

$\lambda_{0}$: exponentially ramped up $0\rightarrow 100$ with temperature
distributed as Log uniform in the range $[10^{-4},0.1]$ for group
$\ell_{0}-\ell_{2}$ with DSL,

  * •

$\lambda_{2}$: Log uniform over the range $[10^{-2},10^{2}]$ for group
$\ell_{0}-\ell_{2}$,

  * •

Learning rates, $lr$: Log uniform over the range $[10^{-3},10^{-1}]$.

XGBoost, Random Forests

  * •

Number of trees: Discrete uniform with range $[1,100]$ for XGBoost and Random
Forests,

  * •

Depths: Discrete uniform with range $[1,10]$ for XGBoost and Random Forests,

  * •

Learning rates: Discrete uniform with range $[10^{-4},1.0]$ for XGBoost,

  * •

Feature importance threshold: Log uniform over the range $[10^{-7},10^{-1}]$
for XGBoost and Random Forests.

  * •

Subsample: Uniform over the range $[0.5,1.0]$.

##  Appendix S5 Appendix for Sec. 7

###  S5.1 Datasets, Computing Setup and Tuning Setup

Table S3: Classification datasets Dataset | N | p | C
---|---|---|---
Churn | 5,000 | 20 | 2
Satimage | 6,435 | 36 | 6
Texture | 5,500 | 40 | 11
Mice-protein | 1,080 | 77 | 8
Isolet | 7,797 | 617 | 26
Madelon | 2,600 | 500 | 2
Lung | 203 | 3,312 | 5
Gisette | 7,000 | 5,000 | 2
Tox | 171 | 5,748 | 4
Arcene | 200 | 10,000 | 2
CLL | 111 | 11,340 | 3
SMK | 187 | 19,993 | 2
GLI | 85 | 22,283 | 2
Dorothea | 1,150 | 100,000 | 2

#### Datasets

We consider 5 classification datasets from the Penn Machine Learning
Benchmarks (PMLB) repository [51]. These are Churn, Satimage, Mice-protein,
Isolet and Texture. ; Mice-protein and Isolet were used in prior feature
selection literature [43]. We used fetch_openml in Sklearn to download the
full datasets. We used 4 datasets from NIPS2003 feature selection challenge
[25]. These are Arcene, Madelon, Gisette, Dorothea. We used 5 datasets (Smk,
Cll, Gli, Lung, Tox) from the feature selection datasets given in this
repo111https://jundongl.github.io/scikit-feature/datasets.html. The datasets
contain continuous, categorical and binary features. The metadata was used to
identify the type of features and categorical features were one-hot encoded.
For first non-NIPS2003 datasets, we randomly split each of the dataset into
64% training, 16% validation and 20% testing sets. For the datasets from
NIPS2003, we split the training set into 80% training and 20% validation and
treated the original validation set as the test set. The labels for original
test set were unavailable, hence we discarded the original test set. The
continuous features in all datasets were z-normalized based on training set
statistics. A summary of the 14 datasets considered is in Table S3.

Computing Setup We used a cluster running Ubuntu 7.5.0 and equipped with Intel
Xeon Platinum 8260 CPUs and Nvidia Volta V100 GPUs. For all experiments of
Sec. 6 and 7, each job involving Skinny Trees, Random Forests, XGBoost,
LightGBM, CatBoost and neural networks were restricted to 4 cores and 64GB of
RAM. Jobs involving Skinny Trees and neural networks on larger datasets
(Gisette, Dorothea) were run on Tesla V100 GPUs.

Tuning The tuning was done in parallel over the competing models and datasets.
The number of selected features affects the AUC. Therefore, to treat all the
methods in a fair manner, we tune the hyperparameter that controls the
sparsity level using Optuna [2] which optimizes the overall AUC across
different Ks. A list of all the tuning parameters and their distributions is
given for every experiment below.

###  S5.2 Appendix for Sec. 7.1

Classical tree: TAO Implementation Given that the authors of TAO do not open-
source their implementation, we have written our own implementation of the TAO
algorithm proposed by [13], following the procedure outlined in Sec. 3.1 of
[73]. For TAO, we use oblique (i.e. hyperplane splits) decision trees with
constant leaves. We take as an initial tree a complete binary tree with random
parameters at each node. We perform TAO updates until maximum number of
iterations are reached (i.e. there is no other stopping criterion). TAO uses
an $\ell_{1}$-regularized logistic regression to solve the decision node
optimization (using LIBLINEAR [22]) where $\lambda\geq 0$ parameter
(controlling node-level sparsity of the tree) is used as a regularization
parameter ($C=1/\lambda$). We tune depth in the range $[1,10]$,
$\lambda\in[10^{-5},10^{5}]$ and number of maximum iterations in the range
$[20,100]$. We perform 500 tuning trials. We find the optimal trial that
satisfies $50\%$ sparsity budget.

Soft tree with group lasso We compare group $\ell_{0}$-based method with a
competitive benchmark: the convex group lasso regulurizer, popularly used in
high-dimensional statistics literature [70]. We consider group-lasso
regularization in the context of soft trees:
$(\lambda_{1}/\sqrt{m|\mathcal{I}|}){\sum}_{k\in[p]}\left\lVert\bm{\mathcal{W}}_{k,:,:}\right\rVert_{2}.$
Although group lasso has not been used in soft trees, it has been used for
feature selection in neural networks [61]. However, their proposal to use GD
on a group $\ell_{1}$ regularized objective does not lead to feature selection
as highlighted by [43]. For a fairer comparison, we implement our own proximal
mini-batch GD method for group lasso in the context of soft trees, which
actually leads to feature selection.

We consider the following optimization problem with group-lasso regularization
in the context of soft trees:

$\displaystyle{\min}_{\bm{\mathcal{W}},\bm{\mathcal{O}}}~{}~{}\hat{\mathbb{E}}[L({\mathbf{y}},\bm{f}({\mathbf{x}};\bm{\mathcal{W}},\bm{\mathcal{O}})]+(\lambda_{1}/\sqrt{m|\mathcal{I}|}){\sum}_{k\in[p]}\left\lVert\bm{\mathcal{W}}_{k,:,:}\right\rVert_{2}.$
(S19)

where $\lambda_{1}$ is a non-negative regularization parameter that controls
both shrinkage and sparsity. We solve it with the algorithm presented in Alg.
2.

1:Data: $X,Y$; Hyperparameters: $\lambda_{1}$, epochs, batch-size ($b$),
learning rate ($\eta$);

2:Initialize ensemble with $m$ trees of depth $d$ ($|\mathcal{I}|=2^{d}-1$):
$\bm{\mathcal{W}},\bm{\mathcal{O}}$

3:for $\text{epoch}=1,2,\ldots,\text{epochs}$ do

4: for $batch=1,\ldots,N/b$ do

5: Randomly sample a batch: $\bm{X}_{batch},\bm{Y}_{batch}$

6: Compute gradient of loss $g$ w.r.t. $\bm{\mathcal{O}}$, W, where
$g=\hat{\mathbb{E}}[L(\bm{Y}_{batch},\bm{F}_{batch}]$.

7: Update leaves:
$\bm{\mathcal{O}}\leftarrow\bm{\mathcal{O}}-\eta\nabla_{\bm{\mathcal{O}}}g$

8: Update hyperplanes:
$\bm{\mathcal{W}}\leftarrow\textit{S-Prox}(\bm{\mathcal{W}}-\eta\nabla_{\bm{\mathcal{W}}}g,\eta,\lambda_{1}))$

9: end for

10:end for

Algorithm 2 Proximal Mini-batch Gradient Descent for Optimizing (S19).

S-Prox in Algorithm 2 finds the global minimum of an optimization problem of
the form:

$\displaystyle\bm{\mathcal{W}}^{(t)}=\text{argmin}_{\bm{\mathcal{W}}}\frac{1}{2\eta}\left\lVert\bm{\mathcal{W}}-\bm{\mathcal{Z}}^{(t)}\right\rVert_{2}^{2}+\frac{\lambda_{1}}{\sqrt{m|\mathcal{I}|}}\sum_{k=1}^{p}\left\lVert\bm{\mathcal{W}}_{k,:,:}\right\rVert_{2}$
(S20)

where
$\bm{\mathcal{Z}}^{(t)}=\bm{\mathcal{W}}^{(t-1)}-\eta\nabla_{\bm{\mathcal{W}}}g$
and $g=\hat{\mathbb{E}}[L(\bm{Y}_{batch},\bm{F}_{batch}]$. This leads to a
feature-wise soft-thresholding operator given by:

$\displaystyle
S_{\frac{\eta\lambda_{1}}{\sqrt{m|\mathcal{I}|}}}(\bm{\mathcal{Z}}_{k,:,:}^{(t)})=\begin{cases}\bm{\mathcal{Z}}_{k,:,:}^{(t)}-\frac{\eta\lambda_{1}}{\sqrt{m|\mathcal{I}|}}\frac{\bm{\mathcal{Z}}_{k,:,:}^{(t)}}{\left\lVert\bm{\mathcal{Z}}_{k,:,:}^{(t)}\right\rVert}&\text{if
$\left\lVert\bm{\mathcal{Z}}_{k,:,:}^{(t)}\right\rVert\geq\frac{\eta\lambda_{1}}{\sqrt{m|\mathcal{I}|}}$}\\\
0&\text{otherwise}\end{cases}$ (S21)

In these experiments, we used a single tree and mini-batch PGD without any
dense-to-sparse scheduler for both models i.e, (i) Soft tree with group lasso
(ii) Skinny Tree. We tune the key hyperparameters, which are given below.

  * •

Batch sizes: $64*b$ with $b$ uniform over the range $\\{1,2,\ldots,64\\}$,

  * •

Learning rates for mini-batch PGD: Log uniform over the range $[10^{-2},10]$,

  * •

Number of Epochs: Discrete uniform with range $[5,1000]$,

  * •

$\lambda_{0}$: Log uniform over the range $[10^{-3},10]$ for group
$\ell_{0}-\ell_{2}$ for Skinny Tree,

  * •

$\lambda_{2}$: Log uniform over the range $[10^{-2},10^{2}]$ for group
$\ell_{0}-\ell_{2}$ for Skinny Tree,

  * •

$\lambda_{1}$: Log uniform over the range $[10^{-3},10]$ for group lasso.

###  S5.3 Appendix for Sec. 7.2 and 7.3

We describe the tuning grid for these experiments below.

Skinny Trees

  * •

Number of trees: Discrete uniform with range $[1,100]$,

  * •

Depths: Discrete uniform with range $[1,5]$,

  * •

Batch sizes: Uniform over the set $\\{16,64,256,1024\\}$,

  * •

Number of Epochs: Discrete uniform over the range $[5,1000]$ (Note tuning over
epochs is important to achieve some trials with desired sparsity for dense-to-
sparse learning, we do not use any validation loss early stopping as that is
less robust during averages across runs/seeds in terms of feature support.),

  * •

$\lambda_{0}$: exponentially ramped up $0\rightarrow 1$ with temperature
distributed as Log uniform in the range $[10^{-4},1]$ for group
$\ell_{0}-\ell_{2}$ with DSL,

  * •

$\lambda_{2}$: Log uniform over the range $[10^{-3},1]$ for group
$\ell_{0}-\ell_{2}$,

  * •

Learning rate ($lr$) for minibatch PGD: Log uniform over the range
$[10^{-2},10]$.

XGBoost, LightGBM, CatBoost, Random Forests

  * •

Number of trees: Discrete uniform with range $[1,300]$ for XGBoost, LightGBM,
CatBoost and Random Forests,

  * •

Depths: Discrete uniform with range $[1,10]$ for XGBoost, LightGBM, CatBoost
and Random Forests,

  * •

Learning rates: Discrete uniform with range $[10^{-4},1.0]$ for XGBoost,
LightGBM and CatBoost,

  * •

Feature importance threshold: Log uniform over the range $[10^{-7},10^{-1}]$
for XGBoost, LightGBM, CatBoost and Random Forests,

  * •

Subsample: Uniform over the range $[0.5,0.9]$ for XGBoost, LightGBM, CatBoost
and Random Forests.

  * •

Bagging Frequency: Uniform over the set $\\{1,\cdots,7\\}$ for LightGBM.

FASTEL (Dense soft trees)

  * •

Number of trees: Discrete uniform with range $[1,100]$,

  * •

Depths: Discrete uniform with range $[1,5]$,

  * •

Batch sizes: Uniform over the set $\\{16,64,256,1024\\}$,

  * •

Number of Epochs: Discrete uniform over the range $[5,1000]$,

  * •

Learning rates ($lr$) for minibatch SGD: Log uniform over the range
$[10^{-2},10]$.

###  S5.4 Appendix for Sec. 7.4

In this paper, we pursue embedded feature selection methods for tree
ensembles. However, for completeness, we also compare Skinny Trees against
some state-of-the-art embedded feature selection methods from neural networks.

Competing Methods We compare against the following baselines.

  1. 1.

LassoNet [43] is based on a ResNet-like [30] architecture with residual
connections. Feature selection is done using hierarchical group lasso.

  2. 2.

AlgNet [18] is based on adaptive lasso and uses proximal full-batch gradient
descent for feature selection in a tanh-activated feedforward network.

  3. 3.

DFS [17] solves cardinality constrained feature selection problem with an
active-set style method.

We describe the tuning details for these neural network models below. We
perform 2000 tuning trials for each method. For DFS, we capped number of
trials completed in total 200 GPU hours. DFS is really slow for even medium
sized datasets (in terms of feature dimensions).

#### LassoNet [43]

  * •

ResNet architecture with 2-layered relu-activated feedforward network with
(linear) skip connections.

  * •

Number of hidden units: Discrete uniform in the set
$\\{\frac{p}{3},\frac{2}{3}p,p,\frac{4}{3}p\\}$,

  * •

Batch sizes: Discrete uniform in the set $\\{64,128,256,512\\}$,

  * •

$\lambda$: Log uniform over the range $[1,10000]$.

  * •

Tuning protocol: $1000$ for dense training stage and $1000$ for each
successive sparse training stages with early stopping with a patience of 25.
We consider 100 sequential stages (100 values for $\lambda$).

AlgNet [18]

  * •

Tanh-activated 4-layered feedforward neural network with number of hidden
units chosen uniformly from a discrete set
$\\{\frac{p}{3},\frac{2}{3}p,p,\frac{4}{3}p\\}$,

  * •

Learning rates ($lr$) for proximal GD: Log uniform over the range $[0.01,1]$.

  * •

$\lambda$: Log uniform over the range $[0.1,1000]$.

  * •

Number of Epochs: Discrete uniform with range $[5,2000]$,

DFS [17]

  * •

ReLU-activated 4-layered feedforward neural network with number of hidden
units chosen uniformly from a discrete set
$\\{\frac{p}{3},\frac{2}{3}p,p,\frac{4}{3}p\\}$,

  * •

Learning rates ($lr$) for Adam: Log uniform over the range $[0.001,0.1]$.

  * •

Weight decay: $0.0025$ for feature selection layer and $0.005$ for remaining
layers.

  * •

$k$: Number of features selected in the range $[1,0.5p]$.

  * •

Number of Epochs: $500$ with early stopping.

###  S5.5 Appendix for Sec. 7.5

We perform 2000 tuning trials for each method.

Skinny Trees

  * •

Number of trees: Discrete uniform with range $[1,100]$,

  * •

Depths: Discrete uniform with range $[1,5]$,

  * •

Batch sizes: Uniform over the set $\\{16,64,256,1024\\}$,

  * •

Number of Epochs: Discrete uniform with range $[5,1000]$,

  * •

$\lambda_{0}$ (without Dense-to-sparse learning): Log uniform in the range
$[1,10^{4}]$ for group $\ell_{0}-\ell_{2}$,

  * •

$\lambda_{0}$ (with Dense-to-sparse learning): exponentially ramped up
$0\rightarrow 1$ with temperature $s$ distributed as Log uniform in the range
$[10^{-4},1]$ for group $\ell_{0}-\ell_{2}$ with DSL,

  * •

$\lambda_{2}$: Log uniform over the range $[10^{-3},1]$ for group
$\ell_{0}-\ell_{2}$,

  * •

Learning rates ($lr$): Log uniform over the range $[10^{-2},10]$.

##  Appendix S6 Measuring statistical significance

We follow the following procedure to test the significance of all models. For
all models, we tune over hyperparameters for 2000 trials. We select the
optimal trial (within the desired feature sparsity budget) based on validation
set. Next, we train each model for 100 repetitions with the optimal
hyperparameters and report the mean results on test set along with the
standard errors.

Table S4: Comparison of test AUC (%) performance of Skinny Trees against
ControlBurn for $25\%$ feature selection budget on binary classification
tasks.

Dataset |  ControlBurn |  Skinny Trees
---|---|---
Churn |  85.66$\pm$0.25 |  91.38$\pm$0.08
Madelon |  94.23$\pm$0.08 |  94.14$\pm$0.09
Gisette |  99.36$\pm$0.01 |  99.81$\pm$0.002
Arcene |  77.89$\pm$0.32 |  80.80$\pm$0.30
Smk |  79.94$\pm$0.32 |  79.29$\pm$0.22
Gli |  98.62$\pm$0.16 |  99.80$\pm$0.07
Dorothea |  84.85$\pm$0.28 |  90.87$\pm$0.02
Average |  88.65 |  90.87

##  Appendix S7 Comparison with ControlBurn

We also compare against ControlBurn [45] on binary classification tasks.
ControlBurn does not support multiclass classification. ControlBurn is another
post-training feature selection method, which formulates an optimization
problem with group-lasso regularizer to perform feature selection on a pre-
trained forest. It relies on a commercial solver to optimize their
formulation. We report the results for $25\%$ feature selection budget in
Table S4. We can observe that Skinny Trees outperforms ControlBurn across many
datasets, achieving $2\%$ (up to $6\%$) improvement in AUC.

We note the hyperparameter tuning for ControlBurn below:

  * •

Method: bagboost,

  * •

Depths: Discrete uniform with range $[1,10]$,

  * •

$\alpha$: Log uniform over the range $[10^{-7},1.0]$.