# Stable Update of Regression Trees

Morten Blørstad &Berent Å. S. Lunde &Nello Blaser 11footnotemark: 1 Department
of Informatics, University of Bergen, 5006 Bergen, Norway.Equinor, 5254
Sandsli, Norway.

###### Abstract

Updating machine learning models with new information usually improves their
predictive performance, yet, in many applications, it is also desirable to
avoid changing the model predictions too much. This property is called
stability. In most cases when stability matters, so does explainability. We
therefore focus on the stability of an inherently explainable machine learning
method, namely regression trees. We aim to use the notion of empirical
stability and design algorithms for updating regression trees that provide a
way to balance between predictability and empirical stability. To achieve
this, we propose a regularization method, where data points are weighted based
on the uncertainty in the initial model. The balance between predictability
and empirical stability can be adjusted through hyperparameters. This
regularization method is evaluated in terms of loss and stability and assessed
on a broad range of data characteristics. The results show that the proposed
update method improves stability while achieving similar or better predictive
performance. This shows that it is possible to achieve both predictive and
stable results when updating regression trees. 111A preliminary version of
results in this paper appeared in the Master thesis of the first author
(Blørstad, 2023).

##  1 Introduction

Training an initial machine learning (ML) model is just the start of the ML
life cycle. The model can be improved as more data becomes available, leading
to an iterative process of continuous improvement. Thus, in industry settings,
ML models are frequently updated to incorporate new information. The main
objective of such a process has traditionally been to improve the model’s
predictive capability. However, many industries, for instance in finance or
medicine, have additional requirements such as interpretability and stability.
Interpretability is important due to ”algorithmic accountability” (Kaminski &
Malgieri (2021)), whereas stability is crucial to ensure consistent and
reliable predictions on new unseen data. In this manuscript, we focus on the
stability aspect.

Achieving stability involves minimizing model variance as stability is
connected to generalization and inversely related to variance (Bousquet &
Elisseeff (2002)). There have been numerous studies on the stability of
learning algorithms (see, e.g., Devroye & Wagner (1979); Kearns & Ron (1999);
Bousquet & Elisseeff (2002)) focusing on theoretical bounds of their
generalization error. Some previous studies have also presented empirical
stability measures for both classification and regression tasks, which can be
divided into structural and semantic stability (Turney (1995); Briand et al.
(2009); Lin et al. (2002); Lim & Yu (2013); Philipp et al. (2017)). We focus
on semantic stability, (dis-)similarity of the predictions of two models,
which we review in Section 2. While _measuring_ stability has been studied
thoroughly, there are also some studies focusing on improving it for
classifiers. The studies by Milani Fard et al. (2016), Goh et al. (2016),
Cotter et al. (2018) and Liu et al. (2022) address the issue of unnecessary
prediction change when training consecutive classifiers. Milani Fard et al.
(2016) propose a Markov Chain Monte Carlo stabilization method, whereas Goh et
al. (2016), Cotter et al. (2018) introduces different optimization constraints
when training classifiers. Liu et al. (2022) improve model stability of deep-
learning classifiers in the context of NLP systems. In this article, we aim to
improve the empirical stability of regression trees.

Regression trees are appealing due to their inherently interpretable
structure, and their ability to discover interaction effects and select
features. However, due to regression trees’ high variance, small changes in
data can often lead to very different splits, which can cause inconsistent
predictions, i.e., instability (Breiman (1996), Hastie et al. (2001)). This
issue of instability in tree models has been addressed in previous research
for decision trees. For instance, Last et al. (2002) showed using Info-Fuzzy
Networks that semantically stable decision trees could be learned for
classification tasks while maintaining a reasonable level of accuracy.
Additionally, Dannegger (2000) proposed a bootstrap-based node-level
stabilizing approach for decision trees. This method was later improved by
Briand et al. (2009), demonstrating improved structural stability compared to
the classical decision tree algorithm (CART). To the best of our knowledge, no
work focuses on improving the empirical semantic stability of regression trees
in the context of continuous updates.

Contributions Our contributions to stabilizing regression trees are as
follows.

  1. 1.

We incorporate stability measures as a regularization term in the loss
function, providing a way of ensuring a level of stability by tuning
regularization strength. Our novel approach involves dynamically weighting the
regularization for each data point based on the model’s uncertainty.

  2. 2.

We evaluate our approach on several different datasets, showing that this
technique in combination with a base stability regularization achieves a
better balance of the predictability-stability trade-off.

  3. 3.

Our python-package for stable regression tree updates is available at
https://github.com/MortenBlorstad/StableTreeUpdates.

The remainder of the paper unfolds as follows: Section 2 introduces empirical
stability and Section 3 introduces regression trees. Section 4 describes the
proposed method for stable updates of regression trees. Section 5 provides the
experimental setup and results on various datasets. Finally, in Section 6 we
discuss the results and conclude.

##  2 Empirical Stability

Consider the regression task of learning a function (or model) $f$ that maps
the input vector $\mathbf{x}\in\mathbb{R}^{m}$ to the corresponding response
variable $y\in\mathbb{R}$. Learning such a function involves approximating the
optimal function $f^{*}$ from the set of the candidate function $\mathcal{F}$
by minimising the empirical loss given a loss function
$\mathcal{L}(\cdot,\cdot)$, over a training dataset
$\mathcal{D}=\\{\mathbf{x}_{i},y_{i}\\}_{i=1}^{n}$

$f^{*}\approx\hat{f}=\operatorname*{arg\,min}_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}[\mathcal{L}(y_{i},f(\mathbf{x}_{i}))].$

In this context, stability is defined as how much a model’s predictions change
due to changes in the training data, such as the removal or addition of a data
point. It reflects a data point’s impact on the learning process. Most studies
focus on theoretical stability guarantees, but for many ML algorithms, no such
theoretical guarantees exist. We instead take an empirical approach and
consider empirical stability measures that can be assessed using test data.
Similar to how we minimize empirical loss to improve predictability, we
introduce an empirical instability measure that we minimize to improve
stability.

Consider a trained model $f_{0}$ at time $t=0$ with the training data
$\mathcal{D}_{0}$. At $t=1$ more data, $\mathcal{D}_{\Delta t}$ becomes
available. The model is updated using all available data
$\mathcal{D}_{1}=\mathcal{D}_{0}\cup\mathcal{D}_{\Delta t}$ to obtain the
updated model $f_{1}$. The instability of the update is given by the
instability function $\mathcal{S}(\cdot,\cdot)$,

$\text{empirical
instability}=\frac{1}{n}\sum_{i=1}^{n}[\mathcal{S}(f_{0}(\mathbf{x}_{i}),f_{1}(\mathbf{x}_{i}))],$

where $\mathcal{S}(\cdot,\cdot)$ measures the discrepancy between the
predictions of $f_{0}$ and $f_{1}$ over a test dateset. In other words, we are
considering the mean pointwise empirical instability. The choice of
$\mathcal{S}(\cdot,\cdot)$, similar to the choice of
$\mathcal{L}(\cdot,\cdot)$, depends on your assumptions on the error
distribution. Assuming Gaussian and Laplace distributions lead to the squared
error $\mathcal{S}(a,b)=(a-b)^{2}$ and absolute error
$\mathcal{S}(a,b)=|a-b|$, respectively. Another empirical instability measure
that can be written as a mean pointwise empirical instability is the negative
coverage probability $\mathcal{S}(a,b,k)=-\mathbb{1}_{|a-b|\leq k}$ (Lin et
al. (2002)).

##  3 Learning Regression Trees

There are many methods of fitting Regression trees, most notably the
Classification And Regression Tree (CART) algorithm by Breiman et al. (1984)
(see Figure 1). The CART algorithm uses a top-down greedy approach to
partition the feature space $R=\mathbb{R}^{m}$ into $D$ non-overlapping
regions, $R_{1},\ldots,R_{d},\ldots,R_{D}$ and fit a leaf prediction,
$\hat{w}_{d}$ to each region based on the response observations within the
region. Starting from a root node $R_{root}$ containing all the data, the
algorithm seeks the splitting variable $j\in{1,\ldots,m}$ and split value $s$
that creates the half-planes $R_{l}(j,s)=\\{x\in R_{root}|x_{j}\leq s\\}$ and
$R_{r}(j,s)=\\{x\in R_{root}|x_{j}>s\\}$ that maximize the loss reduction
$\mathcal{R}$,

$\mathcal{R}(j,s)=\sum_{i:x_{i}\in
R_{root}}{\mathcal{L}(y_{i},\hat{w}_{root})}-\sum_{i:x_{i}\in
R_{l}(j,s)}{\mathcal{L}(y_{i},\hat{w}_{l})}-\sum_{i:x_{i}\in
R_{r}(j,s)}{\mathcal{L}(y_{i},\hat{w}_{r})}.$

Here $R_{l}$ and $R_{r}$ are the datasets of left and right nodes after the
split and $R_{root}=R_{l}\cup R_{r}$ is the dataset before the split. The
terms $\hat{w}_{root}$, $\hat{w}_{l}$ and $\hat{w}_{r}$ denote the predictions
of the root, left and right node. After finding the partition with the largest
loss reduction, the splitting process is repeated on the resulting regions
until reaching a stopping criterion. Given a learned tree model $f$, the
feature mapping $q(\mathbf{x})$ is a function that maps inputs to their
corresponding region $R$ and its prediction for input $x$ is given as
$f(\mathbf{x})=\hat{w}_{q(\mathbf{x})}$.

