# Linear Iterative Feature Embedding: An Ensemble Framework for Interpretable
Model

Agus Sudjianto, Jinwen Qiu 111Corresponding author<EMAIL_ADDRESS>Miaoqi Li, Jie Chen

(Wells Fargo)

###### Abstract

A new ensemble framework for interpretable model called Linear Iterative
Feature Embedding (LIFE) has been developed to achieve high prediction
accuracy, easy interpretation and efficient computation simultaneously. The
LIFE algorithm is able to fit a wide single-hidden-layer neural network (NN)
accurately with three steps: defining the subsets of a dataset by the linear
projections of neural nodes, creating the features from multiple narrow
single-hidden-layer NNs trained on the different subsets of the data,
combining the features with a linear model. The theoretical rationale behind
LIFE is also provided by the connection to the loss ambiguity decomposition of
stack ensemble methods. Both simulation and empirical experiments confirm that
LIFE consistently outperforms directly trained single-hidden-layer NNs and
also outperforms many other benchmark models, including multi-layers Feed
Forward Neural Network (FFNN), Xgboost, and Random Forest (RF) in many
experiments. As a wide single-hidden-layer NN, LIFE is intrinsically
interpretable. Meanwhile, both variable importance and global main and
interaction effects can be easily created and visualized. In addition, the
parallel nature of the base learner building makes LIFE computationally
efficient by leveraging parallel computing.

Keywords: Linear Iterative Feature Embedding, Ensemble Method, Loss
Decomposition, Variable Importance, Interaction Detection.

##  1 Introduction

Ensemble methods have proved successful in the majority of machine learning
competitions, as they integrate multiple machine learning algorithms into one
predictive model. In particular, there are three main types of ensemble
methods, including bootstrap aggregating (bagging), boosting, and stacking.
Bagging and stacking learns base learners independently and aggregates them
following a deterministic averaging process, while boosting learns
sequentially in a very adaptive way and combines learners using a pre-
specified strategy. The main goal of bagging is to arrive at an ensemble
method with less variance than its base learners, whereas boosting mainly try
to produce a strong model that is less biased than their base learners. Stack
ensemble aims at reduce both variance and bias by combining diversified base
learners. Compared with a single model, all of these ensemble methods can
significantly improve predictive performance either by bias or variance
reduction.

However, the final model produced by ensemble method is still regarded as a
black box model, since the combination of base learners leads to a complicated
model structure and makes inner decision-making process not transparent for
human beings. The ability to explain the rationale behind one’s decisions to
others is an important aspect of human intelligence in either social
interaction or educational context. The interpretability of the results to
enable business owners or regulators to better understand risk management
decision processes and compel companies to meet regulatory requirements. Some
commonly used interpretable models, such as linear model and general additive
model cannot compete with ensemble models including Xgboost [10] and Random
Forest [4] since its simple structure cannot capture complicated dynamic data
patterns.

In recent decades, some research works focus on enhancing interpretability by
developing tools to “open up the black box”. There are, broadly speaking,
three inter-related model–based areas of research: a) global diagnostics
(Sobol & Kucherenko (2009) [16], Kucherenko (2010) [17]); b) local diagnostics
(Sundararajan et al. (2017) [24], Ancona et al. (2018) [2]); and c)
development of approximate or surrogate models that may be easier to
understand and explain. However, (Rudin and Cynthia (2019) [22] ) suggests to
avoid using explainable black-box models in high-stakes decisions since they
are sometimes problematic with several reasons such as unreliable presentation
or lack of details of what the original model delivers. Therefore, other
researchers made efforts to build inherently interpretable model, such as
explainable neural network (xNN) proposed by (Vaughan et al. (2018) [26]),
adaptive explainable neural networks (AxNNs) by (Chen et al. (2020) [8]), and
explainable neural networks with constraint (Yang et al. (2020) [28]).
Following this research direction, we try to build an inherently interpretable
with a predictive performance as strong as some black box models or an even
better performance.

In this paper, our LIFE algorithm fulfills three main goals: competitive
predictive performance, boosted computation efficiency, and interpretable
model. We know that the single-hidden-layer NN has universal approximation
property in theory and is easy to be interpreted due to simple architecture.
However, we need to resort to a wide single-hidden-layer NN with a large
number of neural nodes to obtain a strong predictive performance, which is
numerically difficult to estimate in practice. Therefore, by leveraging the
ensemble method and a simple structure of a single-hidden-layer NN, we
developed an innovative and flexible framework called LIFE to train wide
single-hidden-layer NNs, which can achieve both high accuracy and easy
interpretability in both regression and classification settings.

In this algorithm, we first use a special hierarchical structure of multiple
single layer NNs to perform data sampling based on the linear projection of
neurons, and then train multiple narrow single-hidden-layer NNs with ReLU
activation as base learners on different subsets of the dataset; finally, we
aggregate neural nodes as features from multiple base learners into a wide
single layer NN and do a join estimation via a linear model. Compared with the
traditional training strategy, LIFE effectively avoids directly training a
wide single layer NN by extracting features from multiple narrow single-
hidden-layer NNs trained on different subsets of dataset instead. In this way,
we can introduce diversity among the base learners, which tends to decrease
the total uncertainty after ensemble and, thus, yields better results
empirically. To achieve good diversity among base learners, we build a
hierarchical structure of multiple single neural networks, and leverage the
linear projections of the neurons to define the subset of sampling. In
addition, the algorithm combines the features defined by neurons from the
single layer NN base learner, for which we called neurons flattening, and this
can further improve results compared with traditional ensemble methods. This
technique also allows LIFE to take advantage of parallel computing to improve
computational efficiency through training multiple narrow single-hidden-layer
NNs simultaneously.

Extensive analyses on both simulated data and public real data verify its
effectiveness in both predictive and computational performance. Furthermore,
we provide theoretical foundation for LIFE and prove the importance of the
diversity of base learners by exploring the relationship with two-stage
ensemble stacking and the ambiguity loss decomposition for two-stage ensemble
stacking. The final single-hidden-layer NN architecture obtained from LIFE
allows to visualize the neural network weights and bias and understand the
input and output relationship easily. In particular, a single-hidden-layer NN
with rectified linear unit (ReLU) activation function [26] is equivalent to an
additive index model with linear splines on linear projections. Moreover, it
can be considered a local linear model, where all predictors are easily
visualized by a parallel coordinates plot [3]. The main and interaction
effects can also be identified by aggregating local linear model coefficients.

In general, our main contributions are summarized below:

  1. 1.

LIFE is an innovative and flexible framework for ensemble methods, which
allows different kinds of variants.

  2. 2.

LIFE can achieve a better predictive performance than traditional single-
hidden-layer NN training methods, as demonstrated by the theoretical
background and empirical experiments.

  3. 3.

LIFE can still keep model interpretable, and a new interpretation tool is
introduced to detect main and interaction effects.

  4. 4.

LIFE can improve the computation efficiency via easy parallelization,
rendering wide single layer NN training faster.

  5. 5.

An theoretical foundation for LIFE based on ambiguity loss decomposition and
diversity of base learners is provided.

The rest of the paper is organized as follows: In Section 2, we introduce LIFE
algorithm and theoretical rationale behind LIFE through loss decomposition.
Extensive experiments on simulated and real data are conducted in Section 3 to
test the performance of LIFE under various conditions in comparison with other
benchmark algorithms. In Section 4, we explored interpretation of model such
as main or interaction effect detection. In section 5, we discuss the model
pruning and the extension of LIFE algorithm in more depth. In Section 6, we
provide our conclusions. In section 6, we discuss the model pruning and the
extension of LIFE algorithm in more depth.

##  2 Methodology

In this section, we introduce a new proposed ensemble method called LIFE. The
LIFE algorithm is mainly used to train a wide single-hidden-layer NNs in an
iterative way. First, the general framework for LIFE including three steps is
introduced in Section 2.1. Then, details of LIFE algorithms is provided in
Section 2.2. Finally the theoretical foundation is discussed in Section 2.3.

###  2.1 General framework for LIFE

As shown in the Figure 1, the LIFE algorithm consists of three steps: data
sampling, base learner training and feature extraction, model aggregation and
pruning. Figure 1 presents several options in the white box for each step,
which allow for various combinations of those steps to achieve pre-specified
goals.

  1. 1.

The first step consists in defining subsets of data via active functions of NN
neurons, which are obtained by training multiple single-hidden-layer NNs in a
hierarchal structure as illustrated in Figure 2. The diversified base learners
can be generated by training based on these subsets for data sets. More
theoretical explanation on ensemble with diversified base learners will be
provided in Section 2.3. There are other alternative ways to generate sampling
for base learner training, e.g., bootstrapping in the traditional method or
data splitting via random projection. Through leveraging linear projections
from trained NNs, our sampling method in the supervised setting can more
effectively generate the diversity among base learners, which is demonstrated
by empirical experiments in Section 5.1.

  2. 2.

The second step consists in training base learners on different subsets of
data sampled during the first step. Various options can be considered, e.g.,
single layer NN, multiple layer NN, and regression or decision tree, etc., In
this paper, the single-hidden-layer NN is used as the base learner given its
interpretability. After estimating multiple NN base learners, all activation
functions of the neurons from NN base learners are extracted as new features.

  3. 3.

The third step consists in combining all new extracted features from different
base learners in the second step to construct the final predictive model. In
addition, new features can be pruned through regularization or other methods
to generate a more parsimonious model. We use linear model and elastic–net in
our LIFE algorithm which is simple and straightforward in this paper, but some
other more complicated model aggregating methods e.g., adaptive regression
model screen, and pruning methods, e.g., base learner selection algorithm 3,
will be discussed in Section 5.2.

Figure 1: It is the general framework of LIFE algorithm and options colored in
red in each step are used in the paper.

###  2.2 LIFE Algorithm

LIFE algorithm is an iterative process with multiple single layer NN base
learners trained in each iteration. Assume there are $J$ iterations. The first
$J-1$ iterations is used to define the data sampling through a hierarchical
structure, and the last $J$ iteration is used to build the features from
single layer NN base learners. Moreover, let $[K_{1},\cdots,K_{J}]$ denote the
collection of the number of hidden neurons for single-hidden-layer NNs in all
iterations, where $K_{j}$ is the number of hidden neural nodes for single-
hidden-layer NNs which are trained in the $j^{(th)}$ iteration. Figure 2 below
gives an illustration to LIFE framework, with $[K_{1},K_{2},K_{3}]=[3,3,2]$,
where $\hat{b}_{k}^{(j)}$’s and $\hat{w}_{k}^{(j)}$’s represent biases and
weights from single-hidden-layer NNs respectively, $\hat{\beta}_{k}$’s are
coefficients used to linearly combine all new neurons in the final step, and
the $cp$ is the cutoff point to define subsets by controlling subset size.

Figure 2: An illustration of LIFE framework with
$[K_{1},K_{2},K_{3}]=[3,3,2]$. Regions with a smiling emoji face indicate that
the subsets satisfying the sampling conditions and are used for NN models
training in the next iteration.

The LIFE framework illustration in Figure 2 can be separated into three steps
as displayed in Figure 1. The two iterations in the first step are used to
perform data sampling, in which we first fit a single-hidden-layer NN with
three hidden neurons, then define subsets by
$b_{k}^{(1)}+x_{i}^{T}w_{k}^{(1)}>cp$, where $k=1,2,3;i=1,\cdots,N$, given
estimated bias and weight in each neuron. Further, the entire neuron is
dropped and no longer be used for next iteration if the defined subset in this
neuron is either too large or too small based on pre-specified criteria. For
example, the green neuron in the middle is dropped due to small size of
subset. If the subset size is close to the full data size, the sample is
almost identical to original training data, which is not beneficial for
generating diversified samples. On the other hand, if the subset size is too
small, the sample is not representative of the original data, and the base
learner built on this sample does not have good performance on the entire
training data. We will show in Section 2.3, the diversity and accuracy of the
base learners are the two key elements for the final ensemble performance. As
the result of the first iteration in the first step, we leverage NN linear
hyperplane in the neurons to perform data partition as shown in Figure 3, of
which the idea is similar to oblique trees [14]. Plots (a), (b) and (c) in
Figure 3 have shown how data are partitioned in different ways for these three
neurons based on $w_{k}^{(1)}$’s and $b_{k}^{(1)}$’s, $k=1,2,3$, obtained from
the first iteration. The non-white regions (yellow, green or red) indicate the
subsets, where all observation satisfy $b_{k}^{(j)}+x_{i}^{T}w_{k}^{(j)}>cp$,
and will be used to fit single-hidden-layer NNs in the second iteration.
Notice that the non-white region in plot (b) is so small, which corresponds to
the dropped green neuron in Figure 2. Plot (d) in Figure 3 shows the
combination of plot (a) and (b). As we can see, the two subsets are
overlapped, which is different from the subsets defined in the traditional
regression or decision tree structure, which uses the exclusive partition.

