Individualized and Global Feature Attributions for Gradient Boosted Trees in
the Presence of $\ell_{2}$ Regularization

Qingyao Sun

University of Chicago

###### Abstract

While $\ell_{2}$ regularization is widely used in training gradient boosted
trees, popular individualized feature attribution methods for trees such as
Saabas and TreeSHAP overlook the training procedure. We propose Prediction
Decomposition Attribution (PreDecomp), a novel individualized feature
attribution for gradient boosted trees when they are trained with $\ell_{2}$
regularization. Theoretical analysis shows that the inner product between
PreDecomp and labels on in-sample data is essentially the total gain of a
tree, and that it can faithfully recover additive models in the population
case when features are independent. Inspired by the connection between
PreDecomp and total gain, we also propose TreeInner, a family of debiased
global feature attributions defined in terms of the inner product between any
individualized feature attribution and labels on out-sample data for each
tree. Numerical experiments on a simulated dataset and a genomic ChIP dataset
show that TreeInner has state-of-the-art feature selection performance. Code
reproducing experiments is available at https://github.com/nalzok/TreeInner.

##  1 INTRODUCTION

The interpretability of machine learning algorithms has been receiving
research interest due to its relevance to real-life problems in domains
including genomics (Basu et al.,, 2018; Kumbier et al.,, 2018), healthcare
(Rudin and Ustun,, 2018; Caruana et al.,, 2015), safety (Kim,, 2019), data
privacy (Alekh,, 2018), and criminal justice (Rudin et al.,, 2020). Tree
ensembles such as random forests (Breiman,, 2001) and gradient boosted trees
(Friedman,, 2001) are of particular interest because they tend to have high
predictive accuracy for tabular datasets and the tree structure is more
interpretable compared to neural networks (Lundberg et al.,, 2020). In
addition, tree ensembles are widely available in many open-source ML packages
(Chen and Guestrin,, 2016; Ke et al.,, 2017; Prokhorenkova et al.,, 2018).

In interpretable machine learning (Doshi-Velez and Kim,, 2017), individualized
feature attribution (IFA) and global feature attribution (GFA) are two related
topics (Lundberg et al.,, 2018). IFA attributes the prediction by a model for
a particular input to each feature, and it helps to interpret the decision-
making process for a given input. Examples of IFA methods include the Saabas
method (Saabas,, 2014) which is designed for random forests, and SHAP
(Lundberg et al.,, 2018) which is model-agnostic. On the other hand, GFA
focuses on estimating the overall importance of each feature in model
inference, and it provides a way to understand the high-level relationship
between input features and model prediction. Examples of GFA methods include
MDI (Breiman,, 2001) & MDI-oob (Li et al.,, 2019) which are designed for
random forests, and permutation importance (Breiman,, 2001) which is model-
agnostic.

A common way to construct GFA from IFA is by averaging or summing the absolute
value of the IFA for each feature across all samples. However, this approach
does not take into account the relationship between IFA and the fitting
target, which makes the resultant GFA suffer from feature selection bias —
some features tend to be over-weighted or under-weighted, and thus their GFA
is not reflective of their actual importance in the data. Many researchers
proposed ways to deal with this bias. Sandri and Zuccolotto, (2008) proposed a
bias-corrected impurity based on estimating the bias with pseudo data, for
which the R package ranger provides an efficient algorithm and implementation
(Nembrini et al.,, 2018). Some researchers ascribe the bias to over-fitting
and address it by incorporating validation data. In particular, Loecher, 2022b
and Zhou and Hooker, (2021) developed variants of MDI based on different
impurity functions calculated on both in-bag and out-of-bag data. Li et al.,
(2019) proposed MDI-oob, a debiased GFA calculated on out-of-bag data
motivated by a new characterization of MDI. There is also a separate line of
research that focuses on growing trees without such bias, e.g. cForest
(Hothorn et al.,, 2006) and honest trees (Wager and Athey,, 2018). However,
these papers focus on random forests, but there is no theoretical or empirical
evidence that they also apply to gradient-boosted trees (GBT). Additionally,
these works only apply to trees trained without regularization, whereas it is
common practice to apply $\ell_{2}$ regularization to base learners in the GBT
training process, e.g. the $\ell_{2}$ regularization parameter defaults to 1
in xgboost.

Inspired by the ideas researched for random forests, we construct a novel IFA
called PreDecomp and a family of debiased GFAs for gradient-boosted trees
called TreeInner. The major contributions of this paper are listed as follows:

  * •

We propose a novel IFA for $\ell_{2}$-regularized gradient boosted trees
called Prediction Decomposition Attribution (PreDecomp) by generalizing the
Saabas method.

  * •

We derive an alternative analytical expression for the total gain DFA in
gradient boosted trees and propose a family of debiased GFAs called TreeInner
by calculating the inner product between any IFA and the labels on an out-
sample dataset.

  * •

We show that TreeInner materialized by PreDecomp and TreeSHAP achieves state-
of-the-art AUC scores in noisy feature identification through simulation
studies inspired by real data.

###  1.1 Related Works

##### Individualized Feature Attribution (IFA)

Local surrogate models like LIME (Ribeiro et al.,, 2016) require sampling from
a neighborhood, which is especially tricky for tabular data. The Saabas method
(Saabas,, 2014) is only defined for random forests and does not take
regularization into account. SHAP (Lundberg and Lee,, 2017; Lundberg et al.,,
2018; Aas et al.,, 2019) has many ways of operationalization giving very
different results (Sundararajan and Najmi,, 2020). Specifically, using a
marginal or interventional background distribution could distort the data
manifold, using a conditional background distribution could let credit creep
between correlated features, and using an interventional conditional
background distribution provides a causal explanation but requires a
background dataset (Janzing et al.,, 2020). When it comes to implementation,
KernelSHAP has exponential computational complexity in terms of explained
features, whereas TreeSHAP, as a special case of Conditional Expectations
Shapley, violates multiple axioms of Shapley values (Sundararajan and Najmi,,
2020).

##### Global Feature Attribution (GFA)

There are a number of GFA methods developed for tree ensembles. Permutation-
based methods (Breiman,, 2001) generally have high computational costs and
require multiple replications to get an accurate estimation. They also force
the model to extrapolate to unrealistic samples when features are correlated
(Hooker et al.,, 2021). Split-improvement scores such as total gain (also
called MDI or Gini Importance in the context of random forests) and split
count (Chen and Guestrin,, 2016) bias towards variables with high cardinality
(Strobl et al.,, 2007). To address the issue, multiple variants of MDI have
been proposed to be debiased (Li et al.,, 2019) or unbiased (Zhou and Hooker,,
2021; Loecher,, 2020; Loecher, 2022b, ) where “unbiased” means non-informative
features will receive an importance score of zero in expectation. We argue
that finite-sample unbiasedness is an unnatural requirement as the ground-
truth value of feature importance is not well-defined (Grömping,, 2009): why
should an unpredictive feature receive an importance score of zero instead of
negative infinity? In fact, as we show empirically in the simulation study,
unbiasedness is of debatable desirability and assigning negative scores to
noisy features can potentially improve feature selection performance. On the
other hand, SHAP feature importance, i.e. the average absolute value of SHAP
across the samples for each feature, is merely heuristic. Additionally, it
also suffers from a strong dependence on feature cardinality in both random
forests (Loecher, 2022a, ) and gradient boosted trees (Adler and Painsky,,
2022). Apart from model-based GFAs, there also exists non-parametric
techniques that don’t require a fitted model to estimate GFA (Parr and
Wilson,, 2021; Parr et al.,, 2020).

We recently learned that there is a concurrent work by Loecher, 2022a  which
discusses a similar idea as ours. However, they focus on SHAP and random
forests, whereas we propose a whole family of GFAs parameterized by the choice
of IFA for gradient boosted trees trained with $\ell_{2}$ regularization.

###  1.2 Organization

The remainder of the paper is organized as follows. In Section 2, we introduce
a novel IFA called Prediction Decomposition Attribution (PreDecomp) and study
its theoretical properties. In Section 3, we derive an alternative expression
of total gain in terms of PreDecomp and propose TreeInner, a family of GFA
based on the formula. In Section 4, we study our method’s empirical
performance and demonstrate that it has state-of-the-art feature selection
performance. In Section 5, we conclude our paper.

##  2 INDIVIDUALIZED FEATURE ATTRIBUTIONS

For simplicity, we focus on regression trees trained with mean squared error
(MSE) loss. However, our result also applies to classification trees and Gini
index since variance and Gini index are equivalent with one-hot encoding as
shown in Li et al., (2019). Without loss of generality, we also assume all
input features are numerical.

###  2.1 Background and Notations

Assuming the feature space is $\mathbb{R}^{p}$, an individualized feature
attribution is a function on
$\mathcal{M}\times\mathbb{R}^{p}\to\mathbb{R}^{p}$, where $\mathcal{M}$ is the
set of all possible models in interest. A global feature attribution is a
function on $\mathcal{M}\to\mathbb{R}^{p}$.