1 4 5 2 3 $\mathcal{V}\ =\\{1,2,3,4,5\\}$$\hat{w}_{1}=1.1$$q(\mathbf{x})=4;\
\mathbf{x}\in R_{4}$$q(\mathbf{x})=3;\ \mathbf{x}\in R_{3}$$q(\mathbf{x})=1;\
\mathbf{x}\in R_{1}$$R_{1}$$R_{3}$$R_{2}$$R_{4}$$R_{5}$TrueFalse Figure 1: An
example of a Classification and Regression Tree (CART) with five leaf nodes
$\mathcal{V}$. The vector
$\hat{\mathbf{w}}=(\hat{w}_{1},\hat{w}_{2},\hat{w}_{3},\hat{w}_{4},\hat{w}_{5})$
is the possible predictions the tree can make.

##  4 Stable Regression Tree Update

We propose an update method that incorporates a stability regularization term
into the loss function. The key lies in the determination of the
regularization strength, which we propose to be a composite of base
regularization and uncertainty-weighted regularization, each adjusted through
a hyperparameter. This section first describes the modification made to the
CART algorithm followed by a description of our proposed update method.

We modify the CART by using the second-order approximation of the loss
reduction,

$\tilde{\mathcal{R}}(j,s)=\frac{1}{2}\left[\frac{(\sum_{i:\mathbf{x}_{i}\in
R_{l}}g_{i})^{2}}{\sum_{i:\mathbf{x}_{i}\in
R_{l}}h_{i}}+\frac{(\sum_{i:\mathbf{x}_{i}\in
R_{r}}g_{i})^{2}}{\sum_{i:\mathbf{x}_{i}\in
R_{r}}h_{i}}-\frac{(\sum_{i:\mathbf{x}_{i}\in
R_{root}}g_{i})^{2}}{\sum_{i:\mathbf{x}_{i}\in R_{root}}h_{i}}\right],$

where $g_{i}=\partial_{\hat{y}^{k-1}}\mathcal{L}(y_{i},0)$ and
$h_{i}=\partial^{2}_{\hat{y}^{k-1}}\mathcal{L}(y_{i},0)$. This is equivalent
to a gradient boosting iteration from a model
$f(\mathbf{x})^{k-1}=\hat{y}^{k-1}=0$ predicting zero and a learning rate of
one (Chen & Guestrin (2016)). For squared error loss, this tree will be
identical to a CART tree built using exact loss reduction $\mathcal{R}$.

The second-order approximation allows the use of the adaptive tree complexity
method by Lunde et al. (2020) which adaptively stops the splitting process of
a branch based on an information criterion, alleviating the need for
hyperparameter tuning to optimize the tree complexity. Additionally, it makes
it easy to incorporate additional objectives to the loss function as long as
they are twice differentiable. Note that as the information criterion is
rooted in the loss function, any additional objectives added to the loss are
factored in when determining the tree complexity.

###  4.1 Stable Loss

In our approach to improve stability, we introduce a stable loss function that
incorporates a regularization term. Similar to explicit regularization using a
penalty term in the loss function to encourage simpler solutions and address
the high-variance problem (i.e., overfitting), the Stable Loss (SL) adds a
penalty term to encourage small change in the predictions (i.e., stability).
SL uses the predictions of the current tree $f_{0}$ on the combined dataset
$\mathcal{D}_{1}=\\{\mathbf{x}^{(1)},y^{(1)}\\}$ to regularize the loss
function to prevent the predictions of the updated tree $f_{1}$ from deviating
too much from the predictions of $f_{0}$. In other word, when calculating loss
reduction $\tilde{\mathcal{R}}$, it replaces the loss $\mathcal{L}$ with the
stability loss

$\mathcal{L}_{\text{SL}}\left(y_{i},f_{1}(\mathbf{x}_{i});f_{0}(\mathbf{x}_{i})\right)=\\\
\underbrace{\mathcal{L}(y_{i},f_{1}(\mathbf{x}_{i})}_{\textit{loss}})+\gamma_{i}\underbrace{\mathcal{S}(f_{0}(\mathbf{x}_{i}),f_{1}(\mathbf{x}_{i})}_{\textit{instability}}),$
(1)

where $\gamma_{i}$ determines the strength of regularization (or belief in
$f_{0}$). The choice of $\gamma_{i}$ is key, and we explore two distinct
approaches and a combination of these to determine its value.

####  4.1.1 Constant Regularization

The simplest approach employs a uniform regularization factor $\alpha$, such
that $\gamma_{i}=\alpha$ for every data point. Here $\alpha$ is a
hyperparameter that is tuned to balance the trade-off between predictability
and stability.

####  4.1.2 Uncertainty-Weigthed Regularization (UWR)

Intuitively, to preserve predictability, we want to be more stable for
confident regions of $f_{0}$ compared to uncertain ones. The UWR dynamically
adjusts $\gamma_{i}$ based on each data point contribution to the model’s
uncertainty. Abusing the notations
$\hat{w}(\mathbf{x})=\hat{w}_{q(\mathbf{x})}$, we denote
$\hat{\sigma}^{2}_{\hat{w}_{0}(\mathbf{x})}$ as the variance estimate of the
predictions $\hat{w}_{0}(\mathbf{x})$ under the model $f_{0}$. Dividing the
precision $\frac{1}{\hat{\sigma}^{2}_{\hat{w}_{0}(\mathbf{x})}}$ by the number
of data points $n(\mathbf{x})$ contributing to this variance leads to the idea
of using

$\phi(\mathbf{x}_{i})=\frac{c}{n(\mathbf{x}_{i})\hat{\sigma}^{2}_{\hat{w}_{0}(\mathbf{x}_{i})}+\epsilon}.$

for scaling constant $c>0$ and small $\epsilon>0$ to ensure numerical
stability. We introduce $\beta$, a hyperparameter akin to $\alpha$, and set
$\gamma_{i}=\beta\phi(\mathbf{x}_{i})$. This is reminiscent of relative
weights according to a Gaussian prior on $\hat{w}_{0}(\mathbf{x}_{i})$. The
scaling constant $c$ is chosen as the average response variance across leaf
nodes, $c=n_{0}^{-1}\sum_{i=1}\hat{\sigma}^{2}_{y_{i}|\mathbf{x}_{i}}$.
Assuming the variance is $O(n(\mathbf{x})^{-1})$, this places $\beta$ and
$\alpha$ on a comparable scale.

The greedy splitting in CART complicates estimating
$\sigma_{\hat{w}_{0}(\mathbf{x}_{i})}^{2}$. However, the relationship
$C=\mathbb{E}[h_{i}]\sigma_{\hat{w}_{d}}^{2}$ between estimator variance and
loss-based information criteria, $C$ (see Takeuchi (1976) or Burnham &
Anderson (2004) for an English version), aids in this regard. Specifically,
Lunde et al. (2020, Appendix A) provides a conservative estimate for $C$
correcting for greedy split points. We use this relationship to construct the
estimator for the prediction variance in node $d$:

$\displaystyle\hat{\sigma}_{\hat{w}_{d}}^{2}=\frac{\sum_{i\in
I_{d}}(g_{i}+h_{i}\hat{w}_{d})^{2}}{(\sum_{i\in
I_{d}}h_{i})^{2}}\cdot\frac{1+\mathbb{E}[\max_{j}B_{j,root}]}{2}.$

Here $B_{j,root}=\max_{\pi}S_{j}(\pi)\geq 1$ represents the maximum of a
specific Cox-Ingersoll-Ross process $S_{j}(\pi)$ (Cox et al. (1985)), with
$\pi$ indicating the data fraction passed to node $d$ from its parent with
splitting variable $j$. The initial term is the Huber Sandwich Estimator
(Huber, 1967), suitable if the tree structure $q$ wasn’t learned from
$\mathcal{D}_{0}$. The latter adjustment term grows with the number of
possible splits.

####  4.1.3 Combined Regularization

We can also combine the constant regularization and UWR and define
$\gamma_{i}=\alpha+\beta\phi(\mathbf{x}_{i})$. This combination allows for a
base regularization level ($\alpha$) while still adapting to the uncertainty
connected to individual predictions ($\beta$). The hyperparameter pair
$(\alpha,\beta)$ dictates the regularization strategy, (0,0) corresponds to
the baseline model (i.e., minimizing the squared error with no
regularization), whereas ($\alpha>0$, 0) and (0, $\beta>0$) corresponds to
constant regularization and UWR strategy, respectively. The stable loss given
a hyperparameter pair $(\alpha,\beta)$ results in the following estimate of a
leaf prediction

$\hat{w}_{d}=-\frac{\sum_{i\in I_{d}}g_{i}}{\sum_{i\in I_{d}}h_{i}}.$

For squared error loss, this is equivalent to