Figure 3: That is data partition by NN linear projection after the first
iteration and observations in the colored area from (a), (b), (c) will be
selected in the each subset. The (d) shows two kept data partitions in one
plot.

In the second iteration of the first step, single-hidden-layer NNs with three
hidden neurons are fitted independently using the subsets defined from the
first iteration. After that, the entire training set is evaluated for data
sampling through NN linear projection $b_{k}^{(2)}+x_{i}^{T}w_{k}^{(2)}$,
where $k=1,\cdots,6;i=1,\cdots,N,$ and new subsets are defined satisfying
$b_{k}^{(2)}+x_{i}^{T}w_{k}^{(2)}>cp$. Note that, the new subsets are defined
on the entire training data, not the samples from the previous iteration,
which is also very different from traditional regression or decision tree
structure. Again, neurons with too small or too large subset are dropped
forever, as shown in the first node (in brown) and the last node (in pink) of
the second hidden layer from Figure 2, corresponding to (a1) and (c3) in
Figure 4. In addition, Figure 4 shows data partition obtained from the second
iteration, where plots (a1), (a2) and (a3) represent the partitions generated
from NN trained on the subset of the first node (yellow) in the first
iteration, while plots (c1), (c2) and (c3) correspond to the third node (red)
in the first iteration. Plot (d2) shows that each training data point as least
belongs to one of the six subsets in (a1-a3) and (c1-c2). We do expect all or
most data points are covered by different subsets. The reason is that the
neurons with different active regions in NN are representing different
features and patterns from the data, LIFE trains the same type of base learner
model on part of the dataset but evaluates it on the entire dataset. This
leads to small errors appearing in the region of sampled data and large errors
outside region. Therefore, data sampling with these active regions can
effectively define subsets for generating more diverse representation of data
and producing less correlated prediction errors. The sampling in this
supervised manner is better than sampling in a random way and this is further
discussed in Section 5. In practice, the first step can be reduced to one
iteration or have more than two iterations, displayed in Algorithm 1.

Figure 4: It is data partition by NN linear projection after the second
iteration and observations in the colored area from (a1-a3), (c1-c3) will be
selected in the each subset. The (d1) displays all training data points in
green, and the (d2) shows four kept data partitions in one plot.

Input: $[K_{1},\cdots,K_{J}]$, upper bound: $u$, lower bound: $l$, cutoff
point: $cp$

Output: $\hat{\beta}_{k},\hat{b}_{k}^{(J)},\hat{w}_{k}^{(J)}$,
$k=1,\cdots,m_{j}$, where $m_{j}$ is the total number of neural nodes left in
the $j^{th}$ iteration.

Data: Training set $\\{x_{i},y_{i}\\}_{i=1,\cdots,N}$

1 for  _$j=1,\cdots,J$_ do

2 if  _$j=1$_ then

3 1 .Train a single-hidden-layer NN regressor or classifier with the number of
hidden neurons equal to $K_{j}$ on training set $x,y$ by specified optimizer.

4 2\. Collect estimated parameters $\hat{b}_{k}^{(j)},\hat{w}_{k}^{(j)}$,
where $k=1,\cdots,K_{j}$.

5 3\. $m_{j}=K_{j}$

6 else

7 for  _$k=1,\cdots,m_{j-1}$_ do

8 1\. Initialization: $M=0$ , where $M$ records the number of NNs trained in
the current iteration.

9 2\. Sample observations from training set $\\{x_{i},y_{i}\\}_{i=1,\cdots,n}$
based on $\hat{b}_{k}^{(j-1)}+x_{i}^{T}\hat{w}_{k}^{(j-1)}>cp$, where
$x_{i}\subseteq\\{x_{i}\\}_{i=1,\cdots,n}$.

10 3\. Record size of sampled observations as $s$.

11 4\. if  _$l <s/N<u$ _ then

12 1\. Train a single-hidden-layer NN regressor or classifier with the number
of hidden neurons equal to $K_{j}$ on sampled subset of training set $x,y$ by
specified optimizer.

13 2\. Collect estimated parameters $\hat{b}_{q}^{(j)},\hat{w}_{q}^{(j)}$,
where $q=1,\cdots,K_{j}$.

14 3\. $M=M+1$

15

16 4\. else

17 Skip this step.

18 5\. $m_{j}=M\times K_{j}$

19

20

21From $M$ trained NNs, collect $\hat{b}_{k}^{(J)},\hat{w}_{k}^{(J)}$, where
$k=1,\cdots,m_{J}$.

22 Calculate values of each neural nodes
$\sigma(\hat{b}_{k}^{(J)}+x^{T}\hat{w}_{k}^{(J)})$, which are called new
features.

Fit a linear regression or logistic model with or without regularization using
new features as input to estimate $\hat{\beta}_{k}$, where $k=1,\cdots,m_{J}$.

Algorithm 1 LIFE

In the second step, four single-hidden-layer NNs, each with two hidden neurons
are trained as base learners independently on the subsets defined at the end
of the first step, and then the features are generated by the ReLU activation
functions. Please note that the features are evaluated on the entire training
dataset. In the third step, we combine all these new features together.Since
all features are obtained based on the entire training set, which technically
forms a design matrix with dimension $N\times m_{J}$, where $N$ is the size of
training set and $m_{J}$ is the total number of extracted features from step 2
with $J=8$ shown in Figure 2. Then a linear model is fitted on these features
and $\beta_{i},i=1,\cdots,m_{J},$ is the coefficient of each feature. In the
default setting, the linear regression or logistic model is applied to combine
neural nodes extracted from different base learners and make a final
prediction. However, there may be too many features and some of them can even
be highly correlated, leading to overfitting problem. Therefore, we need to
prune some redundant neural nodes through adding regularization or removing
some base learners. Both methods can not only prevent overfitting, but also
produce a more parsimonious NN model with fewer nodes that is beneficial for
interpretation. As a regularized regression method that linearly combines the
$L_{1}$ and $L_{2}$ penalties, the elastic net is one of the options, and it
is used in the paper for the third step of model aggregation and pruning.
LASSO and ridge regression are treated as special cases of elastic net.

By combining multiple relatively narrow but diversified single-hidden-layer
NNs with $K_{J}$ hidden neurons for each, the LIFE algorithm finally
constructs a wide single-hidden-layer NN with $m_{J}$ hidden neurons, which is
$8$ in Figure reffig:tree. This trained process is completely different from
the traditional methods, with which a NN is optimized as a whole by stochastic
gradient based optimizers. In most cases, It is numerically difficult and
computational expensive to train a single layer NN and achieve good
performance. LIFE can help overcome this difficulty and achieve decent
performance. In addition, LIFE can leverage parallel computation to
significantly reduce the training time. In the end, we provide the detailed
pseudo-code of LIFE in $J$ iterations, described in Algorithm 1, where $m_{j}$
indicates the number of remaining neural nodes after the $j^{th}$ iteration in
the first step, where $j=1,\cdots,J$. Both $u$ and $l$ are maximal and minimal
proportions of training set size, which provides the upper and lower bound for
subset size, respectively. The neuron will be dropped if the proportion of
subset size is beyond the range.

###  2.3 Theoretical Foundation

The ensemble method, such as stacking [27], bagging [4], boosting [9] or
Bayesian model averaging [20, 21] is composed of a multiple independently
trained regressors or classifiers whose predictions are combined to make final
predictions. Empirically, ensembles tend to yield better results than a single
model when there is a significant diversity among the models [18]. For the
past few decades, many studies have been focusing on accuracy and diversity of
ensemble methods in either regression [6] or classification case [1, 11, 7].
(Krogh and Vedelsby (1994) [15]) proposed ambiguity decomposition and a
computable approach to minimize the quadratic error of the ensemble estimator,
while (Ueda and Nakano (1996) [25]) derived a general expression of bias-
variance-covariance decomposition. (Brown et al. (2005) [5]) and (Hansen
(2000) [12]) investigated the connections between ambiguity decomposition and
bias-variance-covariance, and have shown they are identical. Based on
ambiguity decomposition, we establish the theoretical foundation for LIFE and
extend loss decomposition for both regression or classification setting, which
will be discussed in Section 2.3.1 and 2.3.2.

####  2.3.1 Connection to Stacking

Stacking is a type of ensemble method, by which a final model is trained from
the combined predictions of another models. In stacking, the predictions from
different machine learning models are used as new inputs and are combined to
generate a new set of predictions. Those predictions can be used on additional
layers, or the process can stop here with a final result. One important
assumption behind stacking is that different base learners can produce weakly
correlated prediction errors that are complementary. If we use weighted
averages, we might believe that some of the base learner are better or more
accurate and can be assigned higher weights. In the framework of stacking, an
even better approach might be to estimate these weights more intelligently by
using another layer of the learning algorithm, such as the linear model.

Aspect |  Stacking |  LIFE
---|---|---
Base Learner |  different models or same model with different settings |  same model with same setting trained on different subsets of data
Model averaging |  linear or nonlinear combination of prediction from base learners |  linear or nonlinear combination of features flatten from base learners
Table 1: Difference between LIFE and Stacking

(a)

(b)

Figure 5: Relationship between LIFE and Stacking (Regression)

(a)

(b)

Figure 6: Relationship between LIFE and Stacking (Classification)

Some major differences between stacking method and LIFE algorithm as shown in
Table 1. Without neural nodes flattening, LIFE is very similar to
stacking.Despite the differences between the two methods, the minimization of
loss function of LIFE is approximately equivalent to the minimization of loss
function for the two-stage stacking method, which firstly fits multiple base
learners (NN trained on different subsets sampling by LIFE algorithm) and use
predictions from base learners as input to train a model averaging model.
Figure 5 and 6 indicate a strong linear relationship between LIFE and the
stacking method with different hyper-parameter setup in terms of MSE loss or
cross entropy loss in both regression and classification cases. This linear
relationship has be verified by both simulated data (MIM) and real data
(California Housing for regression and Gamma Telescope for classification).
Due to the BLUE (Best linear unbiased prediction) property of OLS estimator in
linear regression, the joint estimation of the coefficients of all the
combined features in the three step makes LIFE always outperform two-stage
stack ensemble method with the same setting. This can be verified in Figure 5
that all the points are below the green diagonal line. For classification
case, LIFE also performs better than stacking method with smaller minimum loss
as shown in the white box of Figure 6.

####  2.3.2 Loss Function Decomposition

(Krogh and Vedelsby (1994) [15]) proposed ambiguity decomposition for
quadratic error of the ensemble estimator which is the sum of the quadratic
loss of individual base learners and the ambiguity measure for diversity. We
extend the ambiguity decomposition to both mean square error and cross entropy
error via Taylor expansion, where the two loss functions corresponds to
regression and classification, respectively. Let an ensemble model with M base
learners be expressed as $f_{ens}=\sum_{j=1}^{M}\beta_{j}f^{(j)}$, where
$\sum_{j=1}^{M}\beta_{j}=1$ and $\beta_{j}\geq 0$.

For any loss function that is twice differentiable, we can expand the loss
function of $j^{th}$ base learner around output of an ensemble model based on
Taylor’s theorem with Peano form of the remainder as follows:

$\begin{gathered}l(y,f^{(j)})=l(y,f_{ens})+l^{\prime}(y,f_{ens})(f^{(j)}-f_{ens})+\frac{1}{2}l^{\prime\prime}(y,f^{(j)\star})(f^{(j)}-f_{ens})^{2},\end{gathered}$
(1)

where the value of $f^{(j)\star}$ is between $f_{ens}$ and $f^{(j)}$.
Multiplying both sides of equation (1) by $w_{j}$ and taking a sum yield:

$\begin{gathered}\sum_{j=1}^{M}\beta_{j}l(y,f^{(j)})=\sum_{j=1}^{M}\beta_{j}l(y,f_{ens})+\sum_{j=1}^{M}\beta_{j}l^{\prime}(y,f_{ens})(f^{(j)}-f_{ens})\\\
+\frac{1}{2}\sum_{j=1}^{M}\beta_{j}l^{\prime\prime}(y,f^{(j)\star})(f^{(j)}-f_{ens})^{2}.\end{gathered}$
(2)

The second term on the right side of (2) is expressed by:

$\begin{gathered}\sum_{j=1}^{M}\beta_{j}l^{\prime}(y,f_{ens})(f^{(j)}-f_{ens})=l^{\prime}(y,f_{ens})\\{\sum_{j=1}^{M}\beta_{j}f^{(j)}-f_{ens}\sum_{j=1}^{M}\beta_{j}\\}\\\
=l^{\prime}(y,f_{ens})\\{f_{ens}-f_{ens}\\}=0.\end{gathered}$ (3)

Since this term is zero, the loss function $l(y,f_{ens})$ of the ensemble can
be decomposed into:

$\begin{gathered}l(y,f_{ens})=\sum_{j=1}^{M}\beta_{j}l(y,f^{(j)})-\frac{1}{2}\sum_{j=1}^{M}\beta_{j}l^{\prime\prime}(y,f^{(j)\star})(f_{ens}-f^{(j)})^{2}.\end{gathered}$
(4)

In regression case, let $f_{i}=\sum_{j=1}^{M}\beta_{j}f_{i}^{(j)}$ be
individual predicted value of an ensemble model for $i^{th}$ observation,
where $M$ is the number of base learners, $f_{i}^{(j)}$ represents individual
predicted value of $j^{th}$ base learner for $i^{th}$ observation and
$\beta_{j}$ denotes regression coefficient for $j^{th}$ base learner.
Basically, the mean squared error (MSE) is commonly used loss function
$l(y,f)=\frac{1}{N}\sum_{i}^{N}(y_{i}-f_{i})^{2}$ for regression problems,
where $N$ is the total number of observations in the entire dataset. Based on
equation (5), MSE of an ensemble model can be written in terms of the
ambiguity decomposition given $x_{i},i=1,\cdots,N$:

$\begin{gathered}MSE=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-f_{i})^{2}=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\sum_{j=1}^{M}\beta_{j}f_{i}^{(j)})^{2}\\\
=\underbrace{\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}(y_{i}-f_{i}^{(j)})^{2}}_{\text{accuracy}}-\underbrace{\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}(f_{i}-f_{i}^{(j)})^{2}}_{\text{diversity}}.\end{gathered}$
(5)

On the right-hand side of equation (5), the first term of this decomposition
is to measure average prediction accuracy of base learners, while the second
term is called ambiguity (hence the name of the decomposition) and can be
easily interpreted in terms of diversity between individual base learners.
Unlike the bias-variance-covariance decomposition, the ambiguity decomposition
highlights a trade-off between the average accuracy of base learners, and
their deviation from the ensemble output.

Regarding LIFE, the base learner is a single-hidden-layer neural network
trained on a subset of all observations. A stronger base learner indicates a
better performance of model, which is reflected by first term of equation (5).
If the subset size is small, the base learner is also weak, which deteriorates
performance. Thus, a lower bound is set up. The power of LIFE framework comes
from second term diversity, which is due to data sampling during iterations.
Creating different subsets through sampling allows the model to be trained on
different aspects of data, which produces diversity deliberately without
resorting to other machine learning algorithms. In general, the more diverse
the subsets, the better the predictive performance of LIFE. Hence, the upper
bound is necessary to ensure diversity of subset since subset contains almost
all observations and its size is very large, making subsets loss diversity.
Another parameter cutoff point is set up to balance accuracy and diversity as
well.

For binary classification purposes, let
$f_{i}=\sum_{j=1}^{M}\beta_{j}f_{i}^{(j)}$ be individual predicted probability
of an ensemble model for $i^{th}$ observation, which is weighted average of
predicted probability or log-odds of base learner $f_{i}^{(j)}$, where
$\sum_{j=1}^{M}\beta_{j}=1$ and $\beta_{j}\geq 0$. The cross-entropy loss is
widely used for classification and it can be written as follows for single
observation:

$\begin{gathered}l(y_{i},f_{i})=-y_{i}log(f_{i})-(1-y_{i})log(1-f_{i}).\end{gathered}$
(6)

By plugging the loss function (14) into equation (4), we can write average
cross entropy loss of an ensemble method on training set
$\\{x_{i},y_{i}\\}_{i=1,\cdots,N}$ in the probability space as follows:

$\begin{gathered}\sum_{i=1}^{N}[-y_{i}log(f_{i})-(1-y_{i})log(1-f_{i})]\\\
=\underbrace{\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}[-y_{i}log(f_{i}^{(j)})-(1-y_{i})log(1-f_{i}^{(j)})]}_{\text{accuracy}}\\\
-\underbrace{\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}\\{\frac{y_{i}-2f_{i}^{(j)\star}y_{i}+(f_{i}^{(j)\star})^{2}}{[f_{i}^{(j)\star}(1-f_{i}^{(j)\star})]^{2}}\\}(f_{i}-f_{i}^{(j)})^{2}}_{\text{diversity}},\end{gathered}$
(7)

(a)

(b)

(c)

(d)

Figure 7: Loss decomposition for regression

(a)

(b)

(c)

(d)

Figure 8: Loss decomposition for classification

where $f_{i}^{(j)\star}$ takes value between $f_{i}^{(j)}$ and $f_{i}$. The
term $(f_{i}-f_{i}^{(j)})^{2}$ is a measure of the differences in value
between base learner and the ensemble. The cross-entropy loss and its
decomposition in the log-odds space is provided in the Appendix Loss
Decomposition in the Log-odds Space. Unlike diversity term in the regression
case, the second term (diversity) in the right-hand side of equation (7) also
includes the true class label $y_{i}$ and parameter with unknown value
$f_{i}^{(j)\star}$. However, the interpretation of decomposition is also
clear. It shows that a lower average accuracy of individual base learner can
be compensated by a higher disagreement with the ensemble, scaled by
$\frac{1}{(f_{i}^{(j)\star})^{2}}$ if $y_{i}=1$ or
$\frac{1}{(1-f_{i}^{(j)\star})^{2}}$ if $y_{i}=0$ in the probability space.
Since $\frac{\beta_{j}}{(f_{i}^{(j)\star})^{2}}$ or
$\frac{\beta_{j}}{(1-f_{i}^{(j)\star})^{2}}$ is positive, the more deviance of
predicted probability between base learner and an ensemble model implies more
diversity.

We implement the loss decomposition on simulated data (MIM) and real data
(California Housing for regression and Gamma Telescope for classification) to
LIFE without neural nodes flattening. It ensembles the predictions from
single-hidden-layer NN base learner directly which is the two-stage stacking
model averaging method discussed above. As subset size in the sampling step of
LIFE impact the strength of diversity, we explore the overall loss, the
weighted sum individual base learner accuracy and the ambiguity measure
against the average subset size over all the single-hidden-layer NNs in the
first step of LIFE. Here we vary the subset size by controlling the cutoff
point $cp$ for linear project. Figures 7 and 10 show the relationship between
average subset size (in terms of the proportion of original train data size)
and loss for accuracy, loss for diversity and MSE loss in both regression and
classification cases. In plot (a) and (b), the blue curves indicating the
ambiguity diversity measure have a decreasing trend when the subset size
increases. This is consistent with the intuition that the larger overlapping
the subsets are, the less diverse the base learners are. On the other hand,
the accuracy of the individual base learners is higher when subset size is
larger, as the sample is more representative of the whole training set.
Similarly, training on smaller subsets has stronger diversity but leads to
lower accuracy. Plot (c) and (d) illustrate the tradeoff between average
accuracy and diversity to minimize total loss, where the dash line shows the
optimal subset size achieving the minimum loss. Combining the results from
Section 2.3.1 and Section 2.3.2, we conclude that competitive predictive
performance of LIFE benefits from diversity due to data sampling in the first
step and the feature “flattening” and joint estimation in the third step.

##  3 Empirical Experiment

In this section, we conduct multiple empirical experiments via both simulated
and real data for regression and classification cases to confirm the
competitive performance of LIFE. We have generated multiple datasets with
normal distribution and heavy-tailed predictor distribution (Laplace
distribution), as well as different function forms to analyze predictive
performance and computational efficiency. Other benchmark models including
single-hidden-layer NN trained by different optimizer (Local linear
approximation and Adam algorithm), and other machine learning algorithms
including multi-layer FFNN, Xgboost, and Random Forest are tested on the same
data for comparison after extensive hyper-parameter tuning. Local Linear
Approximation (LLA) algorithm is a recently proposed method to estimate the
weights and biases of single-hidden-layer NN by iterative linear regression
and linear approximation of the ReLU activation function [29]. The LLA
algorithm is distinguished from existing gradient descent algorithms in that
it utilizes the Hessian matrix in the same spirit of Fisher scoring algorithm
for nonlinear regression models with normal error. The outline of the LLA
algorithm is included in the Appendix Local Linear Approximation (LLA)
Algorithm.

###  3.1 Simulated Data

####  3.1.1 Regression

For the regression scenario, there are three different function forms
including Generalized Additive Model (GAM), Additive Index Model (AIM) and
Multiple Index Model (MIM), which are expressed, as follows:

$\begin{gathered}GAM:y_{i}=\beta_{1}x_{1i}+\beta_{2}\sqrt{|x_{2i}|}+\beta_{3}|x_{3i}|+\beta_{4}exp(x_{4i})+\beta_{5}log(|x_{5i}|)+\beta_{6}max(1,x_{6i})+\epsilon_{i},\\\
\beta=\\{\beta_{1},\cdots,\beta_{6}\\}=\\{1.5,\sqrt{5},2,4e^{(\frac{-1.5}{7})},4log(1.5),-4\\},\epsilon_{i}\sim
N(0,1),i=1,\dots,N,\end{gathered}$ (8)

$\begin{gathered}AIM:y_{i}=2log(|\beta_{1}x_{1i}+\cdots+\beta_{4}x_{4i}|)+exp(\frac{\beta_{3}x_{3i}+\cdots+\beta_{6}x_{6i}}{9})\\\
+max(0,\beta_{5}x_{5i}+\beta_{6}x_{6i})+\epsilon_{i},\\\
\beta=\\{\beta_{1},\cdots,\beta_{6}\\}=\\{3,-2.5,2,-1.5,1.5,-1\\},\end{gathered}$
(9)

$\begin{gathered}MIM:y_{i}=exp(\beta_{1}x_{1i}+\beta_{2}x_{2i})\beta_{3}x_{3i}+\frac{\beta_{4}x_{4i}}{1+\beta_{5}|x_{5i}|}+max(2,\beta_{6}x_{6i})+\epsilon_{i},\\\
\beta=\\{\beta_{1},\cdots,\beta_{6}\\}=\\{0.03,-0.025,1,-3,1.5,-2\\},\end{gathered}$
(10)

where $N=20k$ and all predictors $\\{x_{ji}\\}_{j=1,\cdots,6;i=1,\cdots,N}$
are drawn from Normal or Laplace distribution. For regression, the
experimental results show mean and standard deviation of $RMSE$, $R^{2}$ and
training time $T$ over five replications in Tables 2 and 3, while log-loss and
AUC are used as a performance metric in the classification case as shown in
Tables 4 and 5. Bayesian Optimization allows us to jointly tune more
parameters with fewer experiments and find better values, so we implement it
to perform extensive hyper–parameter tuning on all the algorithm. The
important hyper-parameters for LIFE include the number of iterations, the
number of neurons in each iteration, upper and lower bound. We marked optimal
results that have won the campaign in red.

Model | Metric | Oracle | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---|---
|  |  | 1.046 | 1.062 | 1.063 | 1.149 | 1.104 | 1.065 | 1.906
| $RMSE$ | $1.000$ | (0.012) | (0.040) | (0.016) | (0.023) | (0.025) | (0.020) | (0.037)
GAM |  |  | 0.962 | 0.961 | 0.961 | 0.954 | 0.957 | 0.960 | 0.873
| $R^{2}$ | $0.965$ | (0.001) | (0.003) | (0.001) | (0.002) | (0.002) | (0.002) | (0.005)
| $T$ | N/A | 75s | 120s | 46s | 55s | 81s | 8s | 44s
|  |  | 1.077 | 1.092 | 1.153 | 1.202 | 1.159 | 2.606 | 3.285
| $RMSE$ | $1.000$ | (0.018) | (0.032) | (0.029) | (0.040) | (0.033) | (0.173) | (0.161)
AIM |  |  | 0.981 | 0.980 | 0.978 | 0.975 | 0.978 | 0.886 | 0.828
| $R^{2}$ | $0.984$ | (0.001) | (0.001) | (0.001) | (0.004) | (0.001) | (0.010) | (0.010)
| $T$ | N/A | 7s | 51s | 17s | 40s | 46s | 27s | 48s
|  |  | 1.028 | 1.043 | 1.036 | 1.047 | 1.037 | 1.057 | 1.139
| $RMSE$ | $1.000$ | (0.005) | (0.003) | (0.006) | (0.004) | (0.009) | (0.004) | (0.011)
MIM |  |  | 0.933 | 0.931 | 0.931 | 0.930 | 0.931 | 0.929 | 0.918
| $R^{2}$ | $0.937$ | (0.003) | (0.003) | (0.003) | (0.003) | (0.007) | (0.003) | (0.008)
| $T$ | N/A | 29s | 313s | 81s | 45s | 51s | 14s | 44s

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. The figures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 2: Regression on Simulated Data (Normal Distribution)