Denote $\mathcal{D}_{\text{train}}=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{N}$ to
be a training dataset containing iid samples, where
$\mathbf{x}_{i}=(\mathbf{x}_{i,1},\ldots,\mathbf{x}_{i,p})\in\mathbb{R}^{p}$
is the input features and $y_{i}\in\mathbb{R}$ is the response. We may also
refer to a single specimen with $(X,Y)$ where $X=(X_{1},\ldots,X_{p})$.
Similarly, denote $\mathcal{D}_{\text{valid}}$ to be a validation dataset,
which can be the validation split reserved for hyperparameter tuning for early
stopping. We define noisy features as features that are independent of the
outcome. Non-noisy features are also called relevant features.

Gradient boosted trees (Friedman,, 2001) is an ensemble of decision trees
whose prediction is the sum prediction across all trees. Each $t$ represents a
hyper-rectangle $R_{t}$ in the feature space $\mathbb{R}^{p}$, and all samples
in $R_{t}$ share a common prediction $w_{t}$. For node $t$, we define $v(t)$
to be the feature it split on, $t^{\mathrm{left}}$ to be its left child and
$t^{\mathrm{right}}$ to be its right child. Denote $M$ to be the total number
of trees or boosting rounds. For any $m\in[M]=\\{1,2,\ldots,M\\}$, denote the
prediction of tree $m$ to be $f_{m}(X)$, and the prediction of the ensemble
consisting of the first $m$ trees to be $f_{[m]}(X)=\sum_{q\leq m}f_{q}(X)$.

Gradient boosted trees are trained sequentially and greedily: tree by tree,
and split by split. For each split $t$ in tree $m$, the objective function
consists of a training loss term and a regularization term. Particularly,
denote $L(y,\hat{y})=\frac{1}{2}(y-\hat{y})^{2}$ to be the MSE loss, $\lambda$
to be the $\ell_{2}$ penalty parameter, and the objective function is
$\mathrm{obj}_{t}(w)=\sum_{\mathbf{x}_{i}\in
R_{t}}L(y_{i},f_{[m]}(\mathbf{x}_{i}))+\frac{1}{2}w^{2}$111Unless otherwise
specified, $\sum_{\mathbf{x}_{i}\in R_{t}}\cdot$ means
$\sum_{\mathbf{x}_{i}\in R_{t}\cap\mathcal{D}_{\text{train}}}\cdot$, where
$f_{[m]}(\mathbf{x}_{i})=f_{[m-1]}(\mathbf{x}_{i})+f_{m}(\mathbf{x}_{i})$.

With second-order Taylor expansion, we can construct the gradient at sample
$X$ when growing tree $m$

$G_{i,m}=\frac{\partial
L\left(y_{i},\hat{y}_{i}\right)}{\partial\hat{y}_{i}}\bigg{|}_{\hat{y}_{i}=f_{[m-1]}(\mathbf{x}_{i})}=f_{[m-1]}(\mathbf{x}_{i})-y_{i}.$
(1)

and the hessian

$\displaystyle
H_{i,m}=\frac{\partial^{2}L\left(y_{i},\hat{y}_{i}\right)}{\partial\hat{y}_{i}^{2}}\bigg{|}_{\hat{y}_{i}=f_{[m-1]}(\mathbf{x}_{i})}=1.$
(2)

Define $p_{m}(t)$ to be the prediction at a leaf node $t$ in tree $m$. By
considering the second-order optimality condition (Chen and Guestrin,, 2016),
we know that for a leaf node $t$

$\displaystyle p_{m}(t)=$ $\displaystyle\arg\min_{w}\mathrm{obj}_{t}(w)$ (3)
$\displaystyle=$ $\displaystyle-\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}G_{i,m}}{\sum_{\mathbf{x}_{i}\in R_{t}}H_{i,m}+\lambda}$
$\displaystyle=$ $\displaystyle\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda}.$

In practice, people often consider an additional learning rate parameter
$\alpha$ for shrinkage effects, resulting in

$p_{m}(t)=\alpha\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda}.$ (4)

Throughout this paper, we assume $\lambda$ and $\alpha$ are constant for all
trees.

###  2.2 Defining $p_{m}(t)$ for Inner Nodes

The Saabas method was initially proposed in the context of random forests. In
particular, the original blog post by Saabas, (2014) breaks down the
prediction value of a tree model by tracking value changes along the
prediction path. We will explain the details later in Definition 1. A key step
in the Saabas method is to assign a prediction value $p_{m}(t)$ to all nodes
$t$ in tree $m$, and use the change of prediction values to attribute features
used by each split. Normally for tree models, only leaves are used for
prediction and thus have a well-defined $p_{m}(t)$. In this subsection, we
discuss the definition of $p_{m}(t)$ when $t$ is an inner node.

The vanilla Saabas method calculates the prediction of each inner node by
computing the average response of all training samples falling within it.

$\mathring{p}_{m}(t):=\frac{1}{|R_{t}|}\sum_{\mathbf{x}_{i}\in R_{t}}y_{i}.$
(5)

For random forests, this way of assigning values to an inner node is
equivalent to treating it as a leaf node by ignoring all of its children
nodes. Coincidentally, doing so also aligns with the idea that leaf node
predictions are to minimize the mean squared error without any regularization
term. That is to say, we have

$\arg\min_{w}\mathrm{obj}(w)\stackrel{{\scriptstyle\text{\raisebox{0.9pt}{\raisebox{-.9pt}{1}⃝}}}}{{=}}\mathring{p}_{m}(t)\stackrel{{\scriptstyle\text{\raisebox{0.9pt}{\raisebox{-.9pt}{2}⃝}}}}{{=}}\frac{1}{|R_{t}|}\sum_{\mathbf{x}_{i}\in
R_{t}}f_{m}(\mathbf{x}_{i}).$ (6)

Unfortunately, when leaf values are trained with regularization terms such as
$\ell_{2}$, the average fitted value of all training samples within a node
will be different from treating the inner node as a leaf node. In other words,
\raisebox{-.9pt}{1}⃝ and \raisebox{-.9pt}{2}⃝ in the above equation cannot
hold simultaneously. We choose to maintain \raisebox{-.9pt}{1}⃝ in Eq. (6) and
let

$\displaystyle\hat{p}_{m}(t)=$ $\displaystyle\arg\min_{w}\mathrm{obj}(w)$ (7)
$\displaystyle=$ $\displaystyle-\alpha\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}G_{i,m}}{\sum_{\mathbf{x}_{i}\in R_{t}}H_{i,m}+\lambda}$
$\displaystyle=$ $\displaystyle\alpha\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda}.$

Hereafter, we denote $p_{m}(t):=\hat{p}_{m}(t)$.

##### Alternative Definition

For the sake of completeness, we note that some software packages implement
the Saabas method with an alternative definition $\tilde{p}_{m}(t)$ for inner
nodes as the expected output of the model (Lundberg et al.,, 2020).
Specifically, $\tilde{p}_{m}(t)$ calculates the expected prediction value for
all training samples falling within an inner node by enforcing
\raisebox{-.9pt}{2}⃝ in Eq. (6), i.e.

$\tilde{p}_{m}(t)=\frac{1}{|R_{t}|}\sum_{\mathbf{x}_{i}\in
R_{t}}f_{m}(\mathbf{x}_{i}).$ (8)

While $\tilde{p}_{m}(t)$ and $\hat{p}_{m}(t)$ are asymptotically equivalent,
there is some difference in the finite-sample case. We endorse
$\hat{p}_{m}(t)$ for its nice analytical property to be demonstrated later,
and because it provides a robust tree-level bias regardless of the tree
structure. See Section 6 in the supplementary material for further discussion
and an example.

###  2.3 Generalizing Saabas Method

In this subsection, we propose a generalization of the Saabas method to
gradient boosted trees trained with $\ell^{2}$ regularization in terms of
$p_{m}(t)$.

Now we can define $f_{m,k}(X)$, the contribution of feature $k$ in the
prediction of tree $m$.

######  Definition 1 (Prediction Decomposition Attribution (PreDecomp)).

For tree $m$ and feature $k$,

$\displaystyle f_{m,k}(X)=\sum_{t\in I_{m}:v(t)=k}$
$\displaystyle\left[p_{m}(t^{\mathrm{left}})\mathds{1}(X\in
R_{t^{\mathrm{left}}})\right.$ (9)
$\displaystyle+p_{m}(t^{\mathrm{right}})\mathds{1}(X\in
R_{t^{\mathrm{right}}})$ $\displaystyle-\left.p_{m}(t)\mathds{1}(X\in
R_{t})\right].$

Additionally, for the first $m$ trees,

$f_{[m],k}(X)=\sum_{q\leq m}f_{q,k}(X)$ (10)