$\hat{w}_{d}=\frac{\sum_{i\in
I_{d}}\left(y_{i}+\left(\alpha+\beta\phi(\mathbf{x}_{i})\right)\hat{w}_{0}(\mathbf{x}_{i})\right)}{\sum_{i\in
I_{d}}\left(1+\alpha+\beta\phi(\mathbf{x}_{i})\right)}.$

##  5 Experiments

In this section, we evaluate the Pareto efficiency of the proposed update
method with different $(\alpha,\beta)$ settings on different regression tasks.
Note that because we try to achieve two objectives (loss and instability
minimization) simultaneously, there is a trade-off between the two objectives.
The Pareto front represents all the solutions that have the best-achieved loss
for a given level of stability and vice versa. For our experiments, we use
squared error loss and squared error instability. We first describe the
experimental setups and then present the results.

###  5.1 Experimental Setup

Datasets. We evaluate our proposed method on various benchmark regression
datasets to get a comprehensive predictability-stability assessment across a
broad range of data characteristics, detailed in Table 1.

Dataset $\mathbf{n\times m}$ Description California $20,640\times 8$
Aggregated housing data in the California district using one row per census
block group. Boston $506\times 13$ Housing data for the suburbs of Boston.
Carseats $400\times 11$ Simulated data of child car seat sales. College
$777\times 18$ Statistics for a large number of US Colleges. Hitters
$263\times 20$ Major League Baseball Data from the 1986 and 1987 seasons. Wage
$3000\times 16$ Wage data for male workers in the Mid-Atlantic region.

Table 1: Datasets Description. The California dataset is adapted from the book
Sparse Spatial Autoregressions, Statistics and Probability Letters (Kelley
Pace & Barry (1997)), whereas the other datasets are from the book
”Introduction to Statistical Learning” (ISLR) by James et al. (2022). The
dimensions are after data preprocessing, i.e., one-hot encoding and removal of
any missing values.

We perform a comprehensive grid search across 121 configurations of the
hyperparameters $\alpha$ and $\beta$, where both range from 0 to 2 in
increments of 0.2. For all experiments we set $\epsilon=0.01$. Baseline.
Squared error loss, corresponds to hyperparameter setting ($\alpha=0$,
$\beta=0$). Constant. Corresponds to hyperparameter setting ($\alpha>0$,
$\beta=0$). UWR. Corresponds to hyperparameter setting ($\alpha=0$,
$\beta>0$); Combined. Corresponds to hyperparameter setting ($\alpha>0$,
$\beta>0$).

The setups. The setups aim to reflect how a model needs to adapt when new
information becomes available. The main setup consists of a single update
iteration $t=1$. First, a model is initialized and trained on the currently
available data $\mathcal{D}_{0}$. At $t=1$ new information is available that
is added to $\mathcal{D}_{0}$ to obtain $\mathcal{D}_{1}$ and the model is
updated using SL. The loss and stability are computed using repeated cross-
validation. For each dataset, the dataset is divided into five folds, where
each fold serves as test data once. The remaining four folds represent
$\mathcal{D}_{1}$ and half of $\mathcal{D}_{1}$ represents $\mathcal{D}_{0}$.
This process is repeated ten times and the average of the loss and stability
estimates from test data is used to obtain the final loss and stability
measures.

Figure 2: The loss and stability of the model update strategies, highlighting
the models that form the Pareto frontier with black borders.

Secondary setups. Two additional experiments are conducted using the
California dataset. We select five different hyperparameter configurations
from the Pareto-efficient configurations identified in the main setup, each
representing a different trade-off. The varied sample sizes setup uses 5000
data points as test data to compute loss and stability. The dataset
$\mathcal{D}_{1}$ is obtained by samples $n$ samples without replacement from
the remaining data. Half of $\mathcal{D}_{1}$ represents $\mathcal{D}_{0}$. An
initial model is trained on $\mathcal{D}_{0}$ and updated using
$\mathcal{D}_{1}$. We run the experiment for different values of
$n\in[1000,5000,10000,15000]$ and repeat it 50 times and average loss and
stability estimates on test data. The iterative updates setup divides the
dataset into seven folds, where one fold serves as $\mathcal{D}_{0}$ and one
serves as test data. The remaining five folds represent additional
information, $\mathcal{D}_{\Delta_{t}}$, obtained at different update
iterations. An initial model is trained on $\mathcal{D}_{0}$. For each update
iteration a fold is added to $\mathcal{D}_{0}$ to obtain $\mathcal{D}_{t}$ and
the model is updated and evaluated on the test data. The experiment is
repeated 50 times and averages loss and stability estimates on test data.

###  5.2 Experimental Results

Figure 2 shows the results of the main experiment for the California dataset.
The baseline is the least stable model and is not on the Pareto frontier. As
we introduce regularization, more predictive and/or stable solutions can be
achieved by tuning $\alpha$ and $\beta$. The constant approach shows a
consistent reduction in instability as $\alpha$ increases and up to a certain
$\alpha$ value, it also reduces loss. The UWR approach results in more stable
models than the baseline with lower loss compared to both the baseline and
constant approach. The combination of the constant and the UWR approach leads
to better solutions than the pure constant and UWR solutions, showcased by the
Pareto frontier.

Figure 3 shows and compares the results of the main experiment for all six
datasets. Notably, the baseline does not appear on the Pareto frontier for any
of the datasets, with all ($\alpha$, $\beta$) configurations being more
stable. Furthermore, many configurations also achieve lower loss than the
baseline. Comparing the constant and UWR approaches, the constant approach
tends to lead to more stable solutions, whereas the UWR approach tends to
achieve lower loss, although there is an exception for Wage datasets. The
combined approach frequently finds Pareto-efficient solutions that are
unattainable by either the constant or UWR approaches in isolation.

Figure 3: The figure plots the loss and stability of the model update
strategies and highlights the models that form the Pareto frontier with black
borders. The loss and stability values are relative to the baseline.

The result of the varied sample size experiment (Figure 4(a)) shows that loss
and stability decrease with larger sample sizes for baseline. For the five
selected ($\alpha$, $\beta$) configurations, loss decreases while instability
remains around the same level. Among the six models, all five configurations
are always on the Pareto frontier for all sample sizes. The configurations are
consistently more stable than the baseline at corresponding sample sizes. The
configurations (0.2, 0.6) and (0.4, 0.6) also consistently achieve lower
losses than the baseline.

Figure 4(b) presents the result of five selected ($\alpha$, $\beta$)
configurations over five update iterations. Compared to the baseline, each
configuration achieves improved stability at all update iterations. At the
first update iteration, the configurations (0.2, 0.6) and (0.4, 0.6) achieve a
lower loss than the baseline and at the second update iteration, also
configurations (1, 1) and (1.2, 0.4) achieve a lower loss. The next iteration,
(1, 1) becomes more stable than (1.2, 0.4) but with a higher loss than the
baseline. By the final update iteration, all configurations except (1, 1)
result in a more predictive solution than the baseline. At the first
iteration, all configurations are Pareto-efficient solutions, where the
configurations (0.2, 0.6), (0.4, 0.6), (1.2, 0.4) and (2, 0.2) are Pareto-
efficient models across all update iterations. Note that in practice one would
not stick to one configuration across all updates, but rather do a
hyperparameter search for each update iteration.

(a)  The varied sample sizes setup.

(b) The iterative updates setup.

Figure 4: The loss-instability trade-off of model updates with different
hyperparameter configurations $(\alpha,\beta)$ using the California Housing
dataset, highlighting Pareto efficient solutions with black borders: (a)
Effect of sample size on model updates. (b) Evolution of the loss-instability
trade-off over five update iterations, where each point indicates a model’s
performance at a specific iteration $t$.

##  6 Discussion

The objective of this paper was to provide stable updates for regression trees
that better balance the predictability-stability trade-off. We propose a
regularization method with two tunable hyperparameters. Through a series of
experiments conducted across various datasets and conditions, we demonstrated
the limitations of the baseline approach. Compared to our method, the baseline
consistently fails to be on the Pareto frontier, suggesting the necessity of
regularization to achieve a better trade-off. As we explore the regularization
landscape by adjusting $\alpha$ and $\beta$, we observe a diverse set of
Pareto efficient solutions, each representing a unique balance between
predictability and stability.

The main experiment on the California dataset revealed that the constant
approach tends to be more stability-focused as stability is equally important
for each data point. Incrementally increasing the constant regularization
factor $\alpha$ consistently improves model stability and up to a certain
value of $\alpha$ also improves predictability. Conversely, the UWR approach
tends to be more predictability-focused as stability for predictions in
uncertain regions of $f_{0}$ receives less weight. The fact that UWR often
improves predictability demonstrates the benefit of utilizing prediction
uncertainty for regularization. Combining the two approaches leads us to a
spectrum of Pareto-efficient solutions, where we can strategically choose
($\alpha$, $\beta$) configurations based on the desired predictability-
stability trade-off. These findings showcase that both regularization
approaches are relevant. Furthermore, these results also extend to the other
datasets, demonstrating the generalizability of our approach. Across all
datasets, the baseline is never one of the Pareto-efficient models. All
($\alpha$, $\beta$) configurations led to more stable models with many
configurations also achieving lower loss.