Model | Metric | Oracle | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---|---
|  |  | 1.166 | 1.258 | 1.427 | 1.629 | 1.460 | 1.439 | 3.242
| $RMSE$ | $1.000$ | (0.079) | (0.091) | (0.192) | (0.100) | (0.205) | (0.256) | (0.226)
GAM |  |  | 0.989 | 0.987 | 0.984 | 0.979 | 0.983 | 0.983 | 0.920
| $R^{2}$ | $0.992$ | (0.001) | (0.001) | (0.003) | (0.002) | (0.004) | (0.005) | (0.005)
| $T$ | N/A | 48s | 84s | 19s | 96s | 75s | 5s | 20s
|  |  | 1.057 | 1.084 | 1.122 | 1.192 | 1.166 | 2.031 | 2.491
| $RMSE$ | $1.000$ | (0.026) | (0.040) | (0.055) | (0.084) | (0.077) | (0.099) | (0.204)
AIM |  |  | 0.966 | 0.964 | 0.962 | 0.957 | 0.959 | 0.876 | 0.804
| $R^{2}$ | $0.968$ | (0.004) | (0.005) | (0.003) | (0.004) | (0.004) | (0.012) | (0.035)
| $T$ | N/A | 10s | 61s | 22s | 59s | 26s | 8s | 56s
|  |  | 1.060 | 1.110 | 1.060 | 1.096 | 1.088 | 1.099 | 1.555
| $RMSE$ | $1.000$ | (0.017) | (0.065) | (0.018) | (0.026) | (0.018) | (0.013) | (0.081)
MIM |  |  | 0.965 | 0.961 | 0.965 | 0.962 | 0.963 | 0.962 | 0.925
| $R^{2}$ | $0.969$ | (0.001) | (0.005) | (0.001) | (0.002) | (0.001) | (0.001) | (0.008)
| $T$ | N/A | 20s | 195s | 39s | 42s | 32s | 6s | 43s

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. The figures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 3: Regression on Simulated Data (Laplace Distribution)

As illustrated in Tables 2 and 3, the result is predictable regardless of the
distribution predictors drawn. LIFE algorithm with LLA optimizer achieves
higher accuracy among all methods in terms of predictive performance on the
test set. The values of two metrics from LIFE (LLA) are close to oracle
values, which implies LIFE performs well in the data with a smoothing response
surface. If we compare results from one hidden-layer FFNNs trained by LIFE
algorithm with either LLA or Adam base learners and non-ensemble algorithms of
LLA or Adam, LIFE always outperforms the relevant optimization methods used to
train single hidden-layer FFNN as a whole due to the generated diversity of
data sampling. In addition, the performance of LIFE also depends on the
strength of individual NN base learner, which can be easily spotted in the
Tables 2 and 3 that LIFE (LLA) outperforms LIFE(Adam). From the perspective of
computational efficiency, LIFE algorithm also shows some advantages over other
single-hidden-layer NN training algorithms. In general, LIFE algorithm can not
only boost predictive performance of one-hidden layer NN, but also speed up
training, especially with respect to wide NN with large hidden layer
dimension.

####  3.1.2 Classification

For classification case, the functional forms in simulation setup are similar
to the ones in regression case except that the coefficients are a little bit
different. Detailed information on formulas can be found in Appendix
Functional Forms for Classification Case in Simulation Study. Similar to the
setup in regression case, we choose $N=20k$ and all predictors are drawn from
either Normal or Laplace distribution. The response variable is sampled from
Bernoulli distribution with probability calculated using the logit link
function.

Model | Metric | Oracle | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---|---
|  | $0.984$ | 0.974 | 0.968 | 0.973 | 0.945 | 0.973 | 0.975 | 0.964
| $AUC$ | $(0.002)$ | (0.002) | (0.003) | (0.002) | (0.003) | (0.002) | (0.002) | (0.002)
GAM |  | $0.113$ | 0.146 | 0.158 | 0.148 | 0.247 | 0.152 | 0.141 | 0.166
| $logloss$ | $(0.005)$ | (0.004) | (0.010) | (0.005) | (0.029) | (0.005) | (0.006) | 0.006
|  | $0.836$ | 0.751 | 0.535 | 0.752 | 0.742 | 0.752 | 0.687 | 0.715
| $AUC$ | $(0.008)$ | (0.010) | (0.037) | (0.010) | (0.010) | (0.010) | (0.010) | (0.011)
AIM |  | $0.359$ | 0.422 | 0.485 | 0.428 | 0.448 | 0.431 | 0.457 | 0.447
| $logloss$ | $(0.005)$ | (0.003) | (0.005) | (0.006) | (0.010) | (0.004) | (0.001) | (0.002)
|  | $0.888$ | 0.839 | 0.838 | 0.838 | 0.830 | 0.839 | 0.837 | 0.834
| $AUC$ | $(0.005)$ | (0.008) | (0.007) | (0.008) | (0.008) | (0.008) | (0.008) | (0.007)
MIM |  | $0.331$ | 0.386 | 0.388 | 0.387 | 0.399 | 0.387 | 0.388 | 0.391
| $logloss$ | $(0.010)$ | (0.010) | (0.009) | (0.010) | (0.011) | (0.009) | (0.010) | (0.009)

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. The figures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 4: Classification on Simulated Data (Normal Distribution)

Model | Metric | Oracle | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---|---
|  | $0.996$ | 0.990 | 0.984 | 0.990 | 0.857 | 0.976 | 0.991 | 0.979
| $AUC$ | $(0.000)$ | (0.001) | (0.003) | (0.001) | (0.012) | (0.002) | (0.001) | (0.002)
GAM |  | $0.051$ | 0.081 | 0.100 | 0.081 | 0.252 | 0.133 | 0.079 | 0.120
| $logloss$ | $(0.002)$ | (0.003) | (0.006) | (0.003) | (0.012) | (0.004) | (0.004) | (0.003)
|  | $0.866$ | 0.802 | 0.797 | 0.802 | 0.795 | 0.802 | 0.737 | 0.743
| $AUC$ | $(0.012)$ | (0.018) | (0.020) | (0.018) | (0.021) | (0.017) | (0.014) | (0.015)
AIM |  | $0.230$ | 0.270 | 0.274 | 0.271 | 0.281 | 0.273 | 0.298 | 0.296
| $logloss$ | $(0.006)$ | (0.008) | (0.012) | (0.008) | (0.007) | (0.007) | (0.004) | (0.004)
|  | $0.958$ | 0.939 | 0.938 | 0.940 | 0.920 | 0.937 | 0.938 | 0.934
| $AUC$ | $(0.003)$ | (0.004) | (0.003) | (0.004) | (0.003) | (0.002) | (0.004) | (0.004)
MIM |  | $0.225$ | 0.268 | 0.272 | 0.268 | 0.421 | 0.276 | 0.271 | 0.279
| $logloss$ | $(0.011)$ | (0.012) | (0.011) | (0.012) | (0.075) | (0.010) | (0.011) | (0.010)

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. The figures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 5: Classification on Simulated Data (Laplace Distribution)

Table 4 and 5 show the simulation results from binary scenario, where data are
generated from Normal distribution and Laplace distribution, respectively.
Similar to the results in regression case, LIFE(LLA) has won the campaign in
four out of six functional forms. LIFE(Adam) also performs pretty well
especially in Table 5 for Laplace Distribution. For data drawn from Normal
Distribution shown in Table 4, LIFE(Adam) is also quite close to the optimal
result. Furthermore, there is a strong evidence to show LIFE algorithm does
improve the performance of base learners with larger AUC, smaller log-loss and
smaller standard errors of both metrics. Even if the base learner is not
strong enough, like Adam for GAM and MIM in Laplace Distribution (Table 5),
with which AUC or log-loss or both has large standard error, the ensemble
approach in LIFE(Adam) can dramatically reduce the variance. XGboost ranks at
top for GAM in both distributions, however, the differences between LIFE and
XGboost are negligible with only 0.1% of difference in AUC and 3% of
difference in log-loss.

###  3.2 Real Data

Besides implementing LIFE algorithm on simulated data, we also tested it on 7
public datasets for regression and eight datasets for classification and
compared it with other benchmark models including single-hidden-layer FFNN and
Xgboost. All data sets are split into $80\%$ training data and $20\%$ testing
data with 10 different random seeds, which yield results over $10$
replications. For all the datasets, we transformed categorical variables into
dummy variables and standardized the continuous variables, so that the mean
and the variance of each continuous variable are equal to $0$ and $1$,
respectively. A detailed description of all datasets and corresponding data
preprocessing steps are outlined in the Appendix Real Datasets.

####  3.2.1 Regression

The experiment results averaged over ten replications are reported in Table 6,
including root mean squared error ($RMSE$), R-squared ($R^{2}$) and training
time ($T$).

Model | Metric | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---
|  | 2.112 | 2.162 | 2.137 | 2.179 | 2.195 | 2.169 | 2.170
| $RMSE$ | (0.077) | (0.074) | (0.071) | (0.108) | (0.133) | (0.068) | (0.071)
Abalone |  | 0.570 | 0.549 | 0.560 | 0.543 | 0.534 | 0.546 | 0.546
| $R^{2}$ | (0.031) | (0.031) | (0.029) | (0.047) | (0.057) | (0.028) | (0.029)
| $T$ | 2s | 289s | 1s | 35s | 17s | 1s | 1s
|  | 0.225 | 0.288 | 0.276 | 0.293 | 0.270 | 0.238 | 0.343
| $RMSE$ | (0.020) | (0.021) | (0.026) | (0.018) | (0.034) | (0.016) | (0.023)
Airfoil |  | 0.949 | 0.917 | 0.923 | 0.913 | 0.925 | 0.943 | 0.881
| $R^{2}$ | (0.009) | (0.012) | (0.014) | (0.011) | (0.019) | (0.008) | (0.016)
| $T$ | 1s | 14s | 1s | 8s | 1s | 1s | 1s
|  | 1.193 | 1.220 | 1.210 | 1.223 | 1.255 | 1.214 | 1.204
| $RMSE$ | (0.066) | (0.045) | (0.068) | (0.096) | (0.125) | (0.101) | (0.082)
Aquatic Toxicity |  | 0.484 | 0.461 | 0.469 | 0.455 | 0.425 | 0.464 | 0.473
| $R^{2}$ | (0.058) | (0.039) | (0.060) | (0.088) | (0.117) | (0.089) | (0.071)
| $T$ | 1s | 41s | 1s | 2s | 1s | 1s | 1s
|  | 41.98 | 47.66 | 45.34 | 50.74 | 45.99 | 41.71 | 66.38
| $RMSE$ | (1.299) | (0.946) | (1.103) | (1.601) | (1.737) | (1.184) | (1.257)
Bike Sharing |  | 0.946 | 0.930 | 0.937 | 0.922 | 0.935 | 0.947 | 0.866
| $R^{2}$ | (0.003) | (0.003) | (0.003) | (0.005) | (0.005) | (0.003) | (0.005)
| $T$ | 85s | 2194s | 81s | 266s | 114 | 13s | 12s
|  | 0.500 | 0.527 | 0.521 | 0.535 | 0.519 | 0.479 | 0.486
| $RMSE$ | (0.012) | (0.011) | (0.009) | (0.008) | (0.011) | (0.007) | (0.007)
California Housing |  | 0.812 | 0.791 | 0.796 | 0.784 | 0.797 | 0.825 | 0.822
| $R^{2}$ | (0.009) | (0.008) | (0.007) | (0.007) | (0.008) | (0.005) | (0.005)
| $T$ | 13s | 108s | 8s | 74s | 40s | 12s | 40s
|  | 3.865 | 4.111 | 4.081 | 4.133 | 4.067 | 3.786 | 3.951
| $RMSE$ | (0.034) | (0.049) | (0.047) | (0.071) | (0.037) | (0.022) | (0.016)
CASP |  | 0.601 | 0.548 | 0.555 | 0.543 | 0.558 | 0.617 | 0.583
| $R^{2}$ | (0.007) | (0.011) | (0.010) | (0.016) | (0.008) | (0.004) | (0.004)
| $T$ | 216s | 2881s | 54s | 339s | 83s | 60s | 56s
|  | 0.007 | 0.011 | 0.008 | 0.009 | 0.009 | 0.011 | 0.013
| $RMSE$ | (0.000) | (0.001) | (0.000) | (0.001) | (0.001) | (0.000) | (0.000)
Electrical Grid |  | 0.959 | 0.908 | 0.949 | 0.944 | 0.945 | 0.908 | 0.877
| $R^{2}$ | (0.003) | (0.016) | (0.003) | (0.009) | (0.011) | (0.003) | (0.004)
| $T$ | 25s | 27s | 4s | 27s | 17s | 16s | 18s

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. Thefigures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 6: Regression on Real Data