While PreDecomp looks superficially similar to the vanilla Sabbas method, we
note that the key difference lies in our choice of $p_{m}(t)$. The time
complexity of evaluating $f_{m,k}(X)$ is the same as doing inference with the
tree $m$.

Next, we point out that the prediction of each tree is just a sum of
$f_{m,k}(X)$ and a constant term which can be understood as a tree-level bias.

######  Proposition 1 (Local Accuracy).

For any $m\in[M]$, we have the following equation, where $p_{m}(\text{root})$
is the prediction value for the root node.

$\displaystyle f_{m}(X)=$
$\displaystyle\sum_{k=1}^{p}f_{m,k}(X)+p_{m}(\text{root})$ (11)
$\displaystyle=$
$\displaystyle\sum_{k=1}^{p}f_{m,k}(X)-\alpha\frac{\sum_{i=1}^{n}G_{i,m}}{\sum_{i=1}^{n}H_{i,m}+\lambda},$

###### Proof.

See Section 8.1 in the supplementary material. ∎

In the following theorem, we point out that our characterization can
faithfully recover additive models when features are independent.

######  Proposition 2 (Additive Model Consistency with Infinite Samples).

Assume $\mathbb{E}(Y)=0$ for identifiability, and suppose that

  * •

There are infinite training samples (population case);

  * •

$X_{1},\ldots,X_{p}$ are independent;

  * •

The regression function $f^{*}(X)=\mathbb{E}(Y|X)$ is additive, i.e., there
exists uni-variate functions $h_{1},\ldots,h_{p}$ such that
$\mathbb{E}_{X}(h_{k}(X_{k}))=0$ and for any $X\in\mathbb{R}^{p}$

$f^{*}(X)=\sum_{k=1}^{p}h_{k}(X_{k}).$ (12)

then for any $1\leq k\leq p$, the IFA $f_{[M],k}$ can consistently recover
$h_{k}$ as the number of trees $M\to\infty$:

$\lim_{M\to\infty}\mathbb{E}(f_{[M],k}(X)-h_{k}(X_{j}))^{2}=0.$ (13)

###### Proof.

See Section 8.2 in the supplementary material. Note that without the
assumptions $\mathbb{E}(Y)=0$ and $\mathbb{E}_{X}(h_{k}(X))=0$, the global
bias can be absorbed into an arbitrary $h_{k}(X)$, which leads to
unidentifiablity. ∎

Indeed, each $X_{k}$ can also be interpreted as a group of features. There can
be feature dependency and interaction within each group, as long as the groups
are mutually independent and have no interaction with each other.

##  3 GLOBAL FEATURE ATTRIBUTIONS

Based on our understanding of PreDecomp, we can construct a family of debiased
GFA. We begin with analyzing the relationship between PreDecomp and total
gain, which serves as a motivation for our construction.

###  3.1 Alternative Formula for Total Gain

In this subsection, we provide an alternative analytical expression for total
gain.

######  Definition 2 (Gain).

Gain, or the amount of decreased loss due to the split at node $t$ in tree
$m$, is defined as the sum of all the gains due to splits using feature $k$.

$\Delta_{m}(t)=\begin{aligned} &\frac{(\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}}}H_{i,m}+\lambda}+\frac{(\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}}}H_{i,m}+\lambda}\\\ -&\frac{(\sum_{\mathbf{x}_{i}\in
R_{t}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in
R_{t}}H_{i,m}+\lambda}.\end{aligned}$ (14)

Note again that calculating gain implicitly requires a background dataset
$\mathcal{D}$, and that by summing over $\mathbf{x}_{i}\in R_{t}$ we really
mean $\mathbf{x}_{i}\in R_{t}\cap\mathcal{D}$. In the original definition of
gain, we have $\mathcal{D}:=\mathcal{D}_{\text{train}}$.

Next, we define total gain for individual trees and forests.

######  Definition 3 (Total Gain).

The total gain in tree $m$ with respect to feature $k$ is defined as

$\mathrm{TotalGain}_{m,k}=\sum_{t\in I_{m}:v(t)=k}\Delta_{m}(t).$ (15)

The total gain across the whole forest with respect to feature $k$ is defined
as

$\displaystyle\mathrm{TotalGain}_{k}$
$\displaystyle=\sum_{m=1}^{M}\mathrm{TotalGain}_{m,k}$ (16)
$\displaystyle=\sum_{m=1}^{M}\sum_{t\in I_{m}:v(t)=k}\Delta_{m}(t).$

Then we have the following proposition:

######  Proposition 3.

The total gain feature attribution for tree $m$ and feature $k$ can be written
as:

$\displaystyle\mathrm{TotalGain}_{m,k}=$
$\displaystyle-\alpha^{-1}\sum_{i=1}^{N}f_{m,k}(\mathbf{x}_{i})G_{i,m}$ (17)
$\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}f_{m,k}(\mathbf{x}_{i})\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right).$

###### Proof.

See Section 8.3 in the supplementary material. ∎

As a corollary, we immediately have

$\displaystyle\mathrm{TotalGain}_{k}=$
$\displaystyle\sum_{m=1}^{M}\mathrm{TotalGain}_{m,k}$ (18) $\displaystyle=$
$\displaystyle-\alpha^{-1}\sum_{i=1}^{N}\sum_{m=1}^{M}f_{m,k}(\mathbf{x}_{i})G_{i,m}$
$\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}\sum_{m=1}^{M}f_{m,k}(\mathbf{x}_{i})\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right).$

It was pointed out by Li et al., (2019) that in random forests MDI can be
written as the in-sample covariance, in addition to the inner product, between
the response $y_{i}$ and an IFA similar to $f_{m,k}(\mathbf{x}_{i})$. To
compare our results, we make the following remark

######  Remark 1.

When $\lambda=0$, for any $k$,

$\sum_{\mathbf{x}\in R_{\text{root}}}f_{m}(\mathbf{x})=0.$ (19)

###### Proof.

See Section 8.4 in the supplementary material. ∎

When $\lambda=0$, Proposition 3 can also be interpreted as the in-sample
covariance between the gradient and $f_{m,k}(\mathbf{x}_{i})$ by Remark 1. For
$\lambda>0$, we note that it is still approximately a covariance since the
regularization term becomes negligible with a sufficiently large sample size.

###  3.2 Bridging IFA and GFA

Inspired by Proposition 3 which shows the total gain of each tree is just the
inner product of the label and the IFA for the in-sample data, we generalize
the idea to any IFA and any dataset (possibly out-sample).

######  Definition 4 (TreeInner).

For any IFA $r_{m,k}(\mathbf{x})$ defined for tree $m$ and feature $k$, we can
define a tree-level GFA as follows

$\displaystyle\operatorname{GFA}_{r}(m,k)=$
$\displaystyle-\alpha^{-1}\sum_{i=1}^{N}r_{m,k}(\mathbf{x}_{i})G_{i}$ (20)
$\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}r_{m,k}(\mathbf{x}_{i})(y_{i}-f_{[m-1]}(\mathbf{x}_{i})).$

Based on $\operatorname{GFA}_{r}(m,k)$, we can define a GFA for the entire
forest

$\displaystyle\operatorname{GFA}_{r}(k)=$
$\displaystyle\sum_{m=1}^{M}\operatorname{GFA}_{r}(m,k)$ (21) $\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}\sum_{m=1}^{M}\left(r_{m,k}(\mathbf{x}_{i})(y_{i}-f_{[m-1]}(\mathbf{x}_{i}))\right).$

Note that we can evaluate the importance of a feature on the out-sample data,
which mitigates a common feature selection bias existing for previous GFA
methods (Li et al.,, 2019; Loecher,, 2020). In other words, we sum over
$\mathbf{x}_{i}\in R_{t}\cap\mathcal{D}_{\text{valid}}$ instead of
$\mathbf{x}_{i}\in R_{t}\cap\mathcal{D}_{\text{train}}$.

##### Alternative GFA Family

We have considered another way of inducing GFA. It works by first aggregating
the tree-level IFA into a forest-level IFA, and then calculating its inner
product with the label. In particular,

######  Definition 5 (ForestInner).

$\displaystyle\widetilde{\operatorname{GFA}}_{r}(k)$ (22) $\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}\left(\sum_{m=1}^{M}r_{m,k}(\mathbf{x}_{i})\right)\left(\sum_{m=1}^{M}y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)$
$\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}\left(\sum_{m=1}^{M}r_{m,k}(\mathbf{x}_{i})\right)y_{i}.$

When using PreDecomp as the IFA, $\widetilde{\operatorname{GFA}}_{f}(k)$ can
be shown to be asymptotically unbiased for noisy features in the population
case. However, it empirically has worse performance, so we move the discussion
about it to Section 7 and the relevant experimental result to Section 9 in the
supplementary material.

##  4 SIMULATION STUDY

In this section, we conduct several numerical experiments to evaluate the
empirical performance of our theoretical results.