Our secondary experiments showed that our regularization method gives
consistent results regardless of sample size and that the benefits of our
method extend over multiple updates. When updating multiple times, our method
continues to outperform the baseline in terms of predictability and/or
stability. The models with lower $\alpha$ and $\beta$ values are the most
predictive with steady improvements in stability, while the models with large
$\alpha$ and small $\beta$ are the most stable and remain at the same
stability level across updates, with steady improvements in predictability.
Having large $\alpha$ and $\beta$ values seems to quickly become too
stability-focused, which comes at the cost of predictability.

A good choice of hyperparameters depends on the user’s preference.
Conservative users might prefer stability and therefore opt for larger
$\alpha$ and/or $\beta$ values, while users primarily concerned with
performance would rather opt for lower values. In practice, users should
specify a stability threshold before model selection.

A notable constraint of our approach is the requirement for batch training,
needing all available data to update the model. This is necessary as the CART
algorithm, on which our updating strategy relies, is a batch algorithm. Future
work could explore ways to remove the reliance on previous data to make the
method fit in an incremental learning setting.

The focus of our study was on regression tasks using squared error loss. It
would be interesting to apply stable updates of tree models to other tasks,
such as classification or survival analysis and use other types of loss and
instability functions, such as log loss and Poisson loss. This would also
expand the potential for industry applications.

While we show that our method improves the empirical stability of regression
trees in terms of semantic stability, we do not consider structural stability
in this paper. Structural stability is more directly related to explainability
and therefore of substantial interest. Investigating how our method influences
the structural stability of trees is a valuable direction for future research.

Our results show that iterative updates improve predictability. This leads us
to the following untested hypothesis: Learning a regression tree by starting
with a subset of the data and iteratively incorporating more data using our
update method improves upon the CART algorithm. This raises the questions of
how large the subset should be and how many updates to perform.

Lastly, our focus has been solely on individual regression trees. An
interesting extension would be to investigate how our method works in an
ensemble setting like random forests.

In conclusion, stable regression trees provide a good way of balancing the
predictability-stability trade-off, highlighted by its ability to outperform
the baseline approach across various datasets, data sizes and multiple update
iterations.

#### Acknowledgments

NB was supported in part by the Research Council of Norway grant
“Parameterized Complexity for Practical Computing (PCPC)” (NFR, no. 274526).

## References

  * Blørstad (2023) Morten Blørstad.  Improving Stability of Tree-Based Models.  Master’s thesis, University of Bergen, Bergen, 2023.
  * Bousquet & Elisseeff (2002) Olivier Bousquet and André Elisseeff.  Stability and Generalization.  _J. Mach. Learn. Res._ , 2:499–526, 3 2002.  ISSN 1532-4435.  doi: 10.1162/153244302760200704.  URL https://doi.org/10.1162/153244302760200704.
  * Breiman et al. (1984) L Breiman, Jerome H Friedman, Richard A Olshen, and C J Stone.  _Classification and Regression Trees_.  New York, 1st edition edition, 1984.  doi: https://doi.org/10.1201/9781315139470.
  * Breiman (1996) Leo Breiman.  Bagging predictors.  _Machine Learning_ , 24(2):123–140, 1996.  ISSN 1573-0565.  doi: 10.1007/BF00058655.  URL https://doi.org/10.1007/BF00058655.
  * Briand et al. (2009) Bénédicte Briand, Gilles R Ducharme, Vanessa Parache, and Catherine Mercat-Rommens.  A similarity measure to assess the stability of classification trees.  _Computational Statistics & Data Analysis_, 53(4):1208–1217, 2009.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2008.10.033.  URL https://www.sciencedirect.com/science/article/pii/S0167947308004970.
  * Burnham & Anderson (2004) Kenneth P. Burnham and David R. Anderson.  _Model Selection and Multimodel Inference_.  Springer New York, New York, NY, 2004.  ISBN 978-0-387-95364-9.  doi: 10.1007/b97636.
  * Chen & Guestrin (2016) Tianqi Chen and Carlos Guestrin.  XGBoost: A Scalable Tree Boosting System.  _arXiv_ , 3 2016.  doi: 10.1145/2939672.2939785.  URL https://arxiv.org/abs/1603.02754.
  * Cotter et al. (2018) Andrew Cotter, Heinrich Jiang, Serena Wang, Taman Narayan, Maya Gupta, Seungil You, and Karthik Sridharan.  Optimization with Non-Differentiable Constraints with Applications to Fairness, Recall, Churn, and Other Goals, 2018.
  * Cox et al. (1985) John C Cox, Jonathan E Ingersoll, and Stephen A Ross.  A Theory of the Term Structure of Interest Rates.  _Econometrica_ , 53(2):385–407, 1985.  ISSN 00129682, 14680262.  doi: 10.2307/1911242.  URL http://www.jstor.org/stable/1911242.
  * Dannegger (2000) Felix Dannegger.  Tree stability diagnostics and some remedies for instability.  _Statistics in medicine_ , 19:475–491, 1 2000.  doi: 10.1002/(SICI)1097-0258(20000229)19:43.0.CO;2-V.
  * Devroye & Wagner (1979) L Devroye and T Wagner.  Distribution-free performance bounds for potential function rules.  _IEEE Transactions on Information Theory_ , 25(5):601–604, 1979.  doi: 10.1109/TIT.1979.1056087.
  * Goh et al. (2016) Gabriel Goh, Andrew Cotter, Maya Gupta, and Michael P Friedlander.  Satisfying Real-world Goals with Dataset Constraints.  In D D Lee, M Sugiyama, U V Luxburg, I. Guyon, and R Garnett (eds.), _Advances in Neural Information Processing Systems 29_ , pp. 2415–2423. 2016.
  * Hastie et al. (2001) Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  _The Elements of Statistical Learning_.  Springer New York, New York, NY, 2001.  ISBN 978-0-387-84857-0.  doi: 10.1007/978-0-387-84858-7.
  * Huber (1967) Peter J Huber.  The behavior of maximum likelihood estimates under nonstandard conditions.  1967\.  URL https://api.semanticscholar.org/CorpusID:123642812.
  * James et al. (2022) Gareth James, Daniela Witten, Trevor Hastie, Rob Tibshirani, and Balasubramanian Narasimhan.  ISLR2.  Technical report, packaged 2022-11-20, 2022.
  * Kaminski & Malgieri (2021) Margot E Kaminski and Gianclaudio Malgieri.  Algorithmic impact assessments under the GDPR: producing multi-layered explanations.  _International Data Privacy Law_ , 11(2):125–144, 4 2021.  ISSN 2044-3994.  doi: 10.1093/idpl/ipaa020.  URL https://doi.org/10.1093/idpl/ipaa020.
  * Kearns & Ron (1999) Michael Kearns and Dana Ron.  Algorithmic Stability and Sanity-Check Bounds for Leave-One-Out Cross-Validation.  _Neural Computation_ , 11(6):1427–1453, 7 1999.  ISSN 0899-7667.  doi: 10.1162/089976699300016304.  URL https://doi.org/10.1162/089976699300016304.
  * Kelley Pace & Barry (1997) R Kelley Pace and Ronald Barry.  Sparse spatial autoregressions.  _Statistics & Probability Letters_, 33(3):291–297, 1997.  ISSN 0167-7152.  doi: https://doi.org/10.1016/S0167-7152(96)00140-X.  URL https://www.sciencedirect.com/science/article/pii/S016771529600140X.
  * Last et al. (2002) Mark Last, Oded Maimon, and Einat Minkov.  Improving Stability of Decision Trees.  _International Journal of Pattern Recognition and Artificial Intelligence_ , 16(02):145–159, 2002.  doi: 10.1142/S0218001402001599.  URL https://doi.org/10.1142/S0218001402001599.
  * Lim & Yu (2013) Chinghway Lim and B Yu.  Estimation Stability with Cross Validation (ESCV).  _Journal of Computational and Graphical Statistics_ , 25, 1 2013.  doi: 10.1080/10618600.2015.1020159.
  * Lin et al. (2002) Lawrence Lin, A. S Hedayat, Bikas Sinha, and Min Yang.  Statistical Methods in Assessing Agreement.  _Journal of the American Statistical Association_ , 97(457):257–270, 2002.  doi: 10.1198/016214502753479392.  URL https://doi.org/10.1198/016214502753479392.
  * Liu et al. (2022) Huiting Liu, Avinesh S., Siddharth Patwardhan, Peter Grasch, and Sachin Agarwal.  Model Stability with Continuous Data Updates, 2 2022.
  * Lunde et al. (2020) Berent Ånund Strømnes Lunde, Tore Selland Kleppe, and Hans Julius Skaug.  An information criterion for automatic gradient tree boosting.  _arXiv_ , 8 2020.  doi: 10.48550/arxiv.2008.05926.  URL https://arxiv.org/abs/2008.05926.
  * Milani Fard et al. (2016) Mahdi Milani Fard, Quentin Cormier, Kevin Canini, and Maya Gupta.  Launch and Iterate: Reducing Prediction Churn.  In D Lee, M Sugiyama, U Luxburg, I Guyon, and R Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 29. Curran Associates, Inc., 2016.  URL https://proceedings.neurips.cc/paper_files/paper/2016/file/dc5c768b5dc76a084531934b34601977-Paper.pdf.
  * Philipp et al. (2017) Michel Philipp, Thomas Rusch, Kurt Hornik, and Carolin Strobl.  Measuring the Stability of Results From Supervised Statistical Learning.  _Journal of Computational and Graphical Statistics_ , 131, 1 2017.  doi: 10.1080/10618600.2018.1473779.
  * Takeuchi (1976) K Takeuchi.  Distribution of information statistics and validity criteria of models.  _Mathematical Science_ , 153:12–18, 5 1976.
  * Turney (1995) Peter Turney.  Technical Note: Bias and the Quantification of Stability.  _Machine Learning_ , 20:23–33, 7 1995.  doi: 10.1007/BF00993473.