As observed in Table 6, LIFE (LLA) is ranked as the best algorithm in the four
datasets. For the remaining three datasets, LIFE(LLA) is still the second or
third best algorithm among all models with a close or slightly worse
predictive performance than optimal one (Xgboost or Random Forest), which
implies that the LIFE algorithm is competitive with other state-of-art machine
learning algorithms. In addition, there is an average $4.6\%$ or $1.8\%$
improvement in R-square of all real datasets when single-hidden-layer NN is
trained by LIFE (LLA or Adam) instead of other optimization methods (LLA or
Adam), which is consistent with the conclusion made from experiment in the
simulated data. It is also worth mentioning that computation efficiency of NN
training has been significantly boosted for almost all dataset via LIFE
compared with traditional NN training methods. In particular, if we take a
look at the largest real dataset CASP, training time of NN via LIFE reduces to
$216$ seconds from $2881$ seconds or to $54$ seconds from $339$ seconds when
we use LLA as optimizer or Adam respectively, which is almost more than six
times faster. Although tree-based ensemble methods such as Xgboost and Random
Forest show strong predictive power in some datasets, they are still black-box
models, and they are hard to interpret. The biggest advantage of our proposed
algorithm LIFE is that it preserves the interpretability of model, which is
still single-hidden-layer NN with very strong predictive performance and
boosted computation efficiency.

####  3.2.2 Classification

Data | Metric | LIFE(LLA) | LLA | LIFE(Adam) | Adam | FFNN | Xgboost | RF
---|---|---|---|---|---|---|---|---
|  | 0.796 | 0.790 | 0.796 | 0.793 | 0.794 | 0.799 | 0.800
| $AUC$ | (0.006) | (0.007) | (0.006) | (0.006) | (0.006) | (0.008) | (0.005)
Bank Marketing |  | 0.288 | 0.292 | 0.289 | 0.294 | 0.288 | 0.288 | 0.286
| $logloss$ | (0.003) | (0.003) | (0.004) | (0.003) | (0.003) | (0.003) | (0.003)
|  | 0.997 | 0.991 | 0.996 | 0.994 | 0.996 | 0.992 | 0.993
| $AUC$ | (0.002) | (0.009) | (0.004) | (0.005) | (0.004) | (0.006) | (0.006)
Breast Cancer |  | 0.096 | 0.143 | 0.084 | 0.340 | 0.099 | 0.114 | 0.090
Wisconsin | $logloss$ | (0.015) | (0.088) | (0.027) | (0.207) | (0.014) | (0.042) | (0.019)
|  | 0.912 | 0.904 | 0.914 | 0.891 | 0.912 | 0.911 | 0.911
| $AUC$ | (0.001) | (0.002) | (0.001) | (0.003) | (0.001) | (0.001) | (0.001)
Higgs Boson |  | 0.353 | 0.368 | 0.348 | 0.399 | 0.353 | 0.355 | 0.357
| $logloss$ | (0.001) | (0.004) | (0.001) | (0.005) | (0.001) | (0.001) | (0.001)
|  | 0.853 | 0.841 | 0.860 | 0.848 | 0.858 | 0.858 | 0.854
| $AUC$ | (0.005) | (0.003) | (0.003) | (0.006) | (0.003) | (0.003) | (0.003)
Home Lending |  | 0.046 | 0.047 | 0.046 | 0.047 | 0.046 | 0.046 | 0.046
| $logloss$ | (0.001) | (0.001) | (0.001) | (0.001) | (0.001) | (0.001) | (0.001)
|  | 0.943 | 0.928 | 0.938 | 0.909 | 0.938 | 0.936 | 0.937
| $AUC$ | (0.002) | (0.003) | (0.003) | (0.008) | (0.003) | (0.002) | (0.003)
MAGIC Gamma |  | 0.275 | 0.311 | 0.288 | 0.371 | 0.287 | 0.291 | 0.296
Telescope | $logloss$ | (0.007) | (0.007) | (0.007) | (0.019) | (0.006) | (0.005) | (0.005)
|  | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
| $AUC$ | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000)
Mushroom |  | 0.001 | 0.000 | 0.001 | 0.024 | 0.004 | 0.004 | 0.005
| $logloss$ | (0.000) | (0.000) | (0.001) | (0.036) | (0.003) | (0.001) | (0.000)
|  | 0.998 | 0.997 | 0.998 | 0.995 | 0.997 | 0.998 | 0.995
| $AUC$ | (0.001) | (0.001) | (0.001) | (0.002) | (0.001) | (0.000) | (0.001)
Ringnorm |  | 0.049 | 0.086 | 0.051 | 0.147 | 0.075 | 0.057 | 0.162
| $logloss$ | (0.006) | (0.020) | (0.007) | (0.009) | (0.015) | (0.004) | (0.006)

  * •

Note: On the left side of vertical line, all columns represent just single-
hidden-layer NNs optimized by LIFE (LLA), LLA, LIFE (Adam) and Adam, while
there are three state-of-art machine learning methods including FFNN, Xgboost,
and Random Forest (RF) on the right side. FFNN is multi-hidden-layer feed
forward NN, where range for number of hidden layers is from two to four. In
addition, Xgboost and RF are two tree-based ensemble methods. The figures
inside parenthesis indicate the standard deviation of metrics, and time
represents training time of one replication. The numbers colored in red
represents the optimal results of this metric.

Table 7: Classification on Real Data

In the classification case, original LLA algorithm is not stable, since it
involves matrix inversion. We added a ridge parameter into the matrix
inversion and treat it as a hyper-parameter in LLA algorithm. Experiments have
shown adding ridge parameter in LLA can give better and more stable prediction
than not adding ridge parameter. After testing LLA and LIFE(LLA) with or
without ridge parameter, we further choose the best one for the performance.

Table 7 presents similar patterns in real data analyses as in simulation
studies with LIFE(LLA) and LIFE(Adam) taking turns to occupy the dominant
position for most of the data sets. LIFE performs much better than Xgboost and
Random Forest (RF) in most experiments. For example, with Breast Cancer
Wisconsin data, log-loss in LIFE(Adam) is 26.3% lower than that in Xgboost,
and with MAGIC Gamma Telescope data, log-loss has dropped by 7% from Random
Forest to LIFE(LLA). The performance of some datasets, such as Bank Marketing
data, where LIFE cannot outperform Xgboost or RF, however, the performance is
competitive, with only 0.5% and 0.7% of difference in AUC and log-loss,
respectively. Another aspect worth mentioning is that Higgs Boson data
contains quite a few highly correlated variables. The results show that LIFE
algorithm outperforms all the rest of models, which indicates LIFE really does
an excellent job in predicting on highly correlated structures.

Data | Metric | ResNet18 | ResNet18_LIFE | ResNet18_Xgboost | ResNet18_Adam
---|---|---|---|---|---
|  | 0.934 | 0.947 | 0.944 | 0.935
| $AUC$ | (0.030) | (0.026) | (0.027) | (0.034)
MNIST |  | 0.261 | 0.202 | 0.174 | 0.332
| $logloss$ | (0.046) | (0.046) | (0.045) | (0.027)

Table 8: ResNet18 vs. Its Potential Improvers

We have also investigated whether LIFE algorithm can further improve the
performance of trained deep neural network on image data. We here use ResNet18
proposed by (He et al. 2016 [13]) as an example and apply the algorithms on
MNIST data [19]. The detailed information of data preprocessing can be found
in Case 8 from the description lists of real data sets Classification. After
training MNIST data using ResNet18, the output of final convolutional layer
has been extracted, which has size $8000\times 512$, and it is also the input
of feed forward neural network (FFNN) in the final step of ResNet18. This
$8000\times 512$ data is then treated as the input of LIFE(Adam). Further, we
also attach XGBoost, Adam to ResNet18 and compare the results.

####  3.2.3 Classification on Image Data

Figure 9: ResNet18 vs. Its Potential Improvers

Figure 9 shows the performance of LIFE (orange line) is better than ResNet18
(blue line) with consistently larger AUC and smaller log-loss in all of the 10
replications. LIFE is comparable to XGboost in terms of log-loss in general
with all 10 values below ResNet18. However, in terms of AUC, XGboost is worse
than ResNet18 with significantly lower AUC in two replications (seed=2 and
seed=8). On the other hand, Adam (red line) is not able to further improve the
performance of ResNet18, which is consistent with what we have discovered in
the previous empirical studies. Table 8 also indicates LIFE and Xgboost are
both capable of remarkably enhancing a trained deep NN with similar
performance. One last discovery is since the input of LIFE contains 512
columns, it also indicates that LIFE can handle high dimensionality quite well
in terms of prediction.

##  4 Interpretation

Interpretability is the degree to which one human being can understand the
cause of a decision or predict the result of a model. The higher the
interpretability of a machine learning or deep learning model, the easier it
is for someone to comprehend why certain decisions or predictions have been
made. A key advantage of LIFE is that it is still an inherently interpretable
model. From the perspective of the NN structure, the model is a single-hidden-
layer NN with ReLU activation function where all the weights and bias can be
easily extracted and visualized. Moreover, the single layer NN with ReLU
activation function can be rewritten in the form of local linear model
representation, and be interpreted by exploring the patterns of local linear
model coefficients. Finally, the main and interaction effects can be
identified by exploring and aggregating the local linear coefficients.

We use the bike sharing data result as an example to illustrate the intrinsic
interpretability of LIFE. Bike sharing data is a public dataset hosted on UCI
machine learning repository, where there are around $17,000$ observations on
hourly (and daily) bike rental counts along with weather and time information
between $2011$ and $2012$ in the Capital Bikeshare system. Out of the original
$17$ predictors, we removed some non-meaningful and highly correlated ones,
leaving us with $9$ predictors to predict hourly rental counts. At the tiny
expense of predictive performance, we applied both the base learner selection
method shown in Algorithm 3 in Section 5.2 and elastic net to reduce the
number of base learners and features so that the final single hidden-layer NN
has a small number of signficant neurons and is easier for interpretation.

(a)

(b)

Figure 10: Variable Importance Detection

###  4.1 Explore the weights and bias of single layer NN

As LIFE finally generates a single-hidden-layer NN in the third step, we can
explore the weights $\hat{w}_{k}$s and bias $\hat{b}_{k}$s of the NN directly
and identify which variable is important. For bike sharing data, there are
finally $116$ new features (or neurons) after base learner selection and
elastic net regularization. We measure the neuron importance by
$std(\hat{\beta}_{k}\sigma(\hat{b}_{k}+x^{T}\hat{w}_{k}))/std(\hat{f})$, and
$\hat{f}$, where std is the standard deviation, $\hat{f}$ is the predicted
value of response variable for regression or log-odds for classification,
$\hat{w}_{k}$s is the neuron weight, and $\hat{b}_{k}$ is the coefficient for
neuron. This quantity measures the importance of neurons/feature by comparing
the variation of each feature to the total variance. The histogram on the
neuron importance for the $116$ features in Figure 10a shows that there are
only $13$ neurons whose importance values are greater than $2\%$ of the
maximum importance.

Then we can detect how each variable contributes to each neuron by applying
the following measurement:

$\hat{w}_{k}\hat{\beta}_{k}\frac{std(\sigma(\hat{b}_{k}+x^{T}\hat{w}_{k}))}{std(\hat{f})}$

where we allocate neuron importance to each variable by multiplying
$\hat{w}_{k}$. This contribution measurement can be simply visualized by
heatmap between neurons and original variables in Figure 10b . It shows that
hour(hr) and working day(workingday) are top significant variables with darker
colors for almost each important neuron compared with other variables. Another
variable temperature(temp) can also be considered to relatively important
except hr and workingday.

###  4.2 Treat a single layer NN as a local linear model

As we may have many features in the final wide single layer NN, it is
difficult to visualize and explore the weights and bias of all the neurons.
Hence, we also propose to interpret single layer NN from local linear model
perspective. Single-layer NN with ReLU function can be considered a type of
local linear model. Each linear projection would determine the active or
inactive states of the ReLU neurons at hidden layers, which define the layered
pattern. The activation region is constructed as a combination of those
distinct patterns. Those activation regions are mutually exclusive and
regarded as convex polytopes with closed-form boundaries [23]. A linear
equation can be used in all data points inside the activation region to
represent the relationships between response and independent variables. After
defining the region each observation belongs to, we can easily extract a
linear equation for each region based on estimated weights in the hidden and
output layers. The detailed algorithm that performs a linear equation
extraction the following:

Input: Estimated weights and biases
$\\{\hat{w}_{k},\hat{b}_{k}\hat{\beta}_{k}\\}_{k=1,\cdots,m_{J}}$; $\tau$:
threshold for the number of observations

Output: Coefficients of local linear equations
$\\{\hat{E}_{t}\\}_{t=1,\cdots,M}$

Data: Independent variables $\\{x_{i}\\}_{i=1,\cdots,N}$