###  4.1 Datasets

We adopt the two synthetic datasets introduced by Li et al., (2019), each
containing a regression task and a binary classification task. The authors
replicated each experiment 40 times and report the average result, but due to
constraints on the computational budget, we only replicate the experiments 20
times. The data-generating process is described as follows. Note that this
simulation setting has a low signal-to-noise ratio by design.

##### Simulated Dataset

The first simulated dataset consists of 1000 (training) + 1000 (validation)
samples and 50 features. All features are discrete, and the $j$-th feature
takes values in $0,1,\ldots,j$ with equal probability. A set $S$ of 5 features
is randomly sampled from the first 10 features as relevant features, and the
other 45 are treated as noise. This setup represents a challenging case since
split-improvement scores bias toward features with high cardinality, whereas
the relevant features here all have low cardinality. All samples are i.i.d.
and all features are independent. The labels are generated using the following
rules:

Classification

$P(Y=1|X)=\textrm{Logistic}\left(\frac{2}{5}\sum_{j\in S}X_{j}/j-1\right).$

Regression

$\begin{aligned} Y&=\frac{1}{5}\sum_{j\in S}X_{j}/j+\epsilon\text{, where}\\\
\epsilon&\sim\mathcal{N}\left(0,100\cdot\operatorname*{Var}\left(\frac{1}{5}\sum_{j\in
S}X_{j}/j\right)\right).\end{aligned}$

##### Genomic ChIP Dataset

The second dataset is derived from a ChIP-chip and ChIP-seq dataset measuring
the enhancer status of 80 biomolecules at 3912 (training) + 3897 (validation)
regions of the Drosophila genome. The heterogeneity and dependencies among
features make it a more realistic and more challenging dataset for feature
selection problems. To evaluate feature selection in the ChIP data, all
features are min-max scaled between 0 and 1. Next, a set $S$ of 5 features is
randomly sampled as relevant features, and the other 75 are treated as noise.
Then, all noisy features are permuted to break their dependencies with
relevant features and each other, so that a noisy feature will not be
considered relevant due to its correlation with a relevant feature. The labels
are generated using the following rules:

Classification

$P(Y=1|X)=\textrm{Logistic}\left(\frac{2}{5}\sum_{j\in S}X_{j}-1\right)$.

Regression

$\begin{aligned} Y&=\frac{1}{5}\sum_{j\in S}X_{j}+\epsilon\text{, where}\\\
\epsilon&\sim\mathcal{N}\left(0,100\cdot\operatorname*{Var}\left(\frac{1}{5}\sum_{j\in
S}X_{j}\right)\right).\end{aligned}$

###  4.2 Setup

Throughout the experiments, we use a fork222Our fork will be open-sourced
after the review period. For the time being, please use the pre-built wheel
file in the supplementary material to reproduce our experiment. of xgboost
release 1.6.2 to calculate PreDecomp. The hyperparameters below are described
following the naming convention of the xgboost API.

The experiments were executed in parallel on five TPU v3-8 virtual machines
accessed through Google’s TPU Research Cloud program. Each machine is equipped
with 96 virtual cores and 340GiB main memory.

To explore the performance of each GFA family under different over-fitting
severity, we train the model with different hyperparameters as listed in Table
1. Since a grid search would be prohibitively expensive, we alter the
hyperparameters one at a time and study their marginal effect. In particular,
we assign a standard value to every parameter when it is not the variable in
interest. For example, when considering the impact of eta, we fix max_depth to
4, min_child_weight to 1, num_boost_round to 400, reg_lambda to 1, and vary
eta from $10^{-5}$ all the way to $10^{0}$. The standard values are chosen
such that the model roughly has the best predictive performance for all tasks
on the validation set.

Table 1: Hyperparameters Name | Standard | Alternatives
---|---|---
eta | $10^{-2}$ | $\begin{aligned} (&10^{-5},10^{-4},10^{-3},\\\ &10^{-2},10^{-1},10^{0})\end{aligned}$
max_depth | 4 | $(2,4,6,8,10)$
min_child_weight | 1 | $(0.5,1,2,4,8)$
num_boost_round | 400 | $\begin{aligned} (&200,400,600,\\\ &800,1000)\end{aligned}$
reg_lambda | 1 | $(0,0.1,1,10,100)$

Due to space limitations, we only show the result with varying max_depth in
the main text. Additional plots showing the result of sweeping over other
hyperparameters can be found in the supplementary material.

###  4.3 Verification of the Proposed Formula for Total Gain

To empirically verify the correctness of Proposition 3, we compare the result
obtained with Eq. (17) against the output of xgboost’s built-in method
get_score(importance_type="total_gain"). Since we are not interested in the
absolute magnitude of total gain, the result of each method are normalized
such that they sum to one, i.e. have unit $\ell_{1}$ norm in $\mathbb{R}^{p}$.
Figure 1 illustrates the maximum absolute difference between these two methods
as a measure of the deviation. The error is mostly in the order $10^{-6}$, and
it is likely caused by the inexactness of floating-point arithmetic.

Figure 1: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method.

###  4.4 Noisy Feature Identification

In this experiment, we evaluate the faithfulness of TreeInner to identify
noisy features. Labeling the noisy features as 0 and the relevant features as
1, we can evaluate the scores produced by a GFA in terms of its area under the
receiver operating characteristic curve (AUC). Ideally, the GFA values for
relevant features would be higher than those for noisy features, so we expect
to get an AUC close to 1.

We compare three GFA families: TreeInner is the debiased GFA family we
proposed; Abs is the conventional method of calculating the average magnitude
of an IFA; Permutation is the permutation importance. We evaluate the GFA on
both in-sample data and out-sample data, and each family is materialized with
two IFAs: PreDecomp and TreeSHAP. Additionally, the permutation importance is
included as a baseline. The AUC scores for all methods with a model trained
with standard hyperparameters can be found in Table LABEL:table:standard-auc
in Section 9 in the supplementary material. Note that the table includes an
additional GFA family ForestInner as defined in Eq. (22), which is omitted
here to avoid cluttering the plot.

Figure 2 illustrates the change in AUC scores as the model becomes more and
more complex with increasing max_depth. Our proposed GFA family TreeInner
dominates other GFAs in all tasks except classification on ChIP, in which case
the error bars are very wide. Note that using out-sample data generally gives
a higher AUC, which agrees with our belief that feature selection bias comes
from over-fitting. Finally, observe that when calculating TreeInner on out-
sample data, PreDecomp generally gives the same or higher AUC compared to
TreeSHAP.

Figure 2:  AUC score for noisy feature identification, averaged across 20
replications. The error bars correspond to one standard error. Higher is
better.

We attempt to explain the efficacy of TreeInner by plotting the average
normalized score for noisy features in Figure 3. The scores are normalized by
scaling them to have unit $\ell_{2}$ norm in $\mathbb{R}^{p}$. Note that
scores calculated on in-sample data are generally positive despite being
calculated for noisy features, a consequence of over-fitting. When evaluated
on out-sample data, only TreeInner can consistently produce negative scores,
hence separating the scores for noisy and relevant features. We hypothesize
that it is the lack of finite-sample unbiasedness that makes TreeInner a
better feature selector. Moreover, the trend becomes clearer as the number of
trees increases, suggesting that TreeInner could even benefit from model over-
fitting.

Figure 3:  The average normalized score for the noisy features, averaged
across 20 replications. The error bars correspond to one standard error. Lower
is better.

##  5 DISCUSSION

We proposed PreDecomp, a novel IFA for gradient boosted trees in the presence
of $\ell_{2}$ regularization, and TreeInner, a principled way to induce GFA
from any IFA and background dataset. When partnered together, PreDecomp and
TreeInner show state-of-the-art feature selection performance on two real-life
inspired synthetic datasets.

###  5.1 Limitations

In the presence of $\ell_{1}$ regularization, while one can still define
PreDecomp as $\arg\min_{w}\mathrm{obj}(w)$, we are unable to express total
gain in terms of PreDecomp in that case due to the soft thresholding effect.

Like interventional TreeSHAP, our PreDecomp also requires access to a
background dataset. Additionally, the time complexity of calculating total
gain with Eq. (17) is $O(N)$ where $N$ is the number of training samples,
whereas using the definition gives $O(1)$ complexity.

While the experiment provides some empirical insight into how the choice of
IFA could affect the performance of TreeInner, further analysis is required to
understand the relationship in detail.

###  5.2 Societal Impact

Gradient boosted trees are one of the most popular algorithms for tabular
data. As the winning solution to multiple machine learning challenges in
medicine, fraud detection, travel safety, and business, to name just a few,
they are responsible for making many high-stake decisions (Chen and Guestrin,,
2016). The IFA we proposed, i.e. PreDecomp, makes these decisions more
transparent and trustworthy by identifying the key features which lead to the
judgment. The family of GFA we proposed, i.e. TreeInner, can be used for
feature selection, which can serve as building blocks for more predictive and
robust algorithms.