$N_{\mathrm{replay}}\to\infty$. With probability $\geq 1-3\delta$,

$\displaystyle J(\pi^{E})-J(\pi^{\texttt{RE}})\lesssim$ $\displaystyle\
\mathcal{L}_{1}+\mathcal{L}_{2}+\frac{\log\left(\mathcal{N}_{\mathrm{max}}H/\delta\right)+1}{N_{\mathrm{exp}}}$
(141)

where
$\mathcal{N}_{\mathrm{max}}\triangleq\max_{t\in[H]}\mathcal{N}(\mathcal{F}_{t},1/HN_{\mathrm{exp}},\|\cdot\|_{\infty})$
corresponds to the maximal covering number of the function classes
$\mathcal{F}_{t}$, and,

$\displaystyle\mathcal{L}_{1}\triangleq H^{2}\
\mathbb{E}_{\pi^{E}}\\!\left[\frac{\sum_{t=1}^{H}\mathcal{M}(s_{t},t)\mathsf{TV}\left(\pi^{E}_{t}(\cdot|s_{t}),\pi^{\texttt{BC}}_{t}(\cdot|s_{t})\right)}{H}\right],$
(142) $\displaystyle\mathcal{L}_{2}\triangleq
H^{3/2}\sqrt{\frac{\log\left(\mathcal{N}_{\mathrm{max}}H/\delta\right)}{N_{\mathrm{exp}}}\frac{\sum_{t=1}^{H}\mathbb{E}_{\pi^{E}}\left[1-\mathcal{M}(s_{t},t)\right]}{H}}.$

######  Remark 8.

Note that this line of reasoning can be extended to Theorem 3 and Theorem 8 to
show that the same guarantees as eq. 12 and eq. 131 respectively hold, but
with $F_{\mathrm{max}}$ replaced by $\mathcal{N}_{\mathrm{max}}$.

###  B.4 Bounds on RE in the linear expert setting

In this section, we provide an upper bound on the imitation gap of RE in the
presence of linear function approximation, which is studied in Rajaraman et
al. (2021a). This is a special case of the case of IL under parametric
function approximation with Lipschitzness. To avoid any ambiguity, we formally
define IL with linear function approximation, which is the case when $(i)$ the
expert follows an unknown linear classifier in a known set of features, and
$(ii)$ the reward function admits a linear parameterization.

The goal is to show that there exists a simple choice of the membership oracle
such that the imitation gap of the resulting algorithm grows as $H^{3/2}$ and
decay in the size of the expert dataset as $1/N_{\mathrm{exp}}$ up to
logarithmic factors, breaking the error compounding barrier and achieving the
optimal dependency on these parameters. We first introduce the linear setting
below.

######  Assumption 8 (Linear-expert setting).

For each $(s,a,t)$ tuple, the learner is assumed to have a feature
representation $\phi_{t}(s,a)\in\mathbb{R}^{d}$. For each time $t$, there
exists an unknown vector $\theta^{E}_{t}\in\mathbb{R}^{d}$ such that $\forall
s\in\mathcal{S}$,

$\displaystyle\pi^{E}_{t}(s)=\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\theta^{E}_{t},\phi_{t}(s,a)\rangle.$
(143)

i.e., the expert policy is deterministic and realized by a linear classifier.
We assume that $\theta_{t}^{E}\in\mathbb{S}^{d-1}$ without loss of generality.

######  Definition 7 (Policy induced by a linear classifier).

Consider a set of vectors $\theta=\\{\theta_{1},\cdots,\theta_{H}\\}$ where
each $\theta_{t}\in\mathbb{R}^{d}$. A policy $\pi^{\theta}$ is said to be
induced by the set of linear classifiers defined by $\theta$ if for all
$s\in\mathcal{S}$ and $t\in[H]$,

$\pi^{\theta}_{t}(s)=\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\theta_{t},\phi_{t}(s,a)\rangle.$
(144)

By this definition, $\pi^{E}=\pi^{\theta^{E}}$.

######  Definition 8 (Linear reward setting).

Define $\mathcal{R}_{\mathrm{lin,t}}$ as the family of linear reward functions
(defined at the single time-step $t$) which takes the form of an unknown
linear function of a set of the features,

$\displaystyle\mathcal{R}_{\mathrm{lin,t}}=\left\\{\\{r_{t}(s,a)=\langle\omega,\phi_{t}(s,a)\rangle:s\in\mathcal{S},a\in\mathcal{A}\\}:\omega\in\mathbb{R}^{d},\|\omega\|_{2}\leq
1\right\\}.$ (145)

For the rewards to be $1$-bounded, we assume the features satisfy
$\|\phi_{t}(s,a)\|_{2}\leq 1$. Define
$\mathcal{R}_{\mathrm{lin}}=\otimes_{t=1}^{H}\mathcal{R}_{\mathrm{lin,t}}$.
The linear reward setting assumes the true reward function of the MDP,
$r\in\mathcal{R}_{\mathrm{lin}}$.

######  Remark 9.

Note that our guarantees in Theorem 10 hold even if the set of features in the
definition of $\mathcal{R}_{\mathrm{lin,t}}$ in Definition 8 differ from those
used to define the expert classifier Assumption 8. Regardless, we assume that
both sets of features are known to the learner.

In the case of parametric function approximation with Lipschitzness, note that
we assume both the weak margin condition (Assumption 6), as well as the
existence of a linear classification oracle (Assumption 5). Below, in the
linear expert case, we show a sufficient condition which implies both of these
conditions. In particular, define the positive hemisphere with pole at
$\theta$, i.e. $\\{x:\mathbb{B}_{2}^{d}:\langle\theta,x\rangle\geq 0\\}$ as
$\mathbb{H}^{d}_{\theta}$. We abbreviate $\mathbb{H}^{d}_{\theta^{E}_{t}}$ as
$\mathbb{H}^{d}_{t}$.

######  Assumption 9 (Bounded density assumption).

For each time $t\in[H]$, state $s\in\mathcal{S}$, action $a\in\mathcal{A}$ and
$\theta\in\Theta_{t}$, define
$\overline{\phi}_{t}(s,a)=\phi_{t}(s,a_{s}^{\theta})-\phi_{t}(s,a)$ where
$a_{s}^{\theta}=\operatorname*{arg\,max}_{a^{\prime}\in\mathcal{A}}\langle\theta,\phi_{t}(s,a^{\prime})\rangle$.
Consider the measure $\mathrm{Pr}_{\pi^{E}}\left(\exists a\neq
a_{s_{t}}^{\theta}:\overline{\phi}_{t}(s_{t},a)\in\cdot\right)$. Let
$\overline{d}^{E}_{t}$ represent the Radon-Nikodym derivative of this measure
against the uniform measure on $\mathbb{H}^{d-1}_{t}$. The bounded density
assumption states that for each $t\in[H]$ there are constants $c_{\min}>0$ and
$c_{\max}<\infty$ such that for all $x\in\mathbb{H}^{d}_{t}$,

$\displaystyle c_{\min}\leq\overline{d}^{E}_{t}(x)\leq c_{\max}.$ (146)

We now state the main result we prove for IL in the linear setting.

######  Theorem 10.

Under Assumptions 8 and 9, appropriately instantiating RE ensures that with
probability $\geq 1-\delta$,

$J(\pi^{E})-J(\pi^{\texttt{RE}})\lesssim\sqrt{\frac{c_{\max}}{c_{\min}}}\frac{H^{3/2}d^{5/4}\log^{\frac{3}{2}}(N_{\mathrm{exp}}dH/\delta)}{N_{\mathrm{exp}}}.$

The proof of this result follows by showing that under Assumption 9, both the
weak margin condition (Assumption 6) is satisfied, and the classification
oracle (Assumption 5) can be constructed. We begin by showing the former.

######  Lemma 10.

Under Assumption 9, the $\alpha$-weak margin condition is satisfied with
$\alpha=1$ and $\mu=2c_{\max}\sqrt{d}$. In particular, for all
$\theta\in\mathbb{S}^{d-1}$,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(\langle\theta,\phi_{t}(s,a_{s_{t}}^{\theta})\rangle-\max_{a\neq
a_{s_{t}}^{\theta}}\langle\theta,\phi_{t}(s,a)\rangle\geq\eta\right)\geq
1-\left(2c_{\max}\sqrt{d}\right)\eta.$ (147)

where
$a_{s_{t}}^{\theta}\triangleq\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\theta,\phi_{t}(s,a)\rangle$.

###### Proof.