1 Determine if $i^{th}$ observation is active or not in each neuron based on
linear projection $\hat{b}_{k}+x_{i}^{T}\hat{w}_{k}>0$

2 Combine observations into $M$ homogenous activation local regions
$R_{t},t=1,\cdots,M$, depending on if they have the same set of active neurons
( boundary condition)

3 for  _$t=1,\cdots,M$_ do

4 if  _ $\sum_{x_{i}\in R_{t}}i>\tau$ _ then

5 1\. Construct activation set $A_{t}$ of the local region $R_{t}$ relying on
its boundary condition

6 2\. Calculate coefficients: $\hat{E}_{t}=\sum_{k\in
R_{t}}\hat{w}_{k}*\hat{\beta}_{k}$

7

Algorithm 2 Local Linear Equation Extraction

We can visualize those linear equations by a parallel coordinate plot, which
allows comparing the estimated coefficients of all predictors for different
local linear regions. Through the visualization of local linear equations, we
can not only have an overview of the importance of each predictor in each
region by comparing the magnitude of coefficients, but also check the validity
for effect of each predictor on the response variables. It is worth mentioning
that those coefficients are comparable after standardizing all the predictors.
There are three scenario for a particular independent variable:

  1. 1.

Relatively large coefficients of the variable, compared with others in terms
of absolute values and have the same signs, imply that this variable has a
significant positive or negative effect on the response variable if all
coefficients are positive or negative.

  2. 2.

Relatively large coefficients of the variable with both positive and negative
signs strongly imply that this variable has inconsistent slopes across local
activation regions, which might be due to either its own nonlinear main effect
or the interaction effects with other variables.

  3. 3.

Small and close-to-zero coefficients indicate that this feature is not
important to explain the variation of the response variable and can be removed
from the model.

Furthermore, we were able verify if the sign of estimated coefficients of
predictor in all regions is consistent with domain knowledge or business
sense. Figure 11 displays the estimated coefficients of all predictors in the
local activation regions for bike sharing dataset. There are $116$ neurons
extracted from NN base learners and $47$ local regions created by Algorithm 2
and each local region has at least three data points.

Figure 11: Parallel Corrdiantes Plot for Bikesharing Data

It clearly indicates that hour(hr), working day(workingday) and
temperature(temp) are the three most important predictors with relatively
higher absolute values of their corresponding coefficients in several local
regions, which is pretty consistent with result from Figure 10b. Their
estimated coefficients present different directions across local activation
regions, which is consistent with our second scenario.

This gives us a hint of interactions between those variables Other variables
such as humidity(hum) and wind speed(windspeed) are insignificant based on
their absolute values of estimation coefficients from the plot. Sometimes
there are too many local regions and (Sudjianto et al.(2020))[23] provides two
approaches to simplify and reduce the number of local linear equations—merging
and flattening in their paper, where a variety of other diagnostic tools and
plots for local linear model have also been provided.

###  4.3 Main and Interaction Effect Detection

Even though the parallel coordinates plot provides a guideline about the
variable importance in each local region, we still need a solid technique to
detect nonlinear main effects and interaction effects. To achieve this
purpose, we can treat single-hidden-layer NN as a varying coefficient model
through linear equation extraction shown in Algorithm 2. As all the local
linear equation coefficients are varying over local regions, and region
definition depends on predictors, so the coefficients can be treated as a
function of predictors in equation 11.

$\begin{gathered}\hat{f}_{i}=\alpha_{0i}+\alpha_{1i}x_{1i}+\cdots+\alpha_{mi}x_{mi}+\cdots+\alpha_{pi}x_{pi},\
\ \ i=1,\cdots,n,\\\ \text{where
}\alpha_{mi}=g(x_{1i},\cdots,x_{pi}),\end{gathered}$ (11)

where $p$ is the number of predictors and $n$ is number of observations.
$\hat{f}_{i}$ is predicted value for regression and predicted log-odds for
classification. $\alpha_{mi}$ is the coefficient for $m^{th}$ variable at
$i^{th}$ observation, and could also be a function of all predictors, varying
by different observations. Our goal is to investigate what the functional
forms of the estimated coefficients are. Therefore, we separate $\alpha_{mi}$
into two components representing main and interaction effects in equation 12:

$\begin{gathered}\alpha_{mi}=g(x_{1i},\cdots,x_{pi})=\underbrace{g_{main}(x_{mi})}_{\text{main
effect}}+\underbrace{g_{int}(x_{1i},\cdots,x_{pi})}_{\text{interaction
effect}}.\\\ \end{gathered}$ (12)

The first term in the equation 12 is a function of $x_{mi}$ , including the
intercept of $\alpha_{mi}$, and this term captures the main effect of
$x_{mi}$. If $\alpha_{mi}$ has a significant intercept, then linear main
effect can be identified; while a strong relationship with $x_{mi}$ indicates
a nonlinear main effect. The remaining second term is the function of other
predictors and it may or may not contain $x_{mi}$. This term can be used to
detect interactions between $x_{mi}$ and other predictors. For an
illustration, let us look at a simple example with all estimated coefficients
constant except $\alpha_{1i}=\theta_{0}+\theta_{1}x_{1i}+\theta_{2}x_{2i}$,
then the varying coefficient model can be expressed as follows:

$\begin{gathered}y_{i}=\alpha_{0i}+\theta_{0}x_{1i}+\theta_{1}x_{1i}^{2}+\theta_{2}(x_{2i})x_{1i}+\alpha_{2}x_{2i}+\cdots+\alpha_{p}x_{pi}+\epsilon_{i}.\end{gathered}$
(13)

where $g_{main}(x_{1i})=\theta_{0}+\theta_{1}x_{1i}$ and
$g_{int}(x_{2i})=\theta_{2}x_{2i}$. We can easily identify the interaction
term between $x_{1i}$ and $x_{2i}$, and $x_{1i}$ shows a nonlinear main effect
via its quadratic term. To detect main effects and interaction effects from
$\alpha_{mi}$, we propose the two-stage process below:

  1. 1.

Check nonlinearity: Calculate conditional expectation
$E(\hat{\alpha}_{mi}|x_{mi})$ by smoothing estimated coefficients of
predictors against itself 222The smoothing spline is used for illustration and
other estimation methods can also be applied..
$\hat{\alpha}_{mi}\sim g_{main}(x_{mi})\ \ \ m=1,\cdots,p$.

  2. 2.

Check interactions: Remove main effect from $\alpha_{mi}$, and calculate
conditional expectation $E(\hat{\alpha}_{mi}-\hat{g}_{main}(x_{mi})|x_{ki})$
by smoothing estimated coefficients of predictors against each other variable.
$\hat{\alpha}_{mi}-\hat{g}_{main}(x_{mi})\sim g_{k}^{m}(x_{ki})\ \ \ \ k\neq
m$

Figure 12: Plot Matrix Between $\alpha_{m}$ and $x_{m}$

We choose to use a two-stage process instead of a one-stage process, as we can
estimate its main and interaction effect more accurately in the correlated
predictor case and split two effects effectively. Note that some special
interaction effects may not be identified by one-stage process such as
$y=\alpha_{0}+\alpha_{1}x_{1}$ as an example, where $\alpha_{1}=x_{1}x_{2}$.
In this case, $g_{2}^{1}(x_{2i})$ is zero curve. Fortunately, most common
interaction patterns can be identified by our two stage process. As long as
$g_{k}^{m}(x_{ki})$ has a significant pattern on $x_{ki}$, an interaction
effect can be identified.

For the case of bike sharing data, we visualized all pairs of varying
coefficients and variables $(\alpha_{mi}\ vs\ x_{ki})$ with scatter plots in
Figure 12. Due to $\frac{\partial\hat{f}(\bf{x})}{\partial x_{m}}=\alpha_{m}$
, this is also scattered partial derivative plot for $\hat{f}(\bf{x})$ . On
top of the scatter plot, we also draw $g_{main}(x_{mi})$ against $x_{mi}$ in
the diagonal plots show and $g_{k}^{m}(x_{ki})$ agaist $x_{ki}$ in the $(m,k)$
off-diagonal plots . To further quantify the magnitude of the interaction
effects, we calculated weighted standard deviation of $g_{main}(x_{mi})$ and
$g_{k}^{m}(x_{ki})$ with population density as weight. The heatmap of the
interaction measures for bike sharing data is provided in Figure 13 where
diagonals are masked by zero.

Figure 13: Heatmap for Interaction Measures Figure 14: ALE Plots for
Predictors Based on LIFE

The nonlinear patterns of variables can be clearly spotted in the diagonal
plots in Figure 12. The most important variable hour(hr) displays the drastic
fluctuation compared with others, indicating its nonlinear effect on response.
As evidenced in both Figure 12 and Figure 13, the top three interaction pairs
including hr vs workingday, hr vs weekday and hr vs temp can be easily
identified.

In addition to interaction detection, we can obtain and visualize the main
effect of each predictor directly by aggregating the local linear
coefficients. Due to $\frac{\partial\hat{f}(\bf{x})}{\partial
x_{m}}=\alpha_{m}$ in the varying coefficient setting, we compute the exact
main effect of $x_{m}$ by constructing a relationship between $f(x_{j})$ and
$x_{j}$ based on formula
$\int_{0}^{x}E(\frac{\partial\hat{f}(\bf{x})}{\partial x_{m}}|x_{m})dx_{m}$
from Accumulated Local Effects (ALE) plot, where the variable is transformed
back to original scale as seen in Figure 14. This ALE formulation can be
simplified as
$\int_{0}^{x}E(\alpha_{m}|x_{m})dx_{m}=\int_{0}^{x}g_{main}(x_{m})dx_{m}$ and
its numerical implementation of ALE is achieved by the Midpoint Rule. The main
effect forhr has two peaks and one trough, which is similar to partial
dependence plot from other machine learning algorithms, while the main effect
of temp and hum show a quadratic relationship. More specifically, the peak of
bike rentals happen around 7 am and 5-6 pm, while very few people will rent
bikes around 3-4 am. People usually prefer to rent bikes in a nice day with
moderate temperature and humidity. Both of them are pretty consistent with
common sense.

##  5 Discussion

###  5.1 Different Sampling Schemes

The LIFE algorithm can be considered a general framework with three steps, as
discussed in the methodology section. LIFE is very flexible and allows users
to try different combinations of three steps. The first step presents several
data sampling options. In the paper, we use linear projection inside NN
neurons to split data and select data points from active region for base
learner training, as shown in Figure 2. Those linear projections are obtained
with trained NNs in a supervised setting. Given a fixed hyper-parameter setup
for LIFE, we have implemented our method with different sampling choices
including NN projection, Random projection, Bootstrapping. In Table 9, we can
easily see that all the ensemble methods outperform a single single-hidden-
layer NN model optimized by LLA or Adam. Most importantly, sampling by NN
linear projection is better than other sampling methods that create subsets in
a random way.

Model | Metric | NN projection | Random projection | Bootstrapping | LLA or Adam
---|---|---|---|---|---
|  | 1.060 | 1.095 | 1.138 | 1.110
| $RMSE$ | (0.017) | (0.023) | (0.034) | (0.065)
MIM (Regression) |  | 0.965 | 0.962 | 0.960 | 0.961
| $R^{2}$ | (0.001) | (0.002) | (0.003) | (0.005)
|  | 0.500 | 0.525 | 0.519 | 0.527
| $RMSE$ | (0.012) | (0.020) | (0.018) | (0.011)
California Housing |  | 0.812 | 0.792 | 0.797 | 0.791
| $R^{2}$ | (0.009) | (0.016) | (0.014) | (0.008)
|  | 0.940 | 0.935 | 0.936 | 0.920
| $AUC$ | (0.004) | (0.004) | (0.003) | (0.003)
MIM (Classification) |  | 0.268 | 0.277 | 0.275 | 0.421
| $logloss$ | (0.012) | (0.009) | (0.007) | (0.075)
|  | 0.938 | 0.929 | 0.927 | 0.909
| $AUC$ | (0.002) | (0.003) | (0.003) | (0.008)
Gamma Telescope |  | 0.288 | 0.307 | 0.314 | 0.371
| $logloss$ | (0.007) | (0.007) | (0.007) | (0.019)

  * •

Note: The different sampling methods in the first step of the general
framework of LIFE are tested on two simulated data and two real data. NN
projection performs data partition by linear projection inside neurons of NNs
and then samples data from the selected region, which is the main part of LIFE
algorithm used in the paper. Random projection performs data partition by
random linear projection, where weights and biases are randomly drawn from
standard normal. The bootstrapping selects observations randomly from training
set. The number colored in red represents optimal result for this metric. LLA
is used to optimize single-hidden-layer NN base learner in the regression
case, which Adam is used for classification case.

Table 9: Sampling Method Comparison

###  5.2 Base Learner Selection