Additionally, we can use gradient boosted trees as a lens to understand the
importance of features in a dataset by first fitting a model and then
calculating GFA with TreeInner.

#### Acknowledgements

This paper is a continuation of the work I did with Dr. Karl Kumbier and Dr.
Yu (Hue) Wang from the Yu Group at UC Berkeley. I thank Prof. Bin Yu and
members of the Yu Group for the inspiration and discussion. In particular, I
thank Hue for his guidance and significant contribution to this work.
Additionally, I thank Dr. Sayna Ebriahimi and Dr. Kevin Murphy from Google
Brain for their feedback. Finally, I thank Google Cloud for donating the
computational resources via the TRC program.

## References

  * Aas et al., (2019) Aas, K., Jullum, M., and Løland, A. (2019).  Explaining individual predictions when features are dependent: More accurate approximations to shapley values.
  * Adler and Painsky, (2022) Adler, A. I. and Painsky, A. (2022).  Feature importance in gradient boosting trees with cross-validation feature selection.  Entropy, 24(5).
  * Alekh, (2018) Alekh, S. (2018).  Eu general data protection regulation: A gentle introduction.
  * Basu et al., (2018) Basu, S., Kumbier, K., Brown, J. B., and Yu, B. (2018).  Iterative random forests to discover predictive and stable high-order interactions.  Proceedings of the National Academy of Sciences, 115(8):1943–1948.
  * Breiman, (2001) Breiman, L. (2001).  Random Forests.  Machine Learning, 45:1–33.
  * Caruana et al., (2015) Caruana, R., Lou, Y., Gehrke, J., Koch, P., Sturm, M., and Elhadad, N. (2015).  Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’15, page 1721–1730, New York, NY, USA. Association for Computing Machinery.
  * Chen and Guestrin, (2016) Chen, T. and Guestrin, C. (2016).  Xgboost.  Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.
  * Doshi-Velez and Kim, (2017) Doshi-Velez, F. and Kim, B. (2017).  Towards a rigorous science of interpretable machine learning.  arXiv.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  The Annals of Statistics, 29(5):1189–1232.
  * Grömping, (2009) Grömping, U. (2009).  Variable importance assessment in regression: Linear regression versus random forest.  The American Statistician, 63(4):308–319.
  * Hooker et al., (2021) Hooker, G., Mentch, L., and Zhou, S. (2021).  Unrestricted permutation forces extrapolation: variable importance requires at least one more model, or there is no free variable importance.  Statistics and Computing, 31(6):1–16.
  * Hothorn et al., (2006) Hothorn, T., Hornik, K., and Zeileis, A. (2006).  Unbiased recursive partitioning: A conditional inference framework.  Journal of Computational and Graphical Statistics, 15(3):651–674.
  * Janzing et al., (2020) Janzing, D., Minorics, L., and Bloebaum, P. (2020).  Feature relevance quantification in explainable ai: A causal problem.  In Chiappa, S. and Calandra, R., editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 2907–2916. PMLR.
  * Ke et al., (2017) Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., and Liu, T.-Y. (2017).  Lightgbm: A highly efficient gradient boosting decision tree.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Kim, (2019) Kim, J. (2019).  Explainable and Advisable Learning for Self-driving Vehicles.  PhD thesis, University of California, Berkeley, USA.
  * Kumbier et al., (2018) Kumbier, K., Basu, S., Brown, J., Celniker, S., and Yu, B. (2018).  Refining interaction search through signed iterative random forests.  bioRxiv.
  * Li et al., (2019) Li, X., Wang, Y., Basu, S., Kumbier, K., and Yu, B. (2019).  A debiased mdi feature importance measure for random forests.  Advances in Neural Information Processing Systems, 32.
  * Loecher, (2020) Loecher, M. (2020).  From unbiased mdi feature importance to explainable ai for trees.  arXiv.
  * (19) Loecher, M. (2022a).  Debiasing mdi feature importance and shap values in tree ensembles.  In Holzinger, A., Kieseberg, P., Tjoa, A. M., and Weippl, E., editors, Machine Learning and Knowledge Extraction, pages 114–129, Cham. Springer International Publishing.
  * (20) Loecher, M. (2022b).  Unbiased variable importance for random forests.  Communications in Statistics - Theory and Methods, 51(5):1413–1425.
  * Lundberg et al., (2020) Lundberg, S. M., Erion, G., Chen, H., DeGrave, A., Prutkin, J. M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., and Lee, S.-I. (2020).  From local explanations to global understanding with explainable ai for trees.  Nature machine intelligence, 2(1):56–67.
  * Lundberg et al., (2018) Lundberg, S. M., Erion, G. G., and Lee, S.-I. (2018).  Consistent Individualized Feature Attribution for Tree Ensembles.  ArXiv e-prints arXiv:1802.03888.
  * Lundberg and Lee, (2017) Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems 30, pages 4765–4774. Curran Associates, Inc.
  * Nembrini et al., (2018) Nembrini, S., König, I. R., and Wright, M. N. (2018).  The revival of the Gini importance?  Bioinformatics, 34(21):3711–3718.
  * Parr and Wilson, (2021) Parr, T. and Wilson, J. D. (2021).  Partial dependence through stratification.  Machine Learning with Applications, 6:100146.
  * Parr et al., (2020) Parr, T., Wilson, J. D., and Hamrick, J. (2020).  Nonparametric feature impact and importance.
  * Prokhorenkova et al., (2018) Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., and Gulin, A. (2018).  Catboost: unbiased boosting with categorical features.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc.
  * Ribeiro et al., (2016) Ribeiro, M. T., Singh, S., and Guestrin, C. (2016).  ”why should i trust you?”: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 1135–1144, New York, NY, USA. Association for Computing Machinery.
  * Rudin and Ustun, (2018) Rudin, C. and Ustun, B. (2018).  Optimized scoring systems: Toward trust in machine learning for healthcare and criminal justice.  Interfaces, 48:449–466.
  * Rudin et al., (2020) Rudin, C., Wang, C., and Coker, B. (2020).  The Age of Secrecy and Unfairness in Recidivism Prediction.  Harvard Data Science Review, 2(1).  https://hdsr.mitpress.mit.edu/pub/7z10o269.
  * Saabas, (2014) Saabas, A. (2014).  Interpreting random forests.
  * Sandri and Zuccolotto, (2008) Sandri, M. and Zuccolotto, P. (2008).  A bias correction algorithm for the gini variable importance measure in classification trees.  Journal of Computational and Graphical Statistics, 17(3):611–628.
  * Strobl et al., (2007) Strobl, C., Boulesteix, A.-L., Zeileis, A., and Hothorn, T. (2007).  Bias in random forest variable importance measures: Illustrations, sources and a solution.  BMC bioinformatics, 8(1):1–21.
  * Sundararajan and Najmi, (2020) Sundararajan, M. and Najmi, A. (2020).  The many shapley values for model explanation.  In III, H. D. and Singh, A., editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 9269–9278. PMLR.
  * Wager and Athey, (2018) Wager, S. and Athey, S. (2018).  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, 113(523):1228–1242.
  * Zhou and Hooker, (2021) Zhou, Z. and Hooker, G. (2021).  Unbiased measurement of feature importance in tree-based methods.  ACM Transactions on Knowledge Discovery from Data (TKDD), 15(2):1–21.

Supplementary Materials

##  6 COMPARING $\hat{p}_{m}(t)$ AND $\tilde{p}_{m}(t)$

To provide some insight into the relationship between $\hat{p}_{m}(t)$ and
$\tilde{p}_{m}(t)$, we note that both can be equivalently defined with a
recurrence formula

$\displaystyle\hat{p}_{m}(t)$
$\displaystyle=\begin{cases}\alpha\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda},&\text{$t$
is a leaf node}\\\ \begin{aligned}
&\frac{|R_{\mathrm{left}}|+\lambda}{|R_{t}|+\lambda}\hat{p}_{m}(R_{\mathrm{left}})\\\
&+\frac{|R_{\mathrm{right}}|+\lambda}{|R_{t}|+\lambda}\hat{p}_{m}(R_{\mathrm{right}})\end{aligned},&\text{$t$
is an inner node}\end{cases}.$ (23) $\displaystyle\tilde{p}_{m}(t)$
$\displaystyle=\begin{cases}\alpha\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda},&\text{$t$
is a leaf node}\\\ \begin{aligned}
&\frac{|R_{\mathrm{left}}|}{|R_{t}|}\tilde{p}_{m}(R_{\mathrm{left}})\\\
&+\frac{|R_{\mathrm{right}}|}{|R_{t}|}\tilde{p}_{m}(R_{\mathrm{right}})\end{aligned},&\text{$t$
is an inner node}\end{cases}.$

It is easy to see that $\hat{p}_{m}(t)=\tilde{p}_{m}(t)$ if and only if
$\lambda=0$. However, for any fixed tree structure, their difference converges
to zero as the sample size goes to infinity.