Observe that,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(\exists a\neq
a^{\theta}_{s_{t}}:\langle\theta,\phi_{t}(s_{t},a^{\theta}_{s_{t}})\rangle-\langle\theta,\phi_{t}(s_{t},a)\rangle\leq\eta\right)$
(148) $\displaystyle\quad=\mathrm{Pr}_{\pi^{E}}\left(\exists a\neq
a^{\theta}_{s_{t}}:\phi_{t}(s_{t},a^{\theta}_{s_{t}})-\phi_{t}(s_{t},a)\in\\{x\in\mathbb{H}_{\theta}^{d}:\langle
x,\theta\rangle\leq\eta\\}\right)$ (149)
$\displaystyle\quad\overset{(i)}{=}\mathrm{Pr}_{\pi^{E}}\left(\exists a\neq
a^{\theta}_{s_{t}}:\overline{\phi}_{t}(s_{t},a)\in\\{x\in\mathbb{H}_{\theta}^{d}:\langle
x,\theta\rangle\leq\eta\\}\right)$ (150)
$\displaystyle\quad\overset{(ii)}{\leq}c_{\max}\mathrm{Pr}(\langle
U,\theta\rangle\leq\eta)$ (151)

where in $(i)$, $\overline{\phi}_{t}$ is as defined in Assumption 9 and in
$(ii)$, $U$ is uniformly distributed on the unit hemisphere,
$\mathbb{H}^{d}_{\theta}$. Note that $(ii)$ follows from the bounded density
condition, Assumption 9. Note that the RHS essentially corresponds to the
volume (probability measure) of a disc of height $\eta$ cut out of a sphere
from the center. Up to normalization factors, this can be upper bounded by the
surface area of the base of the disc, multiplied by the height of the disc.
Namely,

$\displaystyle\frac{\eta\times\frac{\pi^{\frac{d-1}{2}}}{\Gamma\left(\frac{d-1}{2}+1\right)}}{\frac{1}{2}\frac{\pi^{\frac{d}{2}}}{\Gamma\left(\frac{d}{2}+1\right)}}$
(152)

Using Gautschi’s inequality, for any $x\geq 0$ and $\ell\in(0,1)$
$x^{1-\ell}\leq\frac{\Gamma(x+1)}{\Gamma(x+\ell)}\leq(1+x)^{1-\ell}$. With
$\ell=\frac{1}{2}$,
$\frac{\Gamma(\frac{d}{2}+1)}{\Gamma(\frac{d+1}{2})}\leq\sqrt{d}$. Combining
with eq. 151 results in,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(\exists a\neq
a^{\theta}_{s_{t}}:\langle\theta,\phi_{t}(s_{t},a^{\theta}_{s_{t}})\rangle-\langle\theta,\phi_{t}(s_{t},a)\rangle\leq\eta\right)\leq
2c_{\max}\sqrt{d}\eta$ (153)

Therefore the probability of the complement event is lower bounded by
$1-2c_{\max}\sqrt{d}\eta$, completing the proof. ∎

The final thing to show is that the bounded density assumption can also be
used to construct a classification oracle in the sense of Assumption 5.

Figure 8: If at any time $t\in[H]$ and state $s$,
$\phi_{t}(s,\pi^{E}_{t}(s))-\phi_{t}(s,a)$ lies in the red shaded region for
some action $a$, then, the action played by $\pi^{\texttt{BC}}$ and $\pi^{E}$
at this state are different.

In particular, as discussed in the main paper, we show that algorithms for
minimizing the generalization error, can be used to construct a classification
oracle. The compression based algorithm of Daniely & Shalev-Shwartz (2014)
provides a guarantee on the generalization error. From Theorem 5 of Daniely &
Shalev-Shwartz (2014), in the realizable setting, for linear classification,
the resulting classifier $\hat{\theta}$ has expected $0$-$1$ loss upper
bounded by $(d+\log(1/\delta))\log(n)/n$, given $n$ classification examples.
Namely, in the notation of Assumption 5, the resulting classifier
$\hat{\theta}$ satisfies with probability $\geq 1-\delta$,