For model aggregation and pruning, we can prune neurons to have a single layer
NN with fewer neurons. We used Elastic net for pruning due to its simplicity,
but pruning methods besides Elastic net can also be considered. Based on
properties of LIFE, we also developed an alternative pruning method called
base learner selection to reduce the number of nodes in the final step. It is
assumed that LIFE works well because the correlations of prediction errors are
not strong among different base learners. Therefore, we can remove base
learner one by one, according to the correlation between its prediction errors
and prediction errors from other base learners. In this way, we can still
maintain diversity and solve overfitting issue by keeping fewer necessary base
learners without sacrificing predictive performance a lot. This method is
thoroughly described by Algorithm 3.

Input: $B_{j}$: trained base learner $j=1,\cdots,M$; $\tau$: threshold

Output: $B_{j}$: selected base learner $j\in\\{1,\cdots,M\\}$

Data: Training set $\\{x_{i},y_{i}\\}_{i=1,\cdots,N}$

1 for  _$j=1,\cdots,M$_ do

2 1\. $\hat{y}_{i}^{(j)}=B_{j}(x_{i})$

3 2\. $r_{i}^{(j)}=y_{i}-\hat{y}_{i}^{(j)}$

4for  _$j=1,\cdots,M$_ do

5 1\. Train a linear model $r_{i}^{(j)}\sim
c_{0}+c_{1}r_{i}^{(1)}+\cdots+c_{j-1}r_{i}^{(j-1)}+c_{j+1}r_{i}^{(j+1)}+\cdots+c_{M}r_{i}^{(M)}$

6 2\. Collect $R^{2}_{j}$ from this linear model

Remove $B_{j}$ based on value of $R^{2}_{j}$ in descending order until $\tau$
is achieved (the one with highest $R^{2}$ will be removed first).

Algorithm 3 Base learner selection

In Algorithm 3, the threshold $\theta$ is the percentage of base learners you
want to retrieved from a pool of candidates. When there is a large number of
neurons in LIFE setting, the elastic-net is usually computational expensive
and the base learner selection is a good alternative by parallel computation.
Moreover, we can combine these two pruning methods to achieve a simpler NN
model from wide NN faster. We also illustrate it using two simulated model
(GAM and MIM). The plots (a) and (b) in Figure 15 show the relationship
between $R^{2}$ and number of hidden neurons for the feature extraction.
Setting different thresholds in the base learner selection Algorithm 3, we can
construct single-hidden-layer NN with different number of neurons. In general,
base learner selection algorithms can effectively reduce number of neurons to
produce a simpler model without sacrificing predictive performance or obtain
even better results.

(a)

(b)

Figure 15: Relationship between $R^{2}$ and number of neurons. The red point
indicates one using linear regression as a final step without base learner
selection, while blue points indicates model aggregation by base learner
selection with different threshold.

##  6 Conclusion

In this paper, we have proposed a novel algorithm that fits single-hidden-
layer NN to achieve three goals: ensuring competitive predictive performance,
boosting computational efficiency, and preserving the interpretability of the
model. Unlike traditional NN training methods, we train it in an iterative way
through multiple NNs layer-by-layer training and then effectively combine them
via neural nodes flattening. We have evaluated the performance of our approach
using simulated and empirical data in terms of predictive accuracy and
computational efficiency and found that it consistently outperforms single-
hidden-layer NN trained directly by LLA or Adam optimizer and achieves
competitive results as those of Xgboost.

This superior performance lies in three reasons: First, as an ensemble method,
the LIFE algorithm performs data sampling through linear projection inside
neural nodes, which creates diversity among the models and contributes to bias
and variance reduction of prediction from combined models. Second, the LIFE
algorithm takes advantage of single-hidden-layer NN structure to combine
multiple narrow single-hidden-layer NNs into a wide one via neural nodes
flattening Third, LIFE algorithm benefits from leveraging parallel computing
to train multiple NNs on subsets of data simultaneously. Moreover, the base
learner selection method is introduced in the paper to help us prune redundant
neural nodes and produce a more parsimonious model after several iterations of
the LIFE algorithm. We have also proposed a new method for main and
interaction detection from the perspective of interpretation.

## References

  * Aksela and Laaksonen [2006] M. Aksela and J. Laaksonen.  Using diversity of errors for selecting members of a committee classifier.  _Pattern Recognition_ , 39(4):608–623, 2006.
  * Ancona et al. [2017] M. Ancona, E. Ceolini, C. Öztireli, and M. Gross.  Towards better understanding of gradient-based attribution methods for deep neural networks.  _arXiv preprint arXiv:1711.06104_ , 2017.
  * Andrienko and Andrienko [2001] G. Andrienko and N. Andrienko.  Constructing parallel coordinates plot for problem solving.  In _1st International Symposium on Smart Graphics_ , pages 9–14, 2001\.
  * Breiman [2001] L. Breiman.  Random forests.  _Machine learning_ , 45(1):5–32, 2001.
  * Brown et al. [2005a] G. Brown, J. Wyatt, R. Harris, and X. Yao.  Diversity creation methods: a survey and categorisation.  _Information Fusion_ , 6(1):5–20, 2005a.
  * Brown et al. [2005b] G. Brown, J. L. Wyatt, and P. Ti o.  Managing diversity in regression ensembles.  _Journal of machine learning research_ , 6(Sep):1621–1650, 2005b.
  * Butler IV et al. [2018] H. K. Butler IV, M. A. Friend, K. W. Bauer Jr, and T. J. Bihl.  The effectiveness of using diversity to select multiple classifier systems with varying classification thresholds.  _Journal of Algorithms & Computational Technology_, 12(3):187–199, 2018.
  * Chen et al. [2020] J. Chen, J. Vaughan, V. Nair, and A. Sudjianto.  Adaptive explainable neural networks (axnns).  _Available at SSRN 3569318_ , 2020.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin.  Xgboost: A scalable tree boosting system.  In _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , pages 785–794, 2016.
  * Chen et al. [2015] T. Chen, T. He, M. Benesty, V. Khotilovich, and Y. Tang.  Xgboost: extreme gradient boosting.  _R package version 0.4-2_ , pages 1–4, 2015.
  * Gacquer et al. [2009] D. Gacquer, V. Delcroix, F. Delmotte, and S. Piechowiak.  On the effectiveness of diversity when training multiple classifier systems.  In _European conference on symbolic and quantitative approaches to reasoning and uncertainty_ , pages 493–504. Springer, 2009.
  * Hansen [2000] J. V. Hansen.  _Combining predictors: Meta machine learning methods and bias/variance & ambiguity decompositions_.  PhD thesis, Aarhus University, Computer Science Department, 2000.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Heath et al. [1993] D. Heath, S. Kasif, and S. Salzberg.  Induction of oblique decision trees.  In _IJCAI_ , volume 1993, pages 1002–1007, 1993.
  * Krogh and Vedelsby [1994] A. Krogh and J. Vedelsby.  Neural network ensembles, cross validation, and active learning.  _Advances in neural information processing systems_ , 7:231–238, 1994.
  * Kucherenko et al. [2009] S. Kucherenko et al.  Derivative based global sensitivity measures and their link with global sensitivity indices.  _Mathematics and Computers in Simulation_ , 79(10):3009–3017, 2009.
  * Kucherenko et al. [2010] S. Kucherenko et al.  A new derivative based importance criterion for groups of variables and its link with the global sensitivity indices.  _Computer Physics Communications_ , 181(7):1212–1217, 2010.
  * Kuncheva and Whitaker [2003] L. I. Kuncheva and C. J. Whitaker.  Measures of diversity in classifier ensembles and their relationship with the ensemble accuracy.  _Machine learning_ , 51(2):181–207, 2003.
  * Lecun et al. [1999] Y. Lecun, C. Cortes, and C. Burges.  The mnist dataset of handwritten digits(images), 1999.
  * Qiu et al. [2018] J. Qiu, S. R. Jammalamadaka, and N. Ning.  Multivariate bayesian structural time series model.  _Journal of Machine Learning Research_ , 19:1–33, 2018\.
  * Qiu et al. [2020] J. Qiu, S. R. Jammalamadaka, and N. Ning.  Multivariate time series analysis from a bayesian machine learning perspective.  _Annals of Mathematics and Artificial Intelligence_ , 88(10):1061–1082, 2020.
  * Rudin [2019] C. Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ , 1(5):206–215, 2019.
  * Sudjianto et al. [2020] A. Sudjianto, W. Knauth, R. Singh, Z. Yang, and A. Zhang.  Unwrapping the black box of deep relu networks: Interpretability, diagnostics, and simplification.  _arXiv preprint arXiv:2011.04041_ , 2020.
  * Sundararajan et al. [2017] M. Sundararajan, A. Taly, and Q. Yan.  Axiomatic attribution for deep networks.  _arXiv preprint arXiv:1703.01365_ , 2017.
  * Ueda and Nakano [1996] N. Ueda and R. Nakano.  Generalization error of ensemble estimators.  In _Proceedings of International Conference on Neural Networks (ICNN’96)_ , volume 1, pages 90–95. IEEE, 1996.
  * Vaughan et al. [2018] J. Vaughan, A. Sudjianto, E. Brahimi, J. Chen, and V. N. Nair.  Explainable neural networks based on additive index models.  _arXiv preprint arXiv:1806.01933_ , 2018.
  * Wolpert [1992] D. H. Wolpert.  Stacked generalization.  _Neural networks_ , 5(2):241–259, 1992.
  * Yang et al. [2020] Z. Yang, A. Zhang, and A. Sudjianto.  Enhancing explainability of neural networks through architecture constraints.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Zeng et al. [2020] M. Zeng, Y. Liao, R. Li, and A. Sudjianto.  Local linear approximation algorithm for neural network.  _Manuscript_ , 2020.

## Appendix

### Loss Decomposition in the Log-odds Space

The cross-entropy loss for one observation in the log-odds space can be
written as:

$\begin{gathered}l(y_{i},f_{i})=y_{i}log(1+e^{-f_{i}})+(1-y_{i})log(1+e^{f_{i}})\end{gathered}$
(14)

Therefore, the average cross entropy loss is expressed as:

$\begin{gathered}\sum_{i=1}^{N}[y_{i}log(1+e^{-f_{i}})+(1-y_{i})log(1+e^{f_{i}})]\\\
=\underbrace{\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}[y_{i}log(1+e^{-f_{i}^{(j)}})+(1-y_{i})log(1+e^{f_{i}^{(j)}})]}_{\text{accuracy}}\\\
-\underbrace{\sum_{i=1}^{N}\sum_{j=1}^{M}\beta_{j}\\{\frac{e^{f_{i}^{(j)\star}}+e^{-f_{i}^{(j)\star}}+2}{[(1+e^{-f_{i}^{(j)\star}})(1+e^{f_{i}^{(j)\star}})]^{2}}\\}(f_{i}-f_{i}^{(j)})^{2}}_{\text{diversity}}\end{gathered}$
(15)

### Local Linear Approximation (LLA) Algorithm

The algorithm is extracted from Zeng, et al (2020) [29]. Note that
$h^{(0)}=x$. Then one hidden-layer FNN with $J_{1}$ nodes can be expressed as
follows
$h^{(1)}=m{\sigma}(Wh^{(0)}+b)=(\sigma(w^{T}_{1}x+b_{1}),\cdots,\sigma(w^{T}_{J_{1}}x+b_{J_{1}}))^{T},$
where we drop the subscripts of $W$ and $b$ for ease of presentation. We set
$\sigma(z)=\max\\{z,0\\}$, the ReLU activation function. To estimate the
weights and biases, we minimize a nonlinear LS function

$l_{1}(\theta)=\sum_{i=1}^{n}\\{y_{i}-\beta_{0}-\sum_{j=1}^{J_{1}}\beta_{j}\sigma(w_{j}^{T}x_{i}+b_{j})\\}^{2}$

Given $w_{j}^{(c)}$ and $b_{j}^{(c)}$ in the current step, we propose to
approximate $\sigma(x^{T}w_{j}+b_{j})$ by a linear function based on the
first-order Taylor expansion of $\sigma(z)$:

$\sigma(x^{T}w_{j}+b_{j})\approx\sigma(x^{T}w_{j}^{(c)}+b_{j}^{(c)})+\\{(x^{T}w_{j}+b_{j})-(x^{T}w_{j}^{(c)}+b_{j}^{(c)})\\}I(x^{T}w_{j}^{(c)}+b_{j}^{(c)}>0)$

for $j=1,\cdots,J_{1}$. Thus,

$\beta_{j}\sigma(x^{T}w_{j}+b_{j})\approx\beta_{j}\sigma(x^{T}w_{j}^{(c)}+b_{j}^{(c)})+\gamma_{j}I(x^{T}w_{j}^{(c)}+b_{j}^{(c)}>0)+m{\eta}_{j}^{T}xI(x^{T}w_{j}^{(c)}+b_{j}^{(c)}>0)$