In xgboost and shap.TreeExplainer, a Saabas-like IFA induced by
$\tilde{p}_{m}(t)$ is used as an approximation of TreeSHAP (Chen and
Guestrin,, 2016; Lundberg et al.,, 2018). While $\tilde{p}_{m}(t)$ has its own
merit, we argue that it is better to define IFA in terms of $\hat{p}_{m}(t)$.
As we showed in Proposition 3, our $\hat{p}_{m}(t)$ has some nice analytical
relationship with total gain. Additionally, $\hat{p}_{m}(t)$ provides a robust
prediction value for the root node, as illustrated by Example 1. Moreover,
$\tilde{p}_{m}(t)$ has to be computed during inference time since the leaf
values are unknown during training, while $\hat{p}_{m}(t)$ can be computed and
cached in training time.

######  Example 1 ($\hat{p}_{m}(t)$ offers robust tree-level bias).

To further illustrate why we favor $\hat{p}_{m}(t)$ over $\tilde{p}_{m}(t)$,
consider the following artificial data set consisting of three specimens A, B,
and C:

Table 2: Sample Dataset Specimen | $X_{1}$ | $X_{2}$ | $Y$
---|---|---|---
A | 0 | 0 | 0
B | 0 | 1 | 1
C | 1 | 0 | -1

If we train a tree with $\lambda=1$ and maximum depth 1, the two possible tree
structures are shown in Figure 4. Since $X_{1}$ and $X_{2}$ are equally
predictive in this example, the tree can split on either variable, so the
training procedure would result in Structure I and Structure II with equal
possibility depending on the random seed. In practice, the instability of tree
structure could stem from perturbations in $Y$ and column sub-sampling.

baseline, for tree=draw, ellipse, align=center, s sep+=15pt, l sep+=15pt [{A,
B, C}
G = 0, H = 3 [{A, B}
G = 1, H = 2, edge label=node[midway,fill=white]$X_{1}\leq 0$] [{C}
G = -1, H = 1, edge label=node[midway,fill=white]$X_{1}>0$] ]

baseline, for tree=draw, ellipse, align=center, s sep+=15pt, l sep+=15pt [{A,
B, C}
G = 0, H = 3 [{A, C}
G = -1, H = 2, edge label=node[midway,fill=white]$X_{2}\leq 0$] [{B}
G = 1, H = 1, edge label=node[midway,fill=white]$X_{2}>0$] ]

Figure 4: Structure I (left) and Structure II (right)

We calculate $\hat{p}_{m}(t)$ and $\tilde{p}_{m}(t)$ for the nodes in both
trees in Table 3. We focus on $\hat{p}_{m}(\text{root})$, since it is the
constant term in Proposition 1, which means it can be interpreted as a tree-
level bias. Observe that when using $\hat{p}_{m}(t)$, the prediction value of
the root node is always 0 regardless of the tree structure, whereas
$\tilde{p}_{m}(t)$ gives less stable result by alternating between
$-\frac{1}{18}$ and $\frac{1}{18}$. Indeed, this is always true for
$\hat{p}_{m}(t)$, as we have
$\hat{p}_{m}(\text{root})=-\alpha\frac{\sum_{i}G_{i,m}}{\sum_{i}H_{i,m}+\lambda}$
regardless of the tree hyperparameters and structure.

Table 3: Values of $p_{m}(t)$ Definition | $\hat{p}_{m}(t)$ | $\tilde{p}_{m}(t)$
---|---|---
Structure | Structure I | Structure II | Structure I | Structure II
root | 0 | 0 | $\frac{1}{18}$ | $-\frac{1}{18}$
left | $\frac{1}{3}$ | $-\frac{1}{3}$ | $\frac{1}{3}$ | $-\frac{1}{3}$
right | $-\frac{1}{2}$ | $\frac{1}{2}$ | $-\frac{1}{2}$ | $\frac{1}{2}$

##  7 ALTERNATIVE GFA FAMILY

We have the following proposition for the
$\widetilde{\operatorname{GFA}}_{f}(j)$ proposed in Eq. (22).

######  Proposition 4 (Asymptotic Unbiasedness).

Without loss of generality, we assume $X_{1},\ldots,X_{k}$ are relevant
features, and $X_{k+1},\ldots,X_{p}$ are noisy features. Suppose that

  * •

We are using PreDecomp as the IFA, i.e.
$r_{m,j}(\mathbf{x}_{i})=f_{m,j}(\mathbf{x}_{i})$;

  * •

There are infinite training samples (population case);

  * •

$(X_{1},\ldots,X_{k}),X_{k+1},\ldots,X_{p}$ are mutually independent, i.e.

$\displaystyle P(X_{1},X_{2},\ldots,X_{k},X_{k+1},\ldots,X_{p})$ (24)
$\displaystyle=$ $\displaystyle P(X_{1},X_{2},\ldots,X_{k})P(X_{k+1})\ldots
P(X_{p})$

Then for any noisy feature $k+1\leq j\leq p$, we have

$\lim_{M\to\infty}\widetilde{\operatorname{GFA}}_{f}(j)=0.$ (25)

###### Proof.

We consider the subset of relevant features collectively as a super-feature
$X_{S}=(X_{1},X_{2},\ldots,X_{k})$. The regression function has the form

$\displaystyle f^{*}(X)=$ $\displaystyle\mathbb{E}(Y|X)$ (26) $\displaystyle=$
$\displaystyle\mathbb{E}(Y|X_{S},X_{k+1},\ldots,X_{p})$ $\displaystyle=$
$\displaystyle\mathbb{E}(Y|X_{S})+\underbrace{0+\ldots+0}_{p-k}$

Note that this is an additive model, so by Proposition 2, for $k<j\leq p$ we
have

$\displaystyle\mathbb{E}(f_{[M],j}(X)-0)^{2}\to$ $\displaystyle 0.$ (27)
$\displaystyle f_{[M],j}(X)\to$ $\displaystyle 0.$

As a result,

$\displaystyle\widetilde{\operatorname{GFA}}_{f}(j)=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}\left(\sum_{m=1}^{M}f_{m,j}(\mathbf{x}_{i})\right)y_{i}$
(28) $\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{i=1}^{N}f_{[M],j}(\mathbf{x}_{i})y_{i}$
$\displaystyle\to$ $\displaystyle\alpha^{-1}\sum_{i=1}^{N}0\cdot y_{i}$
$\displaystyle=$ $\displaystyle 0$

∎

As a corollary, we have shown that total gain is asymptotically unbiased in
the population case. This makes sense since the bias in split-improvement
scores is essentially due to the model struggling to tell signal from noise
and overfits to the latter, but access to infinite samples means there is
effectively no noise.

##  8 PROOFS

###  8.1 Proof for Proposition 1

###### Proof.

$\displaystyle\sum_{k=1}^{p}f_{m,k}(X)=$
$\displaystyle\sum_{k=1}^{p}\sum_{t\in I_{m}:v(t)=k}\begin{aligned}
&\left[p_{m}(t^{\mathrm{left}})\mathds{1}(X\in
R_{t^{\mathrm{left}}})\right.\\\ &+p_{m}(t^{\mathrm{right}})\mathds{1}(X\in
R_{t^{\mathrm{right}}})\\\ &-\left.p_{m}(t)\mathds{1}(X\in
R_{t})\right]\end{aligned}$ (29) $\displaystyle=$ $\displaystyle\sum_{t\in
I_{m}}\begin{aligned} &\left[p_{m}(t^{\mathrm{left}})\mathds{1}(X\in
R_{t^{\mathrm{left}}})\right.\\\ &+p_{m}(t^{\mathrm{right}})\mathds{1}(X\in
R_{t^{\mathrm{right}}})\\\ &-\left.p_{m}(t)\mathds{1}(X\in
R_{t})\right]\end{aligned}$

Since all inner nodes $t\in I_{m}$ except the root node is also the child of
another node, the term $p_{m}(t)\mathds{1}(\mathbf{x}\in R_{t})$ will get
canceled by the corresponding term of its parent. The only remaining terms are
the leaf nodes and the root node. Denote $J_{m}$ to be the set of all leaf
nodes in tree $m$, and then

$\displaystyle\sum_{j=1}^{p}f_{m,j}(X)$ $\displaystyle=\sum_{t\in
J_{m}}p_{m}(t)\mathds{1}(X\in R_{t})-p_{m}(\text{root})$ (30) $\displaystyle
f_{m}(X)+p_{m}(\text{root})$ $\displaystyle=\sum_{t\in
J_{m}}p_{m}(t)\mathds{1}(X\in R_{t})$ $\displaystyle=f_{m}(X).$

That completes the proof. ∎

###  8.2 Proof for Proposition 2

###### Proof.

First of all, we observe that with an infinite sample size and the number of
trees approaching infinity, i.e., $M\to\infty$, the GBT is consistent,