$\displaystyle\mathrm{Pr}_{s\sim\mathcal{D}}\left(\operatorname*{arg\,max}_{a\in\mathcal{A}}f_{\theta^{*}}(s^{i},a)\neq\operatorname*{arg\,max}_{a\in\mathcal{A}}f_{\hat{\theta}}(s^{i},a)\right)\leq\frac{(d+\log(1/\delta)\log(n)}{n}.$
(154)

Next we show that under Assumption 9, this equation can be used to bound the
error in the parameter space, $\|\theta^{*}-\hat{\theta}\|_{2}$. Namely, in
Assumption 5, we may choose $\mathcal{E}_{\mathbb{B}_{2}^{d},n,\delta}$ as
$\asymp\frac{(d+\log(1/\delta)\log(n)}{n}$, up to constants depending on
$c_{\min}$.

######  Lemma 11.

Consider the compression based learner
$\widehat{\theta}^{\textsf{BC}}_{t}=\hat{\theta}_{t}$ of Daniely & Shalev-
Shwartz (2014) for multi-class linear classification. Then, under Assumption
9, with probability $\geq 1-\delta$,

$\displaystyle\|\widehat{\theta}^{\textsf{BC}}_{t}-\theta^{*}_{t}\|_{2}\leq\frac{2\pi}{c_{\min}}\frac{(d+\log(1/\delta)\log(N_{\mathrm{exp}})}{N_{\mathrm{exp}}}$
(155)

###### Proof.

Fix $t\in[H]$. The generalization error of
$\widehat{\theta}^{\textsf{BC}}_{t}=\hat{\theta}_{t}$ can be written as,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\theta^{*}_{t},\phi_{t}(s_{t},a)\rangle\neq\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\widehat{\theta}^{\textsf{BC}}_{t},\phi_{t}(s_{t},a)\rangle\right)$
$\displaystyle\quad=\mathrm{Pr}_{\pi^{E}}\left(\exists
a\neq\pi^{E}_{t}(s_{t}):\phi_{t}(s_{t},\pi^{E}_{t}(s_{t}))-\phi_{t}(s_{t},a)\in\mathcal{C}\right),$
(156)

where $\mathcal{C}$ is illustrated in fig. 8 and is formally defined as,

$\displaystyle\mathcal{C}\triangleq\\{x\in\mathbb{H}^{d}_{t}:\langle
x,\widehat{\theta}^{\textsf{BC}}_{t}\rangle\leq 0\\}.$ (157)

On the states which “belong” to $\mathcal{C}$ (i.e. at those states $s$ where
$\exists
a\neq\pi^{E}_{t}(s_{t}):\phi_{t}(s,\pi^{E}_{t}(s_{t}))-\phi_{t}(s,a)\in\mathcal{C}$),
there exists an action $a$ such that $\widehat{\theta}^{\textsf{BC}}_{t}$ is
better correlated with this action than $a^{*}_{s}$. In other words,
$\widehat{\theta}^{\textsf{BC}}_{t}$ and $\theta^{*}$ play different actions
at this state. Note that $\mathcal{C}$ is essentially the set difference of
two hemispheres with different poles. By the bounded density condition,
Assumption 9, and eq. 156,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(\pi^{E}_{t}(s_{t})\neq\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\widehat{\theta}^{\textsf{BC}}_{t},\phi(s,a)\rangle\right)\geq
c_{\min}\mathrm{Pr}\left(U\in\mathcal{C}\right),$ (158)

where $U$ is uniformly distributed over $\mathbb{H}^{d}_{\theta}$. Referring
to fig. 8, we have that,

$\displaystyle\mathrm{Pr}\left(U\in\mathcal{C}\right)=\frac{\alpha}{\pi}$
(159)

where $\alpha$ is the angle between $\widehat{\theta}^{\textsf{BC}}_{t}$ and
$\theta_{t}^{E}$. In particular, from eq. 158,

$\displaystyle\mathrm{Pr}_{\pi^{E}}\left(a^{*}_{s}\neq\operatorname*{arg\,max}_{a\in\mathcal{A}}\langle\widehat{\theta}^{\textsf{BC}}_{t},\phi(s,a)\rangle\right)\geq
c_{\min}\frac{\alpha}{\pi}\geq
c_{\min}\frac{\|\theta^{*}-\widehat{\theta}^{\textsf{BC}}_{t}\|_{2}}{\pi},$
(160)

where in the last inequality, we use the fact that
$\|\theta^{*}\|_{2}=\|\widehat{\theta}^{\textsf{BC}}_{t}\|_{2}=1$ without loss
of generality. By the generalization error bound on
$\widehat{\theta}^{\textsf{BC}}_{t}=\hat{\theta}_{t}$ in eq. 154, with
probability $\geq 1-\delta$,

$\displaystyle\|\theta^{*}-\widehat{\theta}^{\textsf{BC}}_{t}\|_{2}\leq\frac{\pi}{c_{\min}}\frac{(d+\log(1/\delta)\log(N_{\mathrm{exp}})}{N_{\mathrm{exp}}}$
(161)

∎

Lemma 11 shows that under the bounded density condition Assumption 9, the
compression based learner $\hat{\theta}$ of Daniely & Shalev-Shwartz (2014)
essentially induces a classification oracle for linear classification with
$\mathcal{E}_{\mathbb{B}_{2}^{d},n,\delta}=\frac{\pi}{c_{\min}}\frac{(d+\log(1/\delta)\log(n)}{n}$.
Finally, from Corollary 1, we have a bound on the covering number of linear
families. Putting together all of these results with Theorem 3 (noting that we
can replace $F_{\mathrm{max}}$ by $\mathcal{N}_{\mathrm{max}}$ from Remark 8)
results in Theorem 10.

##  Appendix C Additional experimental results

We include some additional experimental results in this section of the paper.

In fig. 9, we plot the distributions of the prefix weights generated by each
membership oracle on simulated BC rollouts on WalkerBulletEnv. Note that
$\mathcal{M}_{\texttt{VAR}}$ is significantly overconfident in prefix weights
compared to $\mathcal{M}_{\texttt{EXP}}$, as indicated by the heavier right-
tail. On the other hand, $\mathcal{M}_{\texttt{RND}}$ and
$\mathcal{M}_{\texttt{MAX}}$ are less overconfident and better overlap with
the idealized prefix weights induced by $\mathcal{M}_{\texttt{EXP}}$. This
aligns with the correlation plot between the various membership oracles in
Fig. 5. Moreover, in terms of policy performance, this further justifies the
superior behavior of $\mathcal{M}_{\texttt{MAX}}$ compared to
$\mathcal{M}_{\texttt{VAR}}$.

Figure 9: Histogram of prefix weights generated by rolling out trajectories
from BC . The green superimposed histogram represents prefix weights generated
by $\mathcal{M}_{\texttt{EXP}}$

In Fig. 10, we consider how each of the changes we described in the main
section of the paper lead to improved performance of our RE baseline. The
first, using a Wasserstein distance, leads to lower expected return but is
required for solving the full moment-matching problem – see Swamy et al.
(2021) for more details. Switching from PPO to the more sample-efficient SAC
(Haarnoja et al., 2018) leads to fast learning. Adding in gradient penalties
for discriminator stability (Swamy et al., 2021; Gulrajani et al., 2017) also
improves final performance and learning speed. The last change we employ,
using Optimistic Mirror Descent (Daskalakis et al., 2017) in both the
discriminator and RL algorithm also (slightly) improves performance. To our
knowledge, we are the first to utilize this technique in the imitation
learning literature and reccomend it as best practice for future moment-
matching algorithms. We refer interested readers to the work of Syrgkanis et
al. (2015) for theoretical details of why OMD enables superior performance.

Figure 10: We ablate the four key changes we made to off-the-shelf GAIL to
improve performance / theoretical guarantees. We see that each improved
performance, with MM significantly out-performing options with fewer changes.
Our improvements upon MM with the Replay Estimation technique are therefore
improving upon an already strong baseline.

##  Appendix D Experimental Setup

We begin with the hyperparameters for our Standard Bullet and Noisy Bullet
experiments.

###  D.1 Expert

We use the Stable Baselines 3 (Raffin et al., 2021) implementation of PPO
(Schulman et al., 2017) or SAC (Haarnoja et al., 2018) to train experts for
each environment. For the most part, we use already tuned hyperparameters from
(Raffin, 2020) in the implementation. The modifications we used are are shown
in table 1.

Parameter | Value
---|---
buffer size | 300000
batch size | 256
$\gamma$ | 0.98
$\tau$ | 0.02
Training Freq. | 64
Gradient Steps | 64
Learning Rate | Lin. Sched. 7.3e-4
policy architecture | 256 x 2
state-dependent exploration | true
training timesteps | 1e6
Table 1:  Expert hyperparameters for Walker Bullet Task and Hopper Bullet Task

####  D.1.1 Noisy Experts

In addition to the default Bullet Tasks, we test performance of algorithms on
noisy environments. Namely, we generate noisy expert data by re-training
expert policies with Gaussian noise added to the actions of the expert during
the exploration phase while training. We then re-generate expert data by
sampling from the expert policies trained on noisy data to analyze the
performance of our method under stochasticity. Table 2 lists the standard
deviation of the (i.i.d.) noise we applied to the actions in the different
environments.

env. | Noise Distribution.
---|---
hopper | $\mathcal{N}(0,0.1)$
walker | $\mathcal{N}(0,0.5)$
Table 2:  Noise we applied to all policies in each environment.

###  D.2 Baseline

We average over 5 runs and use a common architecture of 256 x 2 with ReLU
activations for both our method and the MM baseline we compare against. For
each datapoint, the cumulative reward is averaged over 10 trajectories. For
all tasks, we train on $\\{6,12,18\\}$ expert trajectories with a maximum of
400k iterations of the optimization procedure. Table 3 shows the
hyperparameters we used for MM. Empirically, smaller learning rates, large
batch sizes, and gradient penalties were critical for the stable convergence
of our method.

Parameter | Value
---|---
Batch Size | 2048*
Learning Rate | Linear Schedule of 8e-3*
$f$ Update Freq. | 5000
$f$ gradient target | 0.4
$f$ gradient penalty weight | 10
Table 3:  Learner hyperparameters for MM . * indicates the parameter was
different for the Hopper Initial State shift experiments (4096 for batch size
and Linear Schedule of 8e-4, respectively.).

We note that MM requires careful tuning of $f$ Update Freq. for strong
performance. We searched over step sizes of $\\{1250,2500,5000\\}$ and
selected the one which achieved the most stable updates. In practice, we
recommend evaluating a trained policy on a validation set to set this
parameter.

We also used similar parameters for training SAC, also from the Stable
Baselines 3 (Raffin et al., 2021) implementation, as we did for training the
expert policy. Table 4 shows the choice of hyperparmeters we used for training
SAC. We directly added in the optimistic mirror descent optimizers (Daskalakis
et al., 2017) for both the critic and actor objectives of SAC.

Parameter | Value
---|---
$\gamma$ | 0.98
$\tau$ | 0.02
Training Freq. | 64
Gradient Steps | 64
Learning Rate | Linear Schedule of 7.3e-4
policy architecture | 256 x 2
Table 4:  Leaning hyperparameters for the SAC component of MM

Table 5 shows the learning hyperparameters for any BC policies used for
generating simulated data for the membership oracles. Table 6 shows the number
of training steps per task we used for both the baseline and our method.

Parameter | Value
---|---
entropy weight | 0
l2 weight | 0
training timesteps | 1e5
Table 5:  Learner hyperparameters for Behavioral Cloning env. | training steps
---|---
walker (no noise) | 400000
walker (with noise) | 400000
hopper (no noise) | 400000
hopper (with noise) | 400000
Table 6:  Number of training steps for the different tasks

###  D.3 Our Algorithm

In this section, we use bold text to highlight sensitive hyperparameters. We
use the same network architecture choices as the MM baseline. For all
environments, we generated 100 trajectories of simulated behavior cloning data
to use with our method.

For all tasks, we rolled out 100 trajectories from a BC trained network to use
with our membership oracle. Table 7 shows how we partitioned our dataset
between the BC training set and the expert membership oracle dataset. We also
use the full dataset for moment matching, not just $D_{2}$, as we found this
lead to slightly better performance.

Expert Size | $D_{1}$ | $D_{2}$
---|---|---
6 trajs | 4 | 2
12 trajs | 10 | 2
18 trajs | 16 | 2
Table 7:  Partition of trajectories into $D_{1}$ and $D_{2}$ based on the
number of expert trajectories provided. For the Noisy Walker experiments, we
used $5,10,14$ trajectories for $D_{1}$ instead of the above.

###  D.4 Membership Oracle Parameters

For both $\mathcal{M}_{\texttt{VAR}}$ and $\mathcal{M}_{\texttt{MAX}}$, we use
5 BC networks in the ensemble. We followed the exact same parameters in Table
5 to train each BC imitator.

env | parameter | $\mathcal{M}_{EXP}$ | $\mathcal{M}_{\texttt{RND}}$ | $\mathcal{M}_{\texttt{VAR}}$ | $\mathcal{M}_{\texttt{MAX}}$
---|---|---|---|---|---
walker | $\beta$ | 0.1 | 0.1 | 0.01 | 0.1
walker | $\mu$ | 0.33 | 0.22 | 0.015 | 0.35
hopper | $\beta$ | 0.8 | 0.25 | 0.08 | 0.1
hopper | $\mu$ | 0.68 | 0.4 | 0.05 | 0.25
Table 8:  Membership oracle hyperparameters across different environments

Table 8 shows the choice of $\mu$ and $\beta$ values we used for each
membership oracle.

###  D.5 Initial State Shift Experiments

For these experiments, we used demonstrations generated by an expert trained
on the standard Bullet tasks but subject the learner (both at train and test
time) to a initial velocity perturbation of a zero-mean Gaussian with variance
($\sigma=1e-7$). We refer interested readers to our code for our precise
method of injecting noise as we believe it might be of interest for future
experiments. We note that in all the demonstrations, we see the expert start
from rest. Despite this relatively small shift, we see BC performance drop
significantly, as is characteristic of real-world problems where it
significantly under-performs on-policy IL methods. All results are averaged
over five seeds.

For all environments, we train BC for 1e5 steps (as well as for the query
policies for RE ).

For RE , we train $5$ policies and use the max-distance approximate membership
oracle. We use the above parameters for MM for our base moment-matcher.

env | parameter | $\mathcal{M}_{\texttt{MAX}}$
---|---|---
walker | $\beta$ | 0.01
walker | $\mu$ | 0.0001
hopper | $\beta$ | 0.01
hopper | $\mu$ | 0.0001
Table 9:  Membership oracle hyperparameters across different initial state
shift environments.


# SPLAIN: Augmenting Cybersecurity Warnings with Reasons and Data

Vera A. Kazakova  Jena D. Hwang  Bonnie J. Dorr  Yorick Wilks
<EMAIL_ADDRESS>

J. Blake Gage Alex Memory Mark A. Clark
<EMAIL_ADDRESS>

###### Abstract

Effective cyber threat recognition and prevention demand comprehensible
forecasting systems, as prior approaches commonly offer limited and,
ultimately, unconvincing information. We introduce Simplified Plaintext
Language for Actionable and Informative Narratives (SPLAIN), a natural
language generator that converts warning data into user-friendly cyber threat
explanations. SPLAIN is designed to generate clear, actionable outputs,
incorporating hierarchically organized explanatory details about input data
and system functionality. Given the inputs of individual sensor-induced
forecasting signals and an overall warning from a fusion module, SPLAIN
queries each signal for information on contributing sensors and data signals.
This collected data is processed into a coherent English explanation,
encompassing forecasting, sensing, and data elements for user review. SPLAIN’s
template-based approach ensures consistent warning structure and vocabulary.
SPLAIN’s hierarchical output structure allows each threat and its components
to be expanded to reveal underlying explanations on demand. Our conclusions
emphasize the need for designers to specify the “how” and “why” behind cyber
warnings, advocate for simple structured templates in generating consistent
explanations, and recognize that direct causal links in Machine Learning
approaches may not always be identifiable, requiring some explanations to
focus on general methodologies, such as model and training data.

##  1 Introduction: Spain NLG

Cyber-focused forecasting systems (e.g., Dalton et al. (2017); Sapienza et al.
(2017)) are becoming increasingly critical, with the growing threat and cost
of cyber crime Sobers (2019); Morgan (2019). Comprehensible, reliable, and
compelling forecasting systems are needed if cyber threats are to be
recognized and cyber-attacks avoided.

Warnings of imminent or existing cyber threats provide a basis and an
opportunity for counter-action. The likelihood of users taking effective
preventative actions depends on the system’s ability to accurately predict,
but also to clearly present and convincingly justify these predictions.
Existing systems, however, provide little information beyond threat levels and
confidence values.

We design Simplified Plaintext Language for Actionable and Informative
Narratives (SPLAIN), a natural language generator that converts complex
warning data into human-readable explanations. SPLAIN generates clear and
informative outputs, with actionable cyber-threat information, backed by
details about input data and system functionality.

SPLAIN’s inputs are individual sensor-induced forecasting signals coupled with
an overall warning produced by a fusion module. SPLAIN queries each input to
obtain information regarding contributing sensors and data signals. The
collected information is combined into a single coherent narrative in English
that conveys the overall warning—alongside forecasting, sensing, and data
elements—to be reviewed by the user at their desired level of detail. SPLAIN
employs a hierarchical template-based approach to produce warnings with
consistent structure and vocabulary. Each SPLAIN output can be further
expanded to reveal the underlying explanation for the overall warning.

We conclude that: (1) explainable outputs require that component designers
provide the specifics of how or why (not just what) behind cyber warnings; (2)
a simple set of structured templates can be developed for generating warnings
and explanations with consistent sentence structure and vocabulary; and (3)
direct causal links between inputs and outputs are not always identifiable
within Machine Learning approaches, necessitating that some explanations
describe general methodology only (e.g., model and training data).

##  2 Output: Explainable Warnings

SPLAIN is a component of a larger cyber-threat prediction system, dubbed
“Exploiting Leading Latent Indicators in Predictive Sensor Environments”
(ELLIPSE) Memory and Mueller (2019). This system processes large quantities of
web data in search of unconventional indicators of potential cyber threats,
i.e. those using NLP or data foraging, as opposed to hardware-based detectors.

SPLAIN is a natural language generator that combines different levels of
ELLIPSE outputs and transforms them into cyber threat predictions explained
through details regarding the system’s inputs and functionality. SPLAIN
adheres to the following tenets:

  * •

transparent system functionality

  * •

consistent & precise terminology

  * •

hierarchically expandable structure

  * •

justified assertions at all levels of explanations

##  3 System: Components and Data Flow

Figure 1 presents the flow of data through the ELLIPSE system, focusing on
elements relevant to SPLAIN.

Figure 1: Overall System’s Data Flow

Starting with data: files, tweets, websites feed into sensors, which count,
score, flag, or filter these inputs. The results are then passed to warning
generators (as well as other Sensors). These consume sensor data (standard
case), raw data, model data (e.g. GloVe in "MPEMEL"), and data funneled
through one-time collectors.

Warning Generator outputs are then passed into warning fusion, which combines
related incoming data into fused warnings. These fused outputs are passed into
SPLAIN, which uses ID tags to navigate from fused, to warning, to sensor, to
raw data, and then uses templates to generate and format the final warning
along with all of its supporting data.

##  4 NLG Approach: SPLAIN

SPLAIN employs a consistent format, using precise terminology and modular
structure: the design is easily modifiable if system components are altered,
expanded, or replaced. A hierarchical data organization is central to the
design, allowing for varying levels of explanation detail, on demand. Sensors
are subdivided by type (counter, scorer, event detector, and repository).

Figure 2 presents the components of the SPLAIN design. To produce warning
output, SPLAIN employs templates for scoring sensors. An example is provided
in Figure 3.

Figure 2: SPLAIN Architecture Figure 3: SPLAIN Template

##  5 Sample Output Warning

ELLIPSE incorporates a variety of unconventional threat indicators, including
an Outrage sensor Hollingshead et al. (2019). An example of SPLAIN-processed
output is the following generated text:

The outrage sensor calculated an average affect value of 53%, an average
intensity value of 48.4%, and an average outrage value of 70.91% toward target
X.

Additional structured representations are provided as output for this
explanation, including the context for the average values presented above:

Trigger #1: term “insanity” found.
Context tweet: "The new policies are pure insanity!" Trigger scores are:
affect value of 46%, intensity value of 55.8%, outrage value of 71.66% The
circumplex sentiment model was used.

Trigger #2: term “attack” found.
Context tweet: "This change is an attack on my wallet." Trigger scores are:
affect value of 60%, intensity value of 41%, outrage value of 70.15% The
circumplex sentiment model was used.

These triggers would be accompanied by sources and explanations of any
relevant models or approaches, such as the circumplex model Russell (1980).

##  6 Conclusions

SPLAIN employs a hierarchical template-based approach to produce expandable
warnings with consistent structure and vocabulary. Each SPLAIN output can be
further expanded to reveal the underlying explanation for the overall warning.
Explainable outputs require component designers to provide specifics of the
“HOW” and “WHY”, in addition to the “WHAT.” It is not always possible to
identify direct causal links between the inputs and outputs of Machine
Learning approaches; thus, some cyber-threat explanations must describe
general methodology only (e.g., model and training data).

## Acknowledgements

This research is supported by ODNI and IARPA via the AFRL contract number
FA8750-16-C-0114.

## References

  * Dalton et al. (2017) Adam Dalton, Bonnie Dorr, Leon Liang, and Kristy Hollingshead. 2017.  Improving cyber-attack predictions through information foraging.  In _Big Data (Big Data), 2017 IEEE International Conference on_ , pages 4642–4647. IEEE.
  * Hollingshead et al. (2019) Kristy Hollingshead, Bonnie J. Dorr, Adam Dalton, and Meg Barton. 2019.  Does outrage signal cyber attacks? predicting Bad Behavior from sentiment in online content.  In _Proceedings of FLAIRS-2019_.
  * Memory and Mueller (2019) Alex Memory and Graham Mueller. 2019.  Sensor fusion and structured prediction for cyberattack event networks.  In _Proceedings of the 15th International Workshop on Mining and Learning with Graphs (MLG)_.
  * Morgan (2019) Steve Morgan. 2019.  Cyber crime costs projected to reach $2 trillion by 2019.  _Forbes Magazine_.  Accessed: 2019-01-14.
  * Russell (1980) James A Russell. 1980.  A circumplex model of affect.  _Journal of personality and social psychology_ , 39(6):1161.
  * Sapienza et al. (2017) Anna Sapienza, Alessandro Bessi, Saranya Damodaran, Paulo Shakarian, Kristina Lerman, and Emilio Ferrara. 2017.  Early warnings of cyber threats in online discussions.  In _Data Mining Workshops (ICDMW), 2017 IEEE International Conference on_ , pages 667–674. IEEE.
  * Sobers (2019) Rob Sobers. 2019.  60 must-know cybersecurity statistics for 2019.  https://www.varonis.com/blog/cybersecurity-statistics/.  Accessed: 2019-01-14.