where $\gamma_{j}=\beta_{j}(b_{j}-b_{j}^{(c)})$ and
$m{\eta}_{j}=\beta_{j}(w_{j}-w_{j}^{(c)})$. Define
$z_{1ij}=\sigma(x_{i}^{T}w_{j}^{(c)}+b_{j}^{(c)})$,
$z_{2ij}=I(x_{i}^{T}w_{j}^{(c)}+b_{j}^{(c)}>0),$ and
$z_{3ij}=x_{i}I(x_{i}^{T}w_{j}^{(c)}+b_{j}^{(c)}>0)\\},$ $j=1\cdots,J_{1}$.
Further define $z_{1i}=[z_{1i1},\cdots,z_{1iJ_{1}}]$,
$z_{2i}=[z_{2i1},\cdots,z_{2iJ_{1}}],$
$z_{3i}=[z_{3i1}^{T},\cdots,z_{3iJ_{1}}^{T}],$ and
$z_{i}=[z_{1i},z_{2i},z_{3i}]^{T}$, which is a $J_{1}(p+2)$-dimensional
vector. the objective function is approximated by

$\sum_{i=1}^{n}\\{y_{i}-\beta_{0}-\sum_{j=1}^{J_{1}}\\{\beta_{j}z_{1ij}+\gamma_{j}z_{2ij}+m{\eta}_{j}^{T}z_{3ij}\\}\\}^{2},$

which is the LS function of linear regression with the response $y_{i}$ and
predictors $z_{i}$. Denote the resulting LS estimate of $\beta_{j}$,
$\gamma_{j}$ and $m{\eta}_{j}$ by $\hat{\beta}_{j}$, $\hat{\gamma}_{j}$ and
$\hat{m{\eta}}_{j}$, respectively. By the definition of $\gamma_{j}$ and
$m{\eta}_{j}$, we can update $b_{j}$ and $w_{j}$ as shown in the step 2 of
algorithm If $|\hat{\beta}_{j}|$ is very close to zero, one may simply set
$b_{j}^{(c+1)}=b_{j}^{(c)}$ and $w_{j}^{(c+1)}=w_{j}^{(c)}$. Thus, we may
estimate $W$ and $b$ by iteratively and regressing $y_{i}$ on the updated
$z_{i}$. The procedure can be summarized as the following algorithm.

  1. 1.

Set initial value for $W^{(0)}=[w_{1}^{(0)},\cdots,w_{J_{1}}^{(0)}]^{T}$ and
$b_{j}^{(0)}$, and let $c=0$.

  2. 2.

Calculate $z_{i}$ defined in the text based on $w_{j}^{(c)}$ and
$b_{j}^{(c)}$, obtain the least squares estimate (LSE) $\hat{\beta}_{j}$s,
$\hat{\gamma}_{j}$s and $\hat{\eta}_{j}$s by running a linear regression
$y_{i}$ on covariate $z_{i}$, and update the biases and weights by

$b_{j}^{(c+1)}=b_{j}^{(c)}+\hat{\gamma}_{j}/\hat{\beta}_{j},\ \text{and}\
w_{j}^{(c+1)}=w_{j}^{(c)}+\hat{\eta}_{j}/\hat{\beta}_{j}.$

if $|\hat{\beta}_{j}|\geq\varepsilon$, where $\varepsilon$ is a constant for
numerical stability and is set to be $10^{-3}$ in our numerical experiment,
and keep the corresponding biases and weights unchanged if
$|\hat{\beta}_{j}|<\varepsilon$.

  3. 3.

Set $c=1,2,\cdots,$ and repeat Step 2 until the criterion of algorithm
convergence meets.

### Functional Forms for Classification Case in Simulation Study

In classification case, we use three function forms including Generalized
Additive Model (GAM), Additive Index Model (AIM) and Multiple Index Model
(MIM) expressed as follows, where $\epsilon_{i}\sim N(0,1),i=1,\dots,N$.

$\begin{gathered}GAM:f(\boldsymbol{x}_{i})=1.5x_{1i}+4\sqrt{|-2.5x_{2i}|}+2|x_{3i}|+4\exp(-\frac{3}{14}x_{4i})+4\log(1.5|x_{5i}|)\\\
-4\max(1,x_{6i})+\epsilon_{i}.\end{gathered}$ (16)
$\begin{gathered}AIM:f(\boldsymbol{x}_{i})=\log(|3x_{1i}-2.5x_{2i}+2x_{3i}-1.5x_{4i}|)+\exp\\{(-1.5x_{4i}+1.5x_{5i}-x_{6i})/11\\}+\epsilon_{i}.\end{gathered}$
(17)
$\begin{gathered}MIM:f(\boldsymbol{x}_{i})=\exp(0.03x_{1i}-0.025x_{2i})x_{3i}-\frac{3x_{4i}}{1+1.5|x_{5i}|}+2\max(1,x_{6i})+\epsilon_{i}.\end{gathered}$
(18)

### Real Datasets

#### Regression

Case 1: Abalone. Original data comes from a Marine Resources Division Marine
Research Laboratories in Australia. The goal is to predict the age of abalone
from physical measurements, which is determined by cutting the shell through
the cone, staining it, and counting the number of rings through a microscope.
From the original data examples with missing values were removed, and the
ranges of the continuous values have been scaled for use. Therefore, there are
$4177$ observations left and $9$ variables including an integer response
variable (number of rings), one categorical predictor (sex) and seven
continuous predictors (Length, Diameter, Height, Whole weight, Shucked weight,
Viscera weight, Shell weight).

Case 2: Airfoil. It is NASA data set, obtained from a series of aerodynamic
and acoustic tests of two and three-dimensional airfoil blade sections
conducted in an anechoic wind tunnel. It comprises different size NACA 0012
airfoils at various wind tunnel speeds and angles of attack. There are $1503$
observations and $6$ variables, where the response variable is scaled sound
pressure level in decibels and there are five independent variables containing
frequency in Hertzs, angle of attack in degrees, chord length in meters, free-
stream velocity, in meters per second, suction side displacement thickness in
meters.

Case 3: Aquatic Toxicity. This dataset was used to develop quantitative
regression QSAR models to predict acute aquatic toxicity towards the fish
Pimephales promelas (fathead minnow) on a set of 908 chemicals. to predict
acute aquatic toxicity towards Daphnia Magna as a response variable called
LC50. Without missing values, it contains values for 8 predictors (molecular
descriptors) of 546 chemicals (TPSA, SAacc, H-050, MLOGP, RDCHI, GATS1p, nN,
C-040).

Case 4: Bike Sharing. This dataset contains the hourly and daily count of
rental bikes between years 2011 and 2012 in Capital bike share system with the
corresponding weather and seasonal information. Specifically, the bike sharing
dataset contains $17379$ data points of $16$ predictors. Out of the $16$
independent variables, we removed two non-meaningful information as well as
two response-related information. This leaves us with $12$ variables including
hour, temperature, feeling temperature, humidity, wind speed, season,
workingday weekday and weather situation, year and month.

Case 5: California Housing. This dataset was derived from the 1990 U.S.
census, using one row per census block group and the target variable is the
median house value for California districts. A block group is the smallest
geographical unit for which the U.S. Census Bureau publishes sample data (a
block group typically has a population of 600 to 3000 people). It consists of
20640 observations and the independent variables include longitude, latitude,
housing median age, medium income, population, total rooms, total bedrooms and
households.

Case 6: CASP. This is a dataset of Physicochemical Properties of Protein
Tertiary Structure, which is taken from CASP 5-9. The goal is to predict RMSD-
size of residue given other physical attributes. There are 45730 observations
and $9$ predictors including total surface area, non-polar exposed area,
fractional area of exposed non polar residue, fractional area of exposed non
polar part of residue, molecular mass weighted exposed area, average deviation
from standard exposed area of residue, Euclidian distance, secondary structure
penalty, special Distribution constraints (N,K Value).

Case 7: Electrical Grid. The local stability analysis of the 4-node star
system (electricity producer is in the center) implementing Decentral Smart
Grid Control concept. There are $10000$ observations and 11 predictors
including $tau[x],x=1,2,3,4$ which are reaction time of participant,
$p[x],x=2,3,4$ which are nominal power consumed (negative) divided by
produced(positive)(real), and $g[x],x=1,2,3,4$ that are coefficient (gamma)
proportional to price elasticity. The continuous response variable is the
maximal real part of the characteristic equation root.

#### Classification

Case 1: Bank Marketing. The data provides information regarding direct
marketing campaigns of a Portuguese banking institution. The classification
goal is to predict if the clients, who were contacted based on at least two
phone calls in general, would like to subscribe a term deposit or not. The
data contains 45211 examples and 16 variables including 6 numerical variables
(age, balance, day, campaign, pdays, previous) and 9 categorical variables
(job, marital, education, default, housing, loan, contact, month, poutcome).
Notice that input variable ‘duration’ is not included in the analyses based on
the suggestions from the provider of this dataset since it highly affects the
output target and thus the variable should be discarded due to the intention
of building a realistic predictive model.

Case 2: Breast Cancer Wisconsin. Data comes from original Wisconsin Breast
Cancer Database, with purpose of detecting if the tissue is benign or
malignant. The original dataset contains 699 instances and 9 attributes. After
deleting rows with missing values, finally 683 instances are used in our
analyses. All predictive variables are numerical with a scale of 1 to 10,
including Clump Thickness, Uniformity of Cell Size, Uniformity of Cell Shape,
Marginal Adhesion, Single Epithelial Cell Size, Bare Nuclei, Bland Chromatin,
Normal Nucleoli and Mitoses.

Case 3: Higgs Boson. This dataset was built from official ATLAS full-detector
simulation in 2014 that mixes ‘Higgs to tautau’ events with different
backgrounds, where the events in which Higgs bosons were produced are
comprised in the signal sample, while other known processes mimicking the
signal are considered as background noise. The dataset contains 818238 events
and 63 variables, among which a few are highly correlated. The objective is to
detect signal from background based on characteristics of events such as mass
(estimated, transverse or invariant), transverse momentum, centrality of the
pseudo rapidity of the lepton, azimuth angle, etc.

Case 4: Home Lending. This is an internal First Lien modeling data from Wells
Fargo, which was constructed from the entire panel dataset of loans provided
from the Consumer Lending Group’s (CLG) data modeling team. The goal is to
predict the probability of troubled loans based on fico score, unemployment
rate, delinquency status, loan to value ratio (LTV), total employment, etc.
The toy dataset used in this paper contains 1 million observations and 43
variables, which was selected from 210 million observations obtained by
stacking. Training and testing data are split based on ‘account id’ to prevent
overfitting.

Case 5: MAGIC Gamma Telescope. The dataset comes from Major Atmospheric Gamma
Imaging Cherenkov Telescope project (MAGIC), which is Monte-Carlo generated to
simulate registration of high energy gamma particles in a ground-based
atmospheric Cherenkov gamma telescope using the imaging technique. The
objective is to discriminate the images caused by signals (primary gammas)
from the images caused by background (hadronic showers initiated by cosmic
rays) based on some characteristics such as the attributes of ellipse
(fLength, fWidth, fM3Long, fM3Trans, fAlpha, fDist) and attributes on pixels
(fSize, fConc, fConc1, fAsym). The dataset contains 19020 observations and 10
continuous predictive variables.

Case 6: Mushroom. This dataset includes descriptions in terms of physical
characteristics of samples related with 23 species of gilled mushrooms. Each
species is identified as poisonous or edible. The purpose is to predict the
probability of a species being poisonous given the attributes of the species
such as cap (shape, surface, color), odor, gill (attachment, spacing, size,
color), stalk (surface and color above or below ring), veil (type, color),
ring (number, type), spore print color, population and habitat, all of which
are categorical variables. The dataset contains 8124 instances and 21
variables.

Case 7: Ringnorm. It is a simulated dataset as an implementation of Leo
Breimans ringnorm example with 7400 rows and 21 columns. Two classes contained
in the response variable, within each class explanatory variables are drawn
from a multivariate normal distribution.

Case 8: MNIST Data. Data are originally obtained from the MNIST database of
handwritten digits collected from hundreds of Census Bureau employees and
high-school students. The response variable contains signal digit from $0$ to
$9$, which are further centered into $28\times 28$ gray-scaled images. Then
the images are translated into regular dataset with dimension 60,000 by 784,
indicating there are 60,000 images and $28\times 28$ positions of the field.
Digit ‘1’ is considered as signal and all the other 9 digits are considered as
noise and are re-encoded as ‘0’. Training, validation and testing data are
randomly selected from these 60,000 observations for 10 times with size 8000,
2000, 2000, respectively. Images in the training dataset are then further
rotated by 20 degrees.