$\lim_{M\to\infty}\mathbb{E}(f_{[M]}(X)-f^{*}(X))^{2}=0,$ (31)

This is because gradient boosted trees are guaranteed to approach zero mean
squared error for the training set and with infinite samples, i.e. when the
training set is the same as the population.

Now we would like to show that for tree $m$, $f_{m,k}(X)$ is only a function
of $X_{k}$. In other words, with slight abuse of notations, we have
$f_{m,k}(X)=f_{m,k}(X_{k})$. We prove this by induction. When $m=0$,
$f_{[m],k}(X)=0$ which holds trivially. Suppose the conclusion holds for $m$:
$f_{[m],k}(X)=f_{[m],k}(X_{k})$. For tree $m+1$, to show $f_{m+1,k}(X)$ only
depends on $X_{k}$, by Definition 1, we only need to show that for each node
$t$ and its parent $t_{\mathrm{parent}}$, if $t_{\mathrm{parent}}$ splits on
feature $j=v(t_{\mathrm{parent}})$, then
$p_{m+1}(t)-p_{m+1}(t_{\mathrm{parent}})$ is a uni-variate function of
$X_{j}$. Recall that in Eq. (3), $p_{m}(t)=\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|+\lambda}$, but we
have an infinite number of samples, so the lambda term in the denominator can
be ignored, resulting in $p_{m+1}(t)=\frac{\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)}{|R_{t}|}=\mathbb{E}(f^{*}(X)-f_{[m]}(X)|X\in
R_{t})$. Since $\mathbb{E}(Y)=0$, we know $p_{m}(\text{root})=0$. By
Proposition 1 and our induction assumption, $f_{[m]}(X)$ is additive, i.e.

$f_{[m]}(X)=\sum_{k=1}^{p}f_{[m],k}(X)+0=\sum_{k=1}^{p}f_{[m],k}(X_{k}).$ (32)

Given the assumption that $f^{*}(X)$ is additive, we also have

$f^{*}(X)=\sum_{k=1}^{p}h_{k}(X_{k}).$ (33)

To simplify the following deduction, let’s assume
$R_{t}=\mathds{1}_{X_{1}<t_{1},X_{2}<t_{2}}$. The general case follows the
same logic. Then we have

$\displaystyle\mathbb{E}(f^{*}(X)-f_{[m]}(X)|X_{1}<t_{1},X_{2}<t_{2})$ (34)
$\displaystyle=$
$\displaystyle\sum_{k=1}^{p}\mathbb{E}(h_{k}(X_{k})-f_{[m],k}(X_{k})|X_{1}<t_{1},X_{2}<t_{2}).$

By the independence of $X_{1},\ldots,X_{p}$, the above equation equals to

$\displaystyle\mathbb{E}(h_{1}(X_{1})-f_{[m],1}(X_{1})|X_{1}<t_{1})$ (35)
$\displaystyle+\mathbb{E}(h_{2}(X_{2})-f_{[m],2}(X_{2})|X_{2}<t_{2})$
$\displaystyle+\underbrace{\sum_{k\geq
3}\mathbb{E}(h_{k}(X_{k})-f_{[m],k}(X_{k})}_{\mathrm{const}}.$

Without loss of generality, let’s assume $t_{\mathrm{parent}}$ is
$\mathds{1}_{X_{1}<t_{1}}$. Following the same deduction, we have

$\displaystyle
p_{m+1}(t_{\mathrm{parent}})=\mathbb{E}(h_{1}(X_{1})-f_{[m],1}(X_{1})|X_{1}<t_{1})+\mathrm{const}.$
(36)

Therefore, we have
$p_{m+1}(t)-p_{m+1}(t_{\mathrm{parent}})=\mathbb{E}(h_{2}(X_{2})-f_{[m],2}(X_{2})|X_{2}<t_{2})$
which is a uni-variable function of $X_{2}$. That concludes our proof that
$f_{m,k}(X)$ only depends on $X_{k}$. By induction, we know
$f_{[m],k}(X)=f_{[m],k}(X_{k})$ for any $k=1,\ldots,p$.

Finally, we are ready to show Eq. (13). Since $X_{j}$ are independent and
$f_{[m],k}(X)=f_{[m],k}(X_{k})$, we know the following equality:

$\mathbb{E}(f_{[M]}(X)-f^{*}(X))^{2}=\sum_{k=1}^{p}\mathbb{E}(f_{[M],k}(X_{k})-h_{k}(X_{k}))^{2}.$
(37)

Since the left-hand side goes to zero as $M\to\infty$, the right-hand side
must go to zero as well. Therefore, for any $k$,

$\mathbb{E}(f_{[M],k}(X_{k})-h_{k}(X_{k}))^{2}\to 0.$ (38)

That completes the proof. ∎

###  8.3 Proof for Proposition 3

###### Proof.

Recall the definition of $p_{m}(t)$ in Eq. (7).

$\displaystyle\alpha\Delta(t)=$
$\displaystyle\alpha\frac{(\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}}}H_{i,m}+\lambda}+\alpha\frac{(\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}}}H_{i,m}+\lambda}$ (39)
$\displaystyle-\alpha\frac{(\sum_{\mathbf{x}_{i}\in
R_{t}}G_{i,m})^{2}}{\sum_{\mathbf{x}_{i}\in R_{t}}H_{i,m}+\lambda}$
$\displaystyle=$ $\displaystyle-
p_{m}(t^{\mathrm{left}})\sum_{i=1}^{N}G_{i,m}\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}})$ $\displaystyle-
p_{m}(t^{\mathrm{right}})\sum_{i=1}^{N}G_{i,m}\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}})$
$\displaystyle+p_{m}(t)\sum_{i=1}^{N}G_{i,m}\mathds{1}(\mathbf{x}_{i}\in
R_{t})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}-G_{i,m}\cdot\begin{aligned}
&\left[p_{m}(t^{\mathrm{left}})\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}})\right.\\\
&+p_{m}(t^{\mathrm{right}})\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}})\\\ &\left.-p_{m}(t)\mathds{1}(\mathbf{x}_{i}\in
R_{t})\right]\end{aligned}.$

Then, by the definition of total gain,

$\displaystyle\mathrm{TotalGain}_{m,k}=$ $\displaystyle\sum_{t\in
I_{m}:v(t)=k}\Delta(t)$ (40) $\displaystyle=$
$\displaystyle\alpha^{-1}\sum_{m=1}^{M}\sum_{t\in
I_{m}:v(t)=k}\sum_{i=1}^{N}-G_{i,m}$ $\displaystyle\times\begin{aligned}
&\left[p_{m}(t^{\mathrm{left}})\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{left}}})\right.\\\
&+p_{m}(t^{\mathrm{right}})\mathds{1}(\mathbf{x}_{i}\in
R_{t^{\mathrm{right}}})\\\ &-\left.p_{m}(t)\mathds{1}(\mathbf{x}_{i}\in
R_{t})\right]\end{aligned}$ $\displaystyle=$
$\displaystyle-\alpha^{-1}\sum_{i=1}^{N}f_{m,k}(\mathbf{x}_{i})G_{i,m}.$

That completes the proof. ∎

###  8.4 Proof for Remark 1

###### Proof.

By a similar argument as the proof for Proposition 1, we have

$\displaystyle\sum_{\mathbf{x}\in R_{\text{root}}}f_{m}(\mathbf{x})$ (41)
$\displaystyle=$ $\displaystyle\sum_{\mathbf{x}\in R_{\text{root}}}\sum_{t\in
I_{m}:v(t)=k}\begin{aligned}
&\left[p_{m}(t^{\mathrm{left}})\mathds{1}(\mathbf{x}\in
R_{t^{\mathrm{left}}})\right.\\\
&+p_{m}(t^{\mathrm{right}})\mathds{1}(\mathbf{x}\in R_{t^{\mathrm{right}}})\\\
&-\left.p_{m}(t)\mathds{1}(\mathbf{x}\in R_{t})\right]\\\ \end{aligned}$
$\displaystyle=$ $\displaystyle\sum_{t\in I_{m}:v(t)=k}\begin{aligned}
&\left(p_{m}(t^{\mathrm{left}})|R_{t^{\mathrm{left}}}|\right.\\\
&+p_{m}(t^{\mathrm{right}})|R_{t^{\mathrm{right}}}|\\\
&-\left.p_{m}(t)|R_{t}|\right)\\\ \end{aligned}$ $\displaystyle=$
$\displaystyle\sum_{t\in I_{m}:v(t)=k}\begin{aligned}
&\left(\sum_{\mathbf{x}_{i}\in
R_{\mathrm{left}}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)\right.\\\
&+\sum_{\mathbf{x}_{i}\in
R_{\mathrm{right}}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)\\\
&-\left.\sum_{\mathbf{x}_{i}\in
R_{t}}\left(y_{i}-f_{[m-1]}(\mathbf{x}_{i})\right)\right)\\\ \end{aligned}$
$\displaystyle=$ $\displaystyle\sum_{t\in I_{m}:v(t)=k}0$ $\displaystyle=$
$\displaystyle 0$

∎

##  9 ADDITIONAL EXPERIMENTAL RESULTS

In this section, we provide experimental results which do not fit within the
main text.

###  9.1 AUC Scores

Table 4: AUC Scores for Model Trained with Standard Hyperparameters |  |  |  |  | AUC_mean | AUC_std | risk_mean | risk_std
---|---|---|---|---|---|---|---|---
Dataset | Task | GFA | Domain | IFA |  |  |  |
ChIP | Classification | Abs | Train | PreDecomp | 0.8581 | 0.1200 | 0.4673 | 0.0168
TreeSHAP | 0.8651 | 0.1181 | 0.4673 | 0.0168
Valid | PreDecomp | 0.8599 | 0.1172 | 0.4673 | 0.0168
TreeSHAP | 0.8643 | 0.1190 | 0.4673 | 0.0168
ForestInner | Train | PreDecomp | 0.7953 | 0.1207 | 0.4673 | 0.0168
TreeSHAP | 0.8640 | 0.1060 | 0.4673 | 0.0168
Valid | PreDecomp | 0.7975 | 0.1619 | 0.4673 | 0.0168
TreeSHAP | 0.7553 | 0.1727 | 0.4673 | 0.0168
Permutation | Train | Permutation | 0.8367 | 0.1091 | 0.4673 | 0.0168
Valid | Permutation | 0.8275 | 0.1245 | 0.4673 | 0.0168
TreeInner | Train | PreDecomp | 0.8043 | 0.1213 | 0.4673 | 0.0168
TreeSHAP | 0.8645 | 0.1081 | 0.4673 | 0.0168
Valid | PreDecomp | 0.8671 | 0.1071 | 0.4673 | 0.0168
TreeSHAP | 0.8756 | 0.1021 | 0.4673 | 0.0168
Regression | Abs | Train | PreDecomp | 0.9399 | 0.0681 | 0.1504 | 0.0867
TreeSHAP | 0.9413 | 0.0688 | 0.1504 | 0.0867
Valid | PreDecomp | 0.9383 | 0.0693 | 0.1504 | 0.0867
TreeSHAP | 0.9405 | 0.0696 | 0.1504 | 0.0867
ForestInner | Train | PreDecomp | 0.9172 | 0.0663 | 0.1504 | 0.0867
TreeSHAP | 0.9461 | 0.0651 | 0.1504 | 0.0867
Valid | PreDecomp | 0.9229 | 0.0917 | 0.1504 | 0.0867
TreeSHAP | 0.8867 | 0.0837 | 0.1504 | 0.0867
Permutation | Train | Permutation | 0.9369 | 0.0725 | 0.1504 | 0.0867
Valid | Permutation | 0.9472 | 0.0587 | 0.1504 | 0.0867
TreeInner | Train | PreDecomp | 0.9235 | 0.0694 | 0.1504 | 0.0867
TreeSHAP | 0.9472 | 0.0661 | 0.1504 | 0.0867
Valid | PreDecomp | 0.9604 | 0.0618 | 0.1504 | 0.0867
TreeSHAP | 0.9528 | 0.0734 | 0.1504 | 0.0867
Simulated | Classification | Abs | Train | PreDecomp | 0.6564 | 0.1334 | 0.4745 | 0.0191
TreeSHAP | 0.6749 | 0.1356 | 0.4745 | 0.0191
Valid | PreDecomp | 0.6500 | 0.1364 | 0.4745 | 0.0191
TreeSHAP | 0.6731 | 0.1361 | 0.4745 | 0.0191
ForestInner | Train | PreDecomp | 0.4876 | 0.1379 | 0.4745 | 0.0191
TreeSHAP | 0.6764 | 0.1337 | 0.4745 | 0.0191
Valid | PreDecomp | 0.7722 | 0.1254 | 0.4745 | 0.0191
TreeSHAP | 0.7191 | 0.1438 | 0.4745 | 0.0191
Permutation | Train | Permutation | 0.6712 | 0.1386 | 0.4745 | 0.0191
Valid | Permutation | 0.6599 | 0.1097 | 0.4745 | 0.0191
TreeInner | Train | PreDecomp | 0.5091 | 0.1429 | 0.4745 | 0.0191
TreeSHAP | 0.6798 | 0.1356 | 0.4745 | 0.0191
Valid | PreDecomp | 0.7856 | 0.1277 | 0.4745 | 0.0191
TreeSHAP | 0.7520 | 0.1541 | 0.4745 | 0.0191
Regression | Abs | Train | PreDecomp | 0.3533 | 0.1853 | 6.8457 | 1.6117
TreeSHAP | 0.3776 | 0.1983 | 6.8457 | 1.6117
Valid | PreDecomp | 0.3513 | 0.1849 | 6.8457 | 1.6117
TreeSHAP | 0.3773 | 0.1988 | 6.8457 | 1.6117
ForestInner | Train | PreDecomp | 0.2147 | 0.1183 | 6.8457 | 1.6117
TreeSHAP | 0.3651 | 0.1993 | 6.8457 | 1.6117
Valid | PreDecomp | 0.6380 | 0.0976 | 6.8457 | 1.6117
TreeSHAP | 0.6429 | 0.0965 | 6.8457 | 1.6117
Permutation | Train | Permutation | 0.3769 | 0.2002 | 6.8457 | 1.6117
Valid | Permutation | 0.5493 | 0.1392 | 6.8457 | 1.6117
TreeInner | Train | PreDecomp | 0.2278 | 0.1272 | 6.8457 | 1.6117
TreeSHAP | 0.3673 | 0.1989 | 6.8457 | 1.6117
Valid | PreDecomp | 0.6384 | 0.0953 | 6.8457 | 1.6117
TreeSHAP | 0.6496 | 0.1015 | 6.8457 | 1.6117
Table 4: AUC Scores for Model Trained with Standard Hyperparameters

###  9.2 Sweeping over eta

Figure 5: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method. Figure 6:
AUC score for noisy feature identification, averaged across 20 replications.
The error bars correspond to one standard error. Higher is better.  Figure 7:
The average normalized score for the noisy features, averaged across 20
replications. The error bars correspond to one standard error. Lower is
better.  Figure 8:  The average normalized score for the relevant features,
averaged across 20 replications. The error bars correspond to one standard
error. Lower is better.  Figure 9:  Risk evaluated on the validation set. The
error bars correspond to one standard error. Lower is better.

###  9.3 Sweeping over max_depth

Figure 10: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method. Figure 11:
AUC score for noisy feature identification, averaged across 20 replications.
The error bars correspond to one standard error. Higher is better.  Figure 12:
The average normalized score for the noisy features, averaged across 20
replications. The error bars correspond to one standard error. Lower is
better.  Figure 13:  The average normalized score for the relevant features,
averaged across 20 replications. The error bars correspond to one standard
error. Lower is better.  Figure 14:  Risk evaluated on the validation set. The
error bars correspond to one standard error. Lower is better.

###  9.4 Sweeping over min_child_weight

Figure 15: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method. Figure 16:
AUC score for noisy feature identification, averaged across 20 replications.
The error bars correspond to one standard error. Higher is better.  Figure 17:
The average normalized score for the noisy features, averaged across 20
replications. The error bars correspond to one standard error. Lower is
better.  Figure 18:  The average normalized score for the relevant features,
averaged across 20 replications. The error bars correspond to one standard
error. Lower is better.  Figure 19:  Risk evaluated on the validation set. The
error bars correspond to one standard error. Lower is better.

###  9.5 Sweeping over num_boost_round

Figure 20: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method. Figure 21:
AUC score for noisy feature identification, averaged across 20 replications.
The error bars correspond to one standard error. Higher is better.  Figure 22:
The average normalized score for the noisy features, averaged across 20
replications. The error bars correspond to one standard error. Lower is
better.  Figure 23:  The average normalized score for the relevant features,
averaged across 20 replications. The error bars correspond to one standard
error. Lower is better.  Figure 24:  Risk evaluated on the validation set. The
error bars correspond to one standard error. Lower is better.

###  9.6 Sweeping over reg_lambda

Figure 25: Logarithm of the maximum absolute difference between the normalized
total gain calculated with our methodology and the built-in method. Figure 26:
AUC score for noisy feature identification, averaged across 20 replications.
The error bars correspond to one standard error. Higher is better.  Figure 27:
The average normalized score for the noisy features, averaged across 20
replications. The error bars correspond to one standard error. Lower is
better.  Figure 28:  The average normalized score for the relevant features,
averaged across 20 replications. The error bars correspond to one standard
error. Lower is better.  Figure 29:  Risk evaluated on the validation set. The
error bars correspond to one standard error. Lower is better.