NRGBoost: Energy-Based Generative Boosted Trees

## Abstract

Abstract Despite the rise to dominance of deep learning in unstructured data domains, tree-based methods such as Random Forests (RF) and Gradient Boosted Decision Trees (GBDT) are still the workhorses for handling discriminative tasks on tabular data. We explore generative extensions of these popular algorithms with a focus on explicitly modeling the data density (up to a normalization constant), thus enabling other applications besides sampling.
As our main contribution we propose an energy-based generative boosting algorithm that is analogous to the second order boosting implemented in popular packages like XGBoost.
We show that, despite producing a generative model capable of handling inference tasks over any input variable,
our proposed algorithm can achieve similar discriminative performance to GBDT on a number of real world tabular datasets, outperforming alternative generative approaches.
At the same time, we show that it is also competitive with neural network based models for sampling.

## 1 Introduction

Generative models have achieved tremendous success in computer vision
and natural language processing,
where the ability to generate synthetic
data guided by user prompts opens up many exciting possibilities.
While generating synthetic table records does not necessarily enjoy the same wide appeal,
this problem has still received considerable attention
as a potential avenue for bypassing privacy concerns when sharing data.
Estimating the data density, $p(\mathbf{x})$, is another typical application of generative models which enables a host of different use cases that can be particularly interesting for tabular data.
Unlike discriminative models which are trained to perform inference over a single target variable,
density models can be used more flexibly for inference over different variables or for out of distribution detection.
They can also handle inference with missing data in a principled way by marginalizing over unobserved variables.

The development of generative models for tabular data has mirrored its progression in computer vision
with many of its Deep Learning (DL) approaches being adapted to the tabular domain
(Jordon et al., ; Xu et al., ; Fan et al., ; Engelmann & Lessmann, ; Zhao et al., ; Kotelnikov et al., ).
Unfortunately, these methods are only useful for sampling as they either don’t model the density explicitly or can’t evaluate it due to untractable marginalization over high dimensional latent variable spaces.
Furthermore, despite growing in popularity, DL has still failed to displace tree-based ensemble methods
as the tool of choice for handling tabular discriminative tasks
with gradient boosting still being
found to outperform neural-network-based methods in many real world datasets (Grinsztajn et al., ; Borisov et al., ).

While there have been recent efforts to extend the success of tree-based models to generative modeling (Correia et al., ; Wen & Hang, ; Nock & Guillame-Bert, ; Watson et al., ; Nock & Guillame-Bert, ; Jolicoeur-Martineau et al., ),
direct extensions of Random Forests (RF) and Gradient Boosted Decision Trees (GBDT) are still missing.
It is this gap that we try to address, seeking to keep the general algorithmic structure of these popular algorithms but replacing the optimization of their discriminative objective with a generative counterpart.
Our main contributions in this regard are:

- •
We propose NRGBoost, a novel energy-based generative boosting model that is trained to maximize a local second order approximation to the likelihood at each boosting round.
- •
We propose an amortized sampling approach that greatly reduces the training time of NRGBoost, which, like other energy-based models, is often dominated by sampling.
- •
We explore the use of bagged ensembles of Density Estimation Trees (DET) (Ram & Gray, ) with feature subsampling as a generative counterpart to Random Forests.

The longstanding popularity of GBDT models in machine learning practice can, in part, be attributed to the strength of its empirical results and the efficiency of its existing implementations.
We therefore focus on an experimental evaluation in real world datasets spanning a range of use cases, number of samples and features.
On smaller datasets, our implementation of NRGBoost
can be trained in a few minutes on a mid-range consumer CPU and achieves similar discriminative performance to a standard GBDT model, while remaining competitive with the best generative models for sampling.

Figure: Figure 1: Downsampled MNIST samples generated by the best generative models on this dataset. Despite this being a simple image dataset, it can be challenging for tabular generative models due to the high dimensionality and complex structure of correlations between features. We find NRGBoost to be the only model that is able to generate passable samples.
Refer to caption: x1.png

## 2 Energy Based Models

An Energy-Based Model (EBM) parametrizes the logarithm of a probability density function directly (up to an unspecified normalizing constant):

$$ $q_{f}(\mathbf{x})=\frac{\exp\left(f(\mathbf{x})\right)}{Z[f]}.$ (1) $$

Here $f(\mathbf{x}):\mathcal{X}\rightarrow\mathbb{R}$ is a real function over the input domain.(^1^11 We will assume that $\mathcal{X}$ is finite and discrete to simplify the notation and exposition but everything is applicable to bounded continuous input spaces, replacing the sums with integrals as appropriate.)
We will avoid introducing any parametrization,
instead treating the function $f\in\mathcal{F}(\mathcal{X})$ lying in an appropriate function space over the input space as our model parameter directly.
$Z[f]=\sum_{\mathbf{x}\in\mathcal{X}}\exp\left(f(\mathbf{x})\right)$, known as the partition function,
is then a functional of $f$ giving us the necessary normalizing constant.

This is the most flexible way one could represent a probability density function making essentially no compromises on its structure.
The downside to this is that for most interesting choices of $\mathcal{F}$, computing or estimating this normalizing constant is untractable which makes training these models difficult.
Their unnormalized nature, however, does not prevent EBMs from being useful in a number of applications besides sampling.
They can be used to perform inference over small enough groups of variables.
Partitioning the input space into a set of observed variables $\mathbf{x}_{o}$ and unobserved variables $\mathbf{x}_{u}$ we have

$$ $q_{f}(\mathbf{x}_{u}|\mathbf{x}_{o})=\frac{\exp\left(f(\mathbf{x}_{u},\mathbf{ x}_{o})\right)}{\sum_{\mathbf{x}^{\prime}_{u}}\exp\left(f(\mathbf{x}^{\prime}_ {u},\mathbf{x}_{o})\right)}\,,$ (2) $$

which only involves normalizing (i.e., computing a softmax) over the possible values of the unobserved variables.
Furthermore, for anomaly or out of distribution detection, knowledge of the normalizing constant is not necessary.

One common way to train an energy-based model to approximate a data generating distribution, $p(\mathbf{x})$, is to minimize the Kullback-Leibler divergence between $p$ and $q_{f}$, or equivalently,
maximize the expected log likelihood functional:

$$ $L[f]=\mathbb{E}_{\mathbf{x}\sim p}\log q_{f}(\mathbf{x})=\mathbb{E}_{\mathbf{x }\sim p}f(\mathbf{x})-\log Z[f]\,.$ (3) $$

This optimization is typically carried out by gradient descent over the parameters of $f$,
but due to the intractability of the partition function, one must rely on Markov Chain Monte Carlo (MCMC) sampling to estimate the gradients (Song & Kingma, ).

## 3 NRGBoost

Expanding the increase in log-likelihood in equation 3 due to a variation $\delta f$ around an energy function $f$ up to second order we have

$$ $L[f+\delta f]-L[f]\approx\mathbb{E}_{\mathbf{x}\sim p}\delta f(\mathbf{x})- \mathbb{E}_{\mathbf{x}\sim q_{f}}\delta f(\mathbf{x})-\frac{1}{2}\textrm{Var}_ {\mathbf{x}\sim q_{f}}\delta f(\mathbf{x})\eqqcolon\Delta L_{f}[\delta f]\,.$ (4) $$

The $\delta f$ that maximizes this quadratic approximation should thus have a large positive difference between the expected value under the data
and under $q_{f}$ while having low variance under $q_{f}$.
We note that just like the original log-likelihood, this Taylor expansion is invariant to adding an overall constant to $\delta f$.
This means that, in maximizing equation 4 we can consider only functions that have zero expectation under $q_{f}$ in which case we can simplify $\Delta L_{f}[\delta f]$ as

$$ $\Delta L_{f}[\delta f]=\mathbb{E}_{\mathbf{x}\sim p}\delta f(\mathbf{x})-\frac {1}{2}\mathbb{E}_{\mathbf{x}\sim q_{f}}\delta f^{2}(\mathbf{x})\,.$ (5) $$

We thus formulate our boosting algorithm as modelling the data density with an additive energy function. At each boosting iteration we improve upon the current energy function $f_{t}$ by finding an optimal step $\delta f^{*}_{t}$ that maximizes $\Delta L_{f_{t}}[\delta f]$:

$$ $\delta f^{*}_{t}=\arg\max_{\delta f\in\mathcal{H}_{t}}\Delta L_{f_{t}}[\delta f ]\,,$ (6) $$

where $\mathcal{H}_{t}$ is an appropriate space of functions (satisfying $\mathbb{E}_{\mathbf{x}\sim q_{f_{t}}}\delta f(\mathbf{x})=0$ if equation 5 is used).
The solution to this problem can be interpreted as a Newton step in the space of energy functions.
Because for an EBM, the Fisher Information matrix with respect to the energy function and the hessian of the expected log-likelihood are the same (see Appendix A),
we can also interpret the solution to equation 6 as a natural gradient step.
This approach is analogous to the second order step implemented in modern gradient boosting libraries such as XGBoost (Chen & Guestrin, )
and LightGBM (Ke et al., ) and which can be traced back to Friedman et al. ().

In updating the current iterate, $f_{t+1}=f_{t}+\alpha_{t}\cdot\delta f^{*}_{t}$, we scale $\delta f^{*}_{t}$ by an additional scalar step-size $\alpha_{t}$.
This can be interpreted as a globalization strategy to account for the fact that the quadratic approximation in equation 4
is not necessarily valid over large steps in function space. A common strategy in nonlinear optimization would be to select $\alpha_{t}$ via a line search based on the original log-likelihood.
Common practice in discriminative boosting however is to interpret this step size as a regularization parameter and to select a fixed value in $]0,1]$ with (more) smaller steps typically
outperforming fewer larger ones when it comes to generalization.
We choose to adopt a hybrid strategy, first selecting an optimal step size by line search and then shrinking it by a fixed factor.
We find that this typically accelerates convergence allowing the algorithm to take comparatively larger steps that increase the likelihood in the initial phase of boosting.

For a starting point, $f_{0}$, we can choose the logarithm of any probability distribution over $\mathcal{X}$ as long as it is easy to evaluate.
Sensible choices are a uniform distribution (i.e., $f\equiv 0$), the product of marginals for the training set, or any mixture distribution between these two.
In Figure 2 we show an example of NRGBoost starting from a uniform distribution and learning a toy 2D data density.

### 3.1 Weak Learners

Figure: Figure 2: Density learned by NRGBoost at different boosting iterations (1, 3, 10 and 100), for a toy dataset (shown right), starting from a uniform distribution. Weak learners are trees with 16 leaves.
Refer to caption: x2.png

As a weak learner we will consider piecewise constant functions defined by binary trees over the input space.
Letting $\bigcup_{j=1}^{J}X_{j}=\mathcal{X}$ be the partitioning of the input space induced by the leaves of a binary tree whose internal nodes represent a split along one dimension into two disjoint partitions,
we take as $\mathcal{H}$ the set of functions such as

$$ $\delta f(\mathbf{x})=\sum_{j=1}^{J}w_{j}\mathbf{1}_{X_{j}}(\mathbf{x})\,,$ (7) $$

where $\mathbf{1}_{X}$ denotes the indicator function of a subset $X$ and $w_{j}$ are values associated with each leaf $j\in[1..J]$.
In a standard decision tree these values would typically encode an estimate of $p(y|\mathbf{x}\in X_{j})$, with $y$ being a special *target* variable that is never considered for splitting.
In our generative approach they encode unconditional densities (or more precisely energies) over each leaf’s support and every variable can be used for splitting.
Note that our functions $\delta f$ are thus parametrized by the values $w_{j}$ as well the structure of the tree and the variables and values for the split at each node which ultimately determine the $X_{j}$.
We omit these dependencies for brevity.

Replacing the definition in equation 7 in our objective (equation 5) we get the following
optimization problem to find the optimal decision tree:

$$ $\displaystyle\max_{w_{1},\ldots,w_{J},X_{1},\ldots,X_{J}}$ $\displaystyle\sum_{j=1}^{J}\left(w_{j}P(X_{j})-\frac{1}{2}w_{j}^{2}Q_{f}(X_{j} )\right)$ (8) s.t. $\displaystyle\sum_{j=1}^{J}w_{j}Q_{f}(X_{j})=0\,,$ $$

where $P(X_{j})$ and $Q_{f}(X_{j})$ denote the probability of the event $\mathbf{x}\in X_{j}$ under the respective distribution and the constraint ensures that $\delta f$ has zero expectation under $q_{f}$.
With respect to the leaf weights this is a quadratic program whose optimal solution and objective values are respectively given by

$$ $\displaystyle w_{j}^{*}=\frac{P(X_{j})}{Q_{f}(X_{j})}-1\,,$ $\displaystyle\Delta L_{f}^{*}\left(X_{1},\ldots,X_{J}\right)=\frac{1}{2}\left( \sum_{j=1}^{J}\frac{P^{2}(X_{j})}{Q_{f}(X_{j})}-1\right)\,.$ (9) $$

Because carrying out the maximization of this optimal value over the tree structure that determines the $X_{j}$ is hard, we approximate its
solution by greedily growing a tree that maximizes it when considering how to split each node individually.
A parent leaf with support $X_{P}$ is thus split into two child leaves, with disjoint support, $X_{L}\cup X_{R}=X_{P}$, so as to maximize over all possible partitionings along a single dimension, $\mathcal{P}\left(X_{P}\right)$, the following objective:

$$ $\max_{X_{L},X_{R}\in\mathcal{P}\left(X_{P}\right)}\frac{P^{2}(X_{L})}{Q_{f}(X_ {L})}+\frac{P^{2}(X_{R})}{Q_{f}(X_{R})}-\frac{P^{2}(X_{P})}{Q_{f}(X_{P})}\,.$ (10) $$

Note that, when using parametric weak learners, computing a second order step would typically involve solving a linear system with a full Hessian.
As we can see, this is not the case when the weak learners are decision trees, where the optimal
value to assign to a leaf $j$ does not depend on any information from other leaves and, likewise, the optimal objective value is a sum
of terms, each depending only on information from a single leaf.
This would have not been the case had we tried to optimize the likelihood functional in Equation 3 directly instead of its quadratic approximation.

### 3.2 Amortized Sampling

To compute the leaf values in equation 9 and the splitting criterion in equation 10 we would have to know $P(X)$ and be able
to compute $Q_{f}(X)$ which is infeasible due to the untractable normalization constant.
We therefore estimate these quantities, with recourse to empirical data for $P(X)$, and to samples approximately drawn from the model with MCMC.
Because even if the input space is not partially discrete, $f$ is still discontinuous and constant almost everywhere
we can’t use gradient based samplers and therefore rely on Gibbs sampling instead.
This only requires evaluating each $f_{t}$ along one dimension at a time, while keeping all others fixed which can be computed efficiently for a tree by traversing it only once.
However, since at boosting iteration $t$ our energy function is a sum of $t$ trees, this computation scales linearly with the iteration number. This makes the overall time spent sampling quadratic in the number of boosting iterations and thus
precluding us from training models with a large number of trees.

In order to reduce the burden associated with this sampling, which dominates the runtime of training the model, we propose a new sampling approach that leverages the cumulative nature of boosting.
The intuition behind this approach is that the set of samples used in the previous boosting round are (approximately) drawn from a distribution that is already close to the new model distribution.
It could therefore be helpful to keep some of those samples, especially those that conform the best to the new model.
Rejection sampling allows us to do just that.
The boosting update in terms of the densities takes the following multiplicative form:

$$ $q_{t}(\mathbf{x})=k_{t}\,q_{t-1}(\mathbf{x})\exp\left(\alpha_{t}\delta f_{t}( \mathbf{x})\right)\,.$ (11) $$

Here, $k_{t}$ is an unknown multiplicative constant and since $\delta f_{t}$ is given by a tree, we can easily bound the exponential factor by finding the leaf with the largest value.
We can therefore use the previous model, $q_{t-1}(\mathbf{x})$, as a proposal distribution
for which we already have a set of samples and keep each sample, $\mathbf{x}$, with an acceptance probability of:

$$ $p_{accept}(\mathbf{x})=\exp\left[\alpha_{t}\left(\delta f_{t}(\mathbf{x})-\max _{\mathbf{x}}\delta f_{t}(\mathbf{x})\right)\right]\,.$ (12) $$

We note that knowledge of the constant $k_{t}$ is not necessary to compute this probability.

Our proposed sampling strategy is to maintain a fixed size pool of approximate samples from the model. At the end of each boosting round, we use rejection sampling to remove samples from the pool
and draw new samples from the model using Gibbs sampling.
Note that $q_{0}$ is typically a simple model for which we can both directly evaluate the desired quantities (i.e., $Q_{0}(X)$ for a given partition $X$) and cheaply draw exact samples from.
As such, no samples are required for the first iteration of boosting and for the second we can initialize the sample pool by drawing exact samples from $q_{1}$ with rejection sampling using $q_{0}$ as a proposal distribution.

This approach works better when the range of $\delta f_{t}$ and/or the step sizes $\alpha_{t}$ are small as this leads to larger acceptance probabilities.
Note that in practice it can be helpful to independently refresh a fixed fraction of samples, $p_{refresh}$, at each round of boosting in order to encourage more diverse samples between rounds.
This can be accomplished by keeping each sample with a probability $p_{accept}(\mathbf{x})(1-p_{refresh})$ instead.

### 3.3 Regularization

The simplest way to regularize a boosting model is to stop training when overfitting is detected by monitoring a suitable performance metric on a validation set.
For NRGBoost this could be the increase in log-likelihood at each boosting round.
However, estimating this quantity would require drawing additional validation samples from the model (see Appendix A.1).
An alternative viable validation strategy which needs no additional samples is to simply monitor a discriminative performance metric.
This amounts to monitoring the quality of $q_{f}(x_{i}|\mathbf{x}_{-i})$ instead of the full $q_{f}(\mathbf{x})$.

Besides early stopping, the decision trees themselves can be regularized by limiting the depth or total number of leaves in each tree.
Additionally, we can rely on other strategies
such as disregarding splits that would result in a leaf with too little training data, $P(X)$, model data, $Q_{f}(X)$, volume $V(X)$ or
too high of a ratio between training and model data $\nicefrac{{P(X)}}{{Q_{f}(X)}}$.
We found the latter to be the most effective of these, not only yielding better generalization performance than other approaches, but also
having the added benefit of allowing us to lower bound the acceptance probability of our rejection sampling scheme.
Furthermore, as we show in Appendix A.2, limiting this ratio guarantees that a small enough step size produces
an increase in likelihood at each round of boosting.

## 4 Density Estimation Trees and Density Estimation Forests

Density Estimation Trees (DET)
were proposed by Ram & Gray () as an alternative to histograms and kernel density estimation but have received little attention as generative models for sampling or other applications.
They model the density function as a constant value over the support of each leaf in a binary tree,

$$ $q(\mathbf{x})=\sum_{j=1}^{J}\frac{\hat{P}(X_{j})}{V(X_{j})}\mathbf{1}_{X_{j}}( \mathbf{x})\;,$ (13) $$

with
$\hat{P}(X)$ being an empirical estimate of probability of the event $\mathbf{x}\in X$ and $V(X)$ denoting the volume of $X$.
Note that it is possible to draw an exact sample from this type of model by randomly selecting a leaf, $j\in[1..J]$, given probabilities $\hat{P}(X_{j})$, and then drawing a sample from a uniform distribution over $X_{j}$.

To fit a DET, Ram & Gray () propose optimizing the Integrated Squared Error (ISE) between the data and model distributions which, following a similar approach to Section 3.1, leads the following optimization problem when considering how to split a leaf node:

$$ $\max_{X_{L},X_{R}\in\mathcal{P}\left(X_{P}\right)}D(P(X_{L}),V(X_{L}))+D(P(X_{ R}),V(X_{R}))-D(P(X_{P}),V(X_{P}))\,.$ (14) $$

For the ISE, $D$ should be taken as the function $D_{ISE}(P,V)=\nicefrac{{P^{2}}}{{V}}$
which leads to a similar splitting criterion to Equation 14 but replacing the previous model’s distribution with the volume measure $V$ which can be interpreted as the uniform distribution on $\mathcal{X}$ (up to a multiplicative constant).

##### Maximum Likelihood

Often generative models are trained to
maximize the likelihood of the observed data.
This was left for future work in Ram & Gray () but, as we show in Appendix B, can be accomplished by replacing $D_{ISE}$ with $D_{KL}(P,V)=P\,\log\left(\nicefrac{{P}}{{V}}\right)$.
This choice of minimization criterion can be seen as analogous to the choice between Gini impurity and Shannon entropy in the computation of the information gain in decision trees.

##### Bagging and Feature Subsampling

Following the common approach in decision trees, Ram & Gray () suggest the use of pruning for regularization of DET models.
Practice has however evolved to prefer bagging as a form of regularization rather than relying on single decision trees.
We employ same principle to DETs by fitting many trees on bootstrap samples of the data. We also adopt the common practice from Random Forests
of randomly sampling a subset of features to consider when splitting any leaf node in order to encourage independence between the different trees in the ensemble.
The ensemble model, which we call *Density Estimation Forests* (DEF) in the sequel, is thus an additive mixture of DETs with uniform weights, therefore still allowing for normalized density computation and exact sampling.

## 5 Related Work

##### Generative Boosting

Most prior work on generative boosting focuses on unstructured data and the use of parametric weak learners and is split between two approaches:
(i) Additive methods that model the density function as an additive mixture of weak learners such as Rosset & Segal (); Tolstikhin et al. ().
(ii) Those that take a multiplicative approach modeling the density function as an unnormalized product of weak learners.
The latter is equivalent to the energy based approach that writes the energy function (log density) as an additive sum of weak learners.
Welling et al. () in particular also approach boosting from the point of view of functional optimization of the likelihood or the logistic loss of an energy-based model.
However, they rely on a first order local approximation of the objective since they focus on parametric weak learners such as restricted Boltzmann machines for which a second order step would be impractical.

##### Greedy Multiplicative Boosting

Another more direct multiplicative boosting framework was first proposed by Tu ().
At each boosting round a discriminative classifier is trained to distinguish between empirical data and data generated by the current model by estimating
the likelihood ratio $\nicefrac{{p(\mathbf{x})}}{{q_{t}(\mathbf{x})}}$.
This estimated ratio is used as a direct multiplicative factor to update the current model $q_{t}$ (after being raised to an appropriate step size).
In ideal conditions
this greedy procedure would converge in a single iteration if a step size of 1 would be used.
While Tu () does not prescribe a particular choice of classifier,
Grover & Ermon () proposes a similar concept where the ratio is estimated based on an adversarial bound for an $f$-divergence and Cranko & Nock () provides additional analysis on this method.
We note that the main difference between this greedy approach and NRGBoost is that the latter attempts to update the current density proportionally to an exponential of the likelihood ratio, $\exp\left(\alpha_{t}\cdot\nicefrac{{p(x)}}{{q_{t}(x)}}\right)$, instead of $\left(\nicefrac{{p(\mathbf{x})}}{{q_{t}(\mathbf{x})}}\right)^{\alpha_{t}}$ directly.
In Appendix C we explore the differences between NRGBoost and this approach when it is adapted to use trees as weak learners.

##### Tree-Based Density Modelling

Other authors have proposed tree-based density models similar to DET (Nock & Guillame-Bert, ) or additive mixtures of tree-based models (Correia et al., ; Wen & Hang, ; Watson et al., ) but
perhaps surprisingly, the natural idea of creating an ensemble of DET models through bagging has not been explored before as far as we are aware.
Distinguishing features of some of these alternative approaches are:
(i) Not relying on a density estimation goal to drive the partitioning of the input space.
Correia et al. () leverages a standard discriminative Random Forest, therefore giving special treatment to a particular input variable whose conditional estimation drives the choice of partitions
and Wen & Hang () proposes using a mid-point random tree partitioning.
(ii) Leveraging more complex models for the data in the leaf of a tree instead of a uniform density (Correia et al., ; Watson et al., ).
This can allow for the use of trees that are more representative with a smaller number of leaves.
(iii) Relying on an adversarial framework where the generator and discriminator are both a tree (Nock & Guillame-Bert, ) or an ensemble of trees (Watson et al., ).
While this also leads to an iterative approach, unlike with boosting or bagging, the new model trained at each round doesn’t add to the previous one but replaces it instead.

##### Other Recent Tree-Based approaches

Nock & Guillame-Bert () proposes a different ensemble approach where each tree does not have its own leaf values that get added or multiplied to produce the final density,
but instead serve to collectively define the partitioning of the input space.
To train such models the authors propose a framework where, rather than adding a new tree to the ensemble at every iteration, the model is initialized with a fixed number of tree root nodes and each iteration adds a split to an existing leaf node.
Finally Jolicoeur-Martineau et al. () propose a diffusion model where a tree-based model (e.g., GBDT) is used to regress the score function. Being a diffusion model, however, means that computing densities is untractable.

## 6 Experiments

We evaluate the our proposed methods on 5 tabular datasets from the UCI Machine Learning Repository (Dheeru & Karra Taniskidou, ): Abalone (AB), Physicochemical Properties of Protein Tertiary Structure (PR), Adult (AD), MiniBooNE (MBNE) and Covertype (CT)
as well as the California Housing (CH) dataset available through scikit-learn (Pedregosa et al., ).
We also include a downsampled version of MNIST (by 2x along each dimension) which allows us to visually assess the quality of individual samples, something that is generally difficult with structured tabular data.
More details about these datasets are given in Appendix D.1.

We split our experiments into two sections, the first to evaluate the quality of density models directly on
a single variable inference task and the second to investigate the performance of our proposed models when used for sampling, comparing them against more specialized models.

### 6.1 Single Variable Inference

We test the ability of a generative model, trained to learn the density over all input variables, $q(\mathbf{x})$, to infer the value of a single one (i.e., we evaluate the quality of its estimate of $q(x_{i}|\mathbf{x}_{-i})$).
For this purpose we pick $x_{i}=y$ as the original target of the dataset, noting that the models that we train do not treat this variable in any special way, except for the selection of the best model in validation.
As such, we would expect that the model’s performance in inference over this particular variable is indicative of its strength on any other single variable inference task
and also indicative of the quality of the full $q(\mathbf{x})$ from which the conditional probability estimate is derived.

We use XGBoost (Chen & Guestrin, ) as a baseline for what should be achievable by a strong discriminative model,
noting that it is trained to maximize the discriminative likelihood, $\mathbb{E}_{\mathbf{x}\sim p}\log q(x_{i}|\mathbf{x}_{-i})$, directly,
not wasting model capacity in learning other aspects of the full data distribution.
We also compare to two other tree-based generative baselines: RFDE (Wen & Hang, ) and ARF (FORDE) (Watson et al., ).
The former allows us to gauge the impact of the guided partitioning used in DEF models over a random partitioning of the input space.

We use random search to tune the hyperparameters of the XGBoost model and a grid search to tune the most important hyperparameters of each generative density model.
We employ 5-fold cross-validation, repeating the hyperparameter tuning on each fold.
For the full details of the experimental protocol please refer to Appendix D.

**Table 1: Discriminative performance of different methods at inferring the value of a target variable. We use $R^{2}$ for regression tasks, AUC for binary classification and accuracy for multiclass classification. The reported values are means and standard errors over 5 cross-validation folds. The best *generative method* for each dataset is highlighted in bold and other methods that are not significantly worse (as determined by a paired $t$-test at a 95% confidence level) are underlined.**
|  | $R^{2}$ $\uparrow$ | AUC $\uparrow$ | Accuracy $\uparrow$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | AB | CH | PR | AD | MBNE | MNIST | CT |
| XGBoost | 0.552 ±0.035 | 0.849 ±0.009 | 0.678 ±0.004 | 0.927 ±0.000 | 0.987 ±0.000 | 0.976 ±0.002 | 0.971 ±0.001 |
| RFDE | 0.071 ±0.096 | 0.340 ±0.004 | 0.059 ±0.007 | 0.862 ±0.002 | 0.668 ±0.008 | 0.302 ±0.010 | 0.679 ±0.002 |
| ARF | 0.531 ±0.032 | 0.758 ±0.009 | 0.591 ±0.007 | 0.893 ±0.002 | 0.968 ±0.001 | -^∗ | 0.938 ±0.005 |
| DEF (ISE) | 0.467 ±0.037 | 0.737 ±0.008 | 0.566 ±0.002 | 0.854 ±0.003 | 0.653 ±0.011 | 0.206 ±0.011 | 0.790 ±0.003 |
| DEF (KL) | 0.482 ±0.027 | 0.801 ±0.008 | 0.639 ±0.004 | 0.892 ±0.001 | 0.939 ±0.001 | 0.487 ±0.007 | 0.852 ±0.002 |
| NRGBoost | 0.547 ±0.036 | 0.850 ±0.011 | 0.676 ±0.009 | 0.920 ±0.001 | 0.974 ±0.001 | 0.966 ±0.001 | 0.948 ±0.001 |

We find that NRGBoost outperforms the remaining generative models (see Table 1), even achieving comparable performance to XGBoost on the smaller datasets and with a small gap on the three larger ones (MBNE, MNIST, CT).
We note also that for the regression datasets, the generative models provide an estimate of the full conditional distribution over the target variable rather than a point estimate like XGBoost.
While there are other variants of boosting that can do the same (Duan et al., ), they rely on a parametric assumption about $p(y|\mathbf{x})$ that needs to hold for any $\mathbf{x}$.

### 6.2 Sampling

We evaluate two different aspects of the quality of generated samples: their utility for training a machine learning model and how distinguishable they are from real data.

Besides ARF, we compare to TVAE (Xu et al., ) and TabDDPM (Kotelnikov et al., ), two neural-network based generative models, as well as Forest-Flow (Jolicoeur-Martineau et al., ), a tree-based diffusion model.
Note that these three methods are not capable of density estimation.

##### Machine Learning Efficiency

Machine learning (ML) efficiency has been a popular way to measure the quality of generative models for sampling (Xu et al., ; Borisov et al., ; Kotelnikov et al., ).
It relies on using samples from a generative model to train a discriminative model which is then evaluated on real data,
and is thus similar to the single variable inference performance that we use to compare density models in Section 6.1.
In fact, if the density model’s support covers that of the full data, one would expect the discriminative model to recover the generator’s $q(y|\mathbf{x})$, and therefore its performance, in the limit where infinite generated data is used to train it.

We use an XGBoost model as the discriminative model and train it on the same number of training samples as the original data.
For the density models, we generate samples from the best model found in the previous section and for TVAE and TabDDM we select their hyperparameters by evaluating the ML Efficiency in the real validation set
(full details of the hyperparameter tuning are provided in Appendix D.2).
Note that this leaves these models at a potential advantage since the hyperparameter selection is based on the metric that is being evaluated.

We repeat all experiments 5 times, with different generated datatsets from each model and report the performance of the discriminative model in Table 2.
We find that NRGBoost and TabDDPM alternate as the best-performing model depending on the dataset (with two inconclusive cases), and NRGBoost is never ranked lower than second on any dataset.
We note also that, despite our best efforts with additional manual tuning, we could not achieve a reasonable model with TabDDM on MNIST.

**Table 2: Performance of an XGBoost model trained on synthetic data and on real data (for reference). For consistency, we use the same discriminative metrics as in Table 1 for evaluating its performance. Reported values are the averages and standard errors over 5 synthetic datasets generated by the same generative model. The best generative method for each dataset is highlighted in bold and methods that are not significantly worse (as determined by a $t$-test at a 95% confidence level) are underlined.**
|  | $R^{2}$ $\uparrow$ | AUC $\uparrow$ | Accuracy $\uparrow$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | AB | CH | PR | AD | MBNE | MNIST | CT |
| Real Data | 0.554 | 0.838 | 0.682 | 0.927 | 0.987 | 0.976 | 0.972 |
| TVAE | 0.483 ±0.006 | 0.758 ±0.005 | 0.365 ±0.005 | 0.898 ±0.001 | 0.975 ±0.000 | 0.770 ±0.009 | 0.750 ±0.002 |
| TabDDPM | 0.539 ±0.018 | 0.807 ±0.005 | 0.596 ±0.007 | 0.910 ±0.001 | 0.984 ±0.000 | - | 0.818 ±0.001 |
| Forest-Flow | 0.418 ±0.019 | 0.716 ±0.003 | 0.412 ±0.009 | 0.879 ±0.003 | 0.964 ±0.001 | 0.224 ±0.010 | 0.705 ±0.002 |
| ARF | 0.504 ±0.020 | 0.739 ±0.003 | 0.524 ±0.003 | 0.901 ±0.001 | 0.971 ±0.001 | 0.908 ±0.002 | 0.848 ±0.002 |
| DEF (KL) | 0.450 ±0.013 | 0.762 ±0.006 | 0.498 ±0.007 | 0.892 ±0.001 | 0.943 ±0.002 | 0.230 ±0.028 | 0.753 ±0.002 |
| NRGBoost | 0.528 ±0.016 | 0.801 ±0.001 | 0.573 ±0.008 | 0.914 ±0.001 | 0.977 ±0.001 | 0.959 ±0.001 | 0.895 ±0.001 |

Figure: Figure 3: Joint histogram for the latitude and longitude for the California Housing dataset.
Refer to caption: x3.png

##### Discriminator Measure

Similar to Borisov et al. () we test the capacity of a discriminative model to distinguish between real and generated data.
We use the original validation set as the real part of the training data for the discriminator in order to avoid benefiting generative methods that overfit their original training set.
A new validation set is carved out of the original test set (20%) and used to tune the hyperparameters of an XGBoost model which we use as our discriminator.
We report the AUC of this model on the remainder of the real test data in Table 3 which shows that NRGBoost outperforms other methods except on the PR (with an inconclusive result) and MBNE datasets.

**Table 3: AUC of an XGBoost model trained to distinguish real from generated data (lower means better). Reported values are the averages and standard errors over 5 synthetic datasets generated by the same model. The best generative method for each dataset is highlighted in bold and methods that are not significantly worse (as determined by a $t$-test at a 95% confidence level) are underlined.**
|  | AB | CH | PR | AD | MBNE | MNIST | CT |
| --- | --- | --- | --- | --- | --- | --- | --- |
| TVAE | 0.971 ±0.004 | 0.834 ±0.006 | 0.940 ±0.002 | 0.898 ±0.001 | 1.000 ±0.000 | 1.000 ±0.000 | 0.999 ±0.000 |
| TabDDPM | 0.818 ±0.015 | 0.667 ±0.005 | 0.628 ±0.004 | 0.604 ±0.002 | 0.789 ±0.002 | - | 0.915 ±0.007 |
| Forest-Flow | 0.987 ±0.002 | 0.926 ±0.002 | 0.885 ±0.002 | 0.932 ±0.002 | 1.000 ±0.000 | 1.000 ±0.000 | 0.985 ±0.001 |
| ARF | 0.975 ±0.005 | 0.973 ±0.004 | 0.795 ±0.008 | 0.992 ±0.000 | 0.998 ±0.000 | 1.000 ±0.000 | 0.989 ±0.001 |
| DEF (KL) | 0.823 ±0.013 | 0.751 ±0.008 | 0.877 ±0.002 | 0.956 ±0.002 | 1.000 ±0.000 | 1.000 ±0.000 | 0.999 ±0.000 |
| NRGBoost | 0.625 ±0.017 | 0.574 ±0.012 | 0.631 ±0.006 | 0.559 ±0.003 | 0.993 ±0.001 | 0.943 ±0.003 | 0.724 ±0.006 |

##### Statistical Distance Measure

As an additional measure of the statistical dissimilarity between data and model distributions we evaluate a Wasserstein distance using a similar setup to Jolicoeur-Martineau et al. ().
Namely, we min-max scale the numerical variables, and one-hot encode the categorical variables and scale them by 1/2.
We use a $L_{1}$ distance in the formulation of the optimal transport problem. Since finding the optimal solution has a worst case cubic scaling with the sample size,
we sub-sample a maximum of 5000 samples from the original train or test sets and use an equal number of synthetic samples.
We repeat the evaluation for each method and dataset 5 times using different synthetic data and also different subsampling seeds for the real data (where applicable).

**Table 4: Wasserstein distance between empirical test distribution and synthetic samples. Reported results are averages and standard errors over 5 repeated experiments. Smaller is better.**
|  | AB | CH | PR | AD | MBNE | MNIST | CT |
| --- | --- | --- | --- | --- | --- | --- | --- |
| TVAE | 0.302 ±0.005 | 0.206 ±0.010 | 0.274 ±0.008 | 1.086 ±0.023 | 0.425 ±0.032 | 20.61 ±0.042 | 0.807 ±0.023 |
| TabDDPM | 0.500 ±0.023 | 0.169 ±0.009 | 0.195 ±0.003 | 0.938 ±0.029 | 109.5 ±93.13 | - | 0.646 ±0.022 |
| Forest-Flow | 0.248 ±0.020 | 0.205 ±0.009 | 0.231 ±0.008 | 1.243 ±0.029 | 4.295 ±3.697 | 20.73 ±0.574 | 0.893 ±0.021 |
| ARF | 0.218 ±0.015 | 0.240 ±0.009 | 0.238 ±0.006 | 1.378 ±0.020 | 0.732 ±0.664 | 18.14 ±0.060 | 0.709 ±0.018 |
| DEF (KL) | 0.337 ±0.006 | 0.210 ±0.012 | 0.264 ±0.011 | 1.696 ±0.023 | 18.29 ±14.90 | 261.0 ±35.01 | 0.935 ±0.024 |
| NRGBoost | 0.502 ±0.045 | 0.178 ±0.011 | 0.225 ±0.010 | 1.051 ±0.013 | 22.33 ±20.11 | 15.70 ±0.188 | 0.602 ±0.015 |

**Table 5: Wasserstein distance between empirical train distribution and synthetic samples. Reported results are averages and standard errors over 5 repeated experiments. Smaller is better.**
|  | AB | CH | PR | AD | MBNE | MNIST | CT |
| --- | --- | --- | --- | --- | --- | --- | --- |
| TVAE | 0.299 ±0.005 | 0.195 ±0.010 | 0.274 ±0.011 | 1.075 ±0.024 | 0.366 ±0.055 | 20.59 ±0.041 | 0.808 ±0.029 |
| TabDDPM | 0.467 ±0.021 | 0.153 ±0.007 | 0.190 ±0.005 | 0.895 ±0.019 | 109.6 ±93.14 | - | 0.640 ±0.016 |
| Forest-Flow | 0.234 ±0.020 | 0.192 ±0.007 | 0.238 ±0.009 | 1.239 ±0.024 | 4.214 ±3.683 | 20.73 ±0.574 | 0.896 ±0.022 |
| ARF | 0.199 ±0.015 | 0.233 ±0.009 | 0.232 ±0.007 | 1.363 ±0.018 | 0.698 ±0.684 | 18.14 ±0.060 | 0.703 ±0.019 |
| DEF (KL) | 0.320 ±0.006 | 0.199 ±0.012 | 0.260 ±0.010 | 1.682 ±0.024 | 18.21 ±14.92 | 261.0 ±15.60 | 0.933 ±0.025 |
| NRGBoost | 0.509 ±0.047 | 0.170 ±0.012 | 0.230 ±0.011 | 1.028 ±0.015 | 22.25 ±20.09 | 15.60 ±0.212 | 0.600 ±0.016 |

Our results for the distance between empirical training and test distribution and synthetic samples from each model are reported in Tables 5 and 4 respectively.
We found these results to be somewhat sensitive to the choice of normalization of the numericals and categoricals, as well as to the choice of distance in this normalized space.
Despite NRGBoost not being as dominant as in the discriminator measure, it still fares well, obtaining the best average rank over all datasets and experiments.
Qualitatively speaking, we find that the samples generated by NRGBoost also look visually similar to the real data in both the MNIST and California datasets (see Figures 1 and 3), being the only generative model that can claim that for the former.

## 7 Discussion

The proposed DEF models are easy to sample from and require no sampling during training. However, we find that,
they require deep trees to model the data well which, in turn, also requires a larger number of trees in the ensemble to regularize.
In our experiments we found that the performance of DEF models was often capped by the maximum number of leaves we allowed ($2^{14}$).

In contrast, NRGBoost was able to model the data better while using shallower trees and in fewer number (the models selected on all datasets except for MNIST and CT had $\leq 2^{10}$ leaves).
Its main downside is that it can only be sampled from approximately using more expensive MCMC and also requires sampling during the training process.
While our fast Gibbs sampling implementation and our proposed sampling approach were able to mitigate the slow training, making these models practical,
they can still be cumbersome to use for sampling, particularly if independent samples are required, due to autocorrelation between samples from the same Markov Chain.
We argue however that unlike in image or text generation where fast sampling is necessary for an interactive user experience,
this can be less of a concern for the task of generating synthetic datasets where the one time cost of sampling is not as important as faithfully
capturing the data generating distribution.

A limitation of tree-based models when compared to DL approaches is that they require keeping the full training data in memory.
While discretization allows this data to be stored in a compact format, this could still limit the applicability of these methods to very large datasets.

## 8 Conclusion

We extend the two most popular tree-based discriminative methods for use in generative modeling.
We find that our boosting approach, in particular, offers generally good discriminative performance and competitive sampling performance to more specialized alternatives.
We hope that these results encourage further research into generative boosting approaches for tabular data,
in particular exploring other applications besides sampling that are enabled by (unnormalized) density models.

## Appendix A Theory

The expected log-likelihood for an energy-based model (EBM),

$$ $q_{f}(\mathbf{x})=\frac{\exp\left(f(\mathbf{x})\right)}{Z[f]}\,,$ (15) $$

is given by

$$ $L[f]=\mathbb{E}_{\mathbf{x}\sim p}\log q_{f}(\mathbf{x})=\mathbb{E}_{\mathbf{x }\sim p}f(\mathbf{x})-\log Z[f]\,.$ (16) $$

The *first variation* of $L$ can be computed as

$$ $\delta L[f;g]\coloneqq\left.\frac{dL[f+\epsilon g]}{d\epsilon}\right|_{ \epsilon=0}=\mathbb{E}_{\mathbf{x}\sim p}\,g(\mathbf{x})-\delta\log Z[f;g]= \mathbb{E}_{\mathbf{x}\sim p}\,g(\mathbf{x})-\mathbb{E}_{\mathbf{x}\sim q_{f}} \,g(\mathbf{x})\,.$ (17) $$

This is a linear functional of its second argument, $g$, and can be regarded as a directional derivative of $L$ at $f$ along a variation $g$.
The last equality comes from the following computation of the first variation of the log-partition function:

$$ $\displaystyle\delta\log Z[f;g]$ $\displaystyle=\frac{\delta Z[f;g]}{Z[f]}$ (18) $\displaystyle=\frac{1}{Z[f]}\sum_{\mathbf{x}}\exp^{\prime}\left(f(\mathbf{x}) \right)g(\mathbf{x})$ (19) $\displaystyle=\sum_{\mathbf{x}}\frac{\exp\left(f(\mathbf{x})\right)}{Z[f]}g( \mathbf{x})$ (20) $\displaystyle=\mathbb{E}_{\mathbf{x}\sim q_{f}}\,g(\mathbf{x})\,.$ (21) $$

Analogous to a Hessian, we can differentiate Equation 17 again along a second independent variation $h$ of $f$ yielding a symmetric bilinear functional which we will write as $\delta^{2}L[f;g,h]$.
Note that the first term in equation 16 is linear in $f$ and thus has no curvature, so we only have to consider the log partition function itself:

$$ $\displaystyle\delta^{2}L[f;g,h]$ $\displaystyle\coloneqq\left.\frac{\partial^{2}L[f+\epsilon g+\varepsilon h]}{ \partial\epsilon\partial\varepsilon}\right|_{(\epsilon,\varepsilon)=0}$ (22) $\displaystyle=-\delta^{2}\log Z[f;g,h]=-\delta\left\{\delta\log Z[f;g]\right\} [f;h]$ (23) $\displaystyle=-\delta\left\{\frac{1}{Z[f]}\sum_{\mathbf{x}}\exp\left(f(\mathbf {x})\right)g(\mathbf{x})\right\}[f;h]$ (24) $\displaystyle=\frac{\delta Z[f;h]}{Z^{2}[f]}\sum_{\mathbf{x}}\exp\left(f( \mathbf{x})\right)g(\mathbf{x})-\frac{1}{{Z[f]}}\sum_{\mathbf{x}}\exp^{\prime} \left(f(\mathbf{x})\right)g(\mathbf{x})h(\mathbf{x})$ (25) $\displaystyle=\frac{\delta Z[f;h]}{Z[f]}\cdot\mathbb{E}_{\mathbf{x}\sim q_{f}} g(\mathbf{x})-\frac{1}{{Z[f]}}\sum_{\mathbf{x}}\exp\left(f(\mathbf{x})\right)g (\mathbf{x})h(\mathbf{x})$ (26) $\displaystyle=\mathbb{E}_{\mathbf{x}\sim q_{f}}h(\mathbf{x})\cdot\mathbb{E}_{ \mathbf{x}\sim q_{f}}g(\mathbf{x})-\mathbb{E}_{\mathbf{x}\sim q_{f}}h(\mathbf{ x})g(\mathbf{x})$ (27) $\displaystyle=-\textrm{Cov}_{\mathbf{x}\sim q_{f}}\left(g(\mathbf{x}),h( \mathbf{x})\right)\,.$ (28) $$

Note that this functional is negative semi-definite for all $f$, i.e. $\delta^{2}L[f;h,h]\leq 0$, meaning that the log-likelihood is a concave functional of $f$.

Using these results, we can now compute the Taylor expansion of the increment in log-likelihood $L$ from a change $f\rightarrow f+\delta f$ up to second order in $\delta f$:

$$ $\displaystyle\Delta L_{f}[\delta f]$ $\displaystyle=\delta L[f;\delta f]+\frac{1}{2}\delta^{2}L[f;\delta f,\delta f]$ (29) $\displaystyle=\mathbb{E}_{\mathbf{x}\sim p}\delta f(\mathbf{x})-\mathbb{E}_{ \mathbf{x}\sim q_{f}}\delta f(\mathbf{x})-\frac{1}{2}\textrm{Var}_{\mathbf{x} \sim q_{f}}\delta f(\mathbf{x})\,.$ (30) $$

As an aside, defining the functional derivative, $\frac{\delta J[f]}{\delta f(\mathbf{x})}$, of a functional $J$ implicitly by:

$$ $\sum_{\mathbf{x}}\frac{\delta J[f]}{\delta f(\mathbf{x})}g(\mathbf{x})=\delta J [f;g]\,,$ (31) $$

we can formally define, by analogy with the parametric case, the Fisher Information "Matrix" (FIM) at $f$ as the following bilinear functional of two independent variations $g$ and $h$:

$$ $\displaystyle F[f;g,h]$ $\displaystyle\coloneqq-\sum_{\mathbf{y},\mathbf{z}}\left[\mathbb{E}_{\mathbf{x }\sim q_{f}}\frac{\delta^{2}\log q_{f}(\mathbf{x})}{\delta f(\mathbf{y})\delta f (\mathbf{z})}\right]g(\mathbf{y})h(\mathbf{z})$ (32) $\displaystyle=\sum_{\mathbf{y},\mathbf{z}}\frac{\delta^{2}\log Z[f]}{\delta f( \mathbf{y})\delta f(\mathbf{z})}g(\mathbf{y})h(\mathbf{z})$ (33) $\displaystyle=\delta^{2}\log Z[f;g,h]\,.$ (34) $$

The only difference to the second-order variation of 16 computed in equation 22 would be that the expectation is taken under the model distribution, $q_{f}$, instead of the data distribution $p$.
However, because the only term in $\log q_{f}(\mathbf{x})$ that is non-linear in $f$ is the log-partition functional, which is not a function of $\mathbf{x}$,
this expectation plays no role in the computation and we get the result that
the Fisher Information is the same as the negative Hessian of the log-likelihood for these models.

### A.1 Application to Piecewise Constant Functions

Considering a weak learner such as

$$ $\delta f(\mathbf{x})=\sum_{j=1}^{J}w_{j}\mathbf{1}_{X_{j}}(\mathbf{x})\,,$ (35) $$

where the subsets $X_{j}$ are disjoint and cover the entire input space, $\mathcal{X}$,
we have that

$$ $\displaystyle\mathbb{E}_{\mathbf{x}\sim q}\delta f(\mathbf{x})$ $\displaystyle=\sum_{\mathbf{x}\in\mathcal{X}}q(\mathbf{x})\sum_{j=1}^{J}w_{j} \mathbf{1}_{X_{j}}(\mathbf{x})$ (36) $\displaystyle=\sum_{j=1}^{J}w_{j}\sum_{\mathbf{x}\in X_{j}}q(\mathbf{x})\ = \sum_{j=1}^{J}w_{j}Q(X_{j})\,.$ (37) $$

Similarly, making use of the fact that $\mathbf{1}_{X_{i}}(\mathbf{x})\mathbf{1}_{X_{j}}(\mathbf{x})=\delta_{ij}
\mathbf{1}_{X_{i}}(\mathbf{x})$, we can compute

$$ $\displaystyle\mathbb{E}_{\mathbf{x}\sim q}\delta f^{2}(\mathbf{x})=\sum_{ \mathbf{x}\in\mathcal{X}}q(\mathbf{x})\left(\sum_{j=1}^{J}w_{j}\mathbf{1}_{X_{ j}}(\mathbf{x})\right)^{2}=\sum_{j=1}^{J}w_{j}^{2}Q(X_{j})\,.$ (38) $$

In fact, we can extend this to any ordinary function of $\delta f$:

$$ $\displaystyle\mathbb{E}_{\mathbf{x}\sim q}\,g\left(\delta f(\mathbf{x})\right)$ $\displaystyle=\sum_{\mathbf{x}\in\mathcal{X}}q(\mathbf{x})\sum_{j=1}^{J} \mathbf{1}_{X_{j}}(\mathbf{x})g\left(\delta f(\mathbf{x})\right)$ (39) $\displaystyle=\sum_{j=1}^{J}\sum_{\mathbf{x}\in X_{j}}q(\mathbf{x})g(w_{j})$ (40) $\displaystyle=\sum_{j=1}^{J}g(w_{j})Q(X_{j})\,,$ (41) $$

where we made use of the fact that the $\mathbf{1}_{X_{j}}$ constitute a partition of unity:

$$ $\displaystyle 1=\sum_{j=1}^{J}\mathbf{1}_{X_{j}}(\mathbf{x})\,.$ (42) $$

Finally, we can compute the increase in log-likelihood from a step $f\rightarrow f+\alpha\cdot\delta f$ as

$$ $\displaystyle L[f+\alpha\cdot\delta f]-L[f]$ $\displaystyle=\mathbb{E}_{\mathbf{x}\sim p}\left[\alpha\cdot\delta f(\mathbf{x })\right]-\log Z[f+\alpha\cdot\delta f]+\log Z[f]$ (43) $\displaystyle=\alpha\mathbb{E}_{\mathbf{x}\sim p}\delta f(\mathbf{x})-\log \mathbb{E}_{\mathbf{x}\sim q_{f}}\exp(\alpha\delta f(\mathbf{x}))$ (44) $\displaystyle=\alpha\sum_{j=1}^{J}w_{j}P\left(X_{j}\right)-\log\sum_{j=1}^{J}Q _{f}\left(X_{j}\right)\exp\left(\alpha w_{j}\right)\,,$ (45) $$

where in equation 44 we made use of the equality:

$$ $\displaystyle\log Z[f+\alpha\cdot\delta f]-\log Z[f]=\log\frac{\sum_{\mathbf{x }}\exp(f(\mathbf{x})+\alpha\delta f(\mathbf{x}))}{Z[f]}=\log\sum_{\mathbf{x}}q _{f}(\mathbf{x})\exp(\alpha\delta f(\mathbf{x}))\,,$ (46) $$

and of the result in equation 41 in the final step.

This result can be used to conduct a line search over the step size using training data and to estimate an increase in likelihood at each round of boosting for the purpose of early stopping, using validation data.

### A.2 A Simple Lower Bound on the Likelihood Increase

We can take the results from the last subsection to show that each boosting round improves the likelihood of the model under ideal conditions.
Substituting the NRGBoost leaf values $w_{j}=\nicefrac{{P(X_{j})}}{{Q_{f}(X_{j})}}-1$ in Equation 45 we get:

$$ $\displaystyle L[f+\alpha\cdot\delta f]-L[f]$ $\displaystyle=\alpha\left[\sum_{j=1}^{J}\frac{P^{2}\left(X_{j}\right)}{Q_{f}(X _{j})}-1\right]-\log\sum_{j=1}^{J}Q_{f}\left(X_{j}\right)\exp\left[\alpha\left (\frac{P\left(X_{j}\right)}{Q_{f}(X_{j})}-1\right)\right]$ (47) $\displaystyle=\alpha\chi^{2}(P\|Q_{f})-\log\mathbb{E}\exp\left(\alpha W\right)\,,$ (48) $$

where, with some abuse of notation, we denote by $\chi^{2}(P\|Q_{f})$ the $\chi^{2}$-divergence between
the two discrete distributions over $[1..J]$ induced by the partitioning obtained in the current boosting round.
This value is always non-negative, and can only be zero if $P\left(X_{j}\right)=Q_{f}(X_{j})$ for all leaves.
It therefore corresponds to an increase in likelihood that is as large as one can make the difference between the induced $P$ and $Q_{f}$ distributions by a judicious choice of $X_{j}$.
Note that this quantity is precisely the objective that NRGBoost tries to greedily maximize with its choice of $X_{j}$ (see Equation 9).

The second term in Equation 48 can be interpreted as the log of a moment generating function
for a centered random variable $W$ that takes values in $\left\{\nicefrac{{P\left(X_{j}\right)}}{{Q_{f}(X_{j})}}-1\right\}_{j\in[1..J]}$ (i.e., the set of the leaf values), with corresponding probabilities $\left\{Q_{f}(X_{j})\right\}_{j\in[1..J]}$.
With our proposed regularization approach (see Section 3.3), we limit the ratio $\nicefrac{{P\left(X_{j}\right)}}{{Q_{f}(X_{j})}}$ in any given $X_{j}$ to a maximum value $R$ when partitioning the domain
and therefore we have that $W\in[-1,R-1]$.

Since $W$ is bounded, by Hoeffding’s Lemma (1963),
it is sub-gaussian with variance proxy $R^{2}/4$ and the log of its moment generating function is thus upper bounded by:

$$ $\displaystyle\log\mathbb{E}\exp\left(\alpha W\right)\leq\frac{\alpha^{2}R^{2}} {8}\,.$ (49) $$

As a result, we have the following lower bound on the log-likelihood increase:

$$ $\displaystyle L[f+\alpha\cdot\delta f]-L[f]$ $\displaystyle\geq\alpha\left[\chi^{2}(P\|Q_{f})-\alpha\frac{R^{2}}{8}\right]\,.$ (50) $$

As long as in the current round of boosting a partitioning of the input space can be found that yields a non-zero $\chi^{2}(P\|Q_{f})$,
the likelihood is guaranteed to increase as long as we choose a small enough step-size: $\alpha\in\left(0,\frac{8}{R^{2}}\chi^{2}(P\|Q_{f})\right)$.
In particular, choosing a step size $\alpha=\frac{4}{R^{2}}\chi^{2}(P\|Q_{f})$ produces an increase in log-likelihood of at least $2\left[\nicefrac{{\chi^{2}(P\|Q_{f})}}{{R}}\right]^{2}$.

We note that this result, while insightful, assumes that one uses the exact probabilities under the model distribution in the NRGBoost update.
In practice, we would only be able to approximately estimate these with MCMC.
Furthermore, it can only guarantee an increase in the training log-likelihood as we also don’t have access to the probabilities under the data generating distribution.
We leave further analysis for future work.

## Appendix B Density Estimation Trees

Density Estimation Trees (DET) (Ram & Gray, ) model the density function as a piecewise constant function,

$$ $q(\mathbf{x})=\sum_{j=1}^{J}v_{j}\mathbf{1}_{X_{j}}(\mathbf{x})\,,$ (51) $$

where $X_{j}$ are given by a partitioning of the input space $\mathcal{X}$ induced by a binary tree
and the $v_{j}$ are the density values associated with each leaf that, for the time being, we will only require to be such that $q(\mathbf{x})$ sums to one.

Ram & Gray () proposes fitting DET models to directly minimize a generative objective, the Integrated Squared Error (ISE) between the data generating distribution, $p(\mathbf{x})$ and the model:

$$ $\min_{q\in\mathcal{Q}}\sum_{\mathbf{x}\in\mathcal{X}}\left(p(\mathbf{x})-q( \mathbf{x})\right)^{2}\,.$ (52) $$

Noting that $q$ is a function as in Equation 51 and that $\bigcup_{j=1}^{J}X_{j}=\mathcal{X}$, we can rewrite this as

$$ $\displaystyle\min_{v_{1},\ldots,v_{J},X_{1},\ldots,X_{J}}$ $\displaystyle\sum_{\mathbf{x}\in\mathcal{X}}p^{2}(\mathbf{x})+\sum_{j=1}^{J} \sum_{\mathbf{x}\in X_{j}}\left(v_{j}^{2}-2v_{j}p(\mathbf{x})\right)$ (53) s.t. $\displaystyle\sum_{j=1}^{J}\sum_{\mathbf{x}\in X_{j}}v_{j}=1\,.$ $$

Since the first term in the objective does not depend on the model this optimization problem can be further simplified as

$$ $\displaystyle\min_{v_{1},\ldots,v_{J},X_{1},\ldots,X_{J}}$ $\displaystyle\sum_{j=1}^{J}\left(v_{j}^{2}V(X_{j})-2v_{j}P(X_{j})\right)$ (54) s.t. $\displaystyle\sum_{j=1}^{J}v_{j}V(X_{j})=1\,,$ $$

where $V(X)$ denotes the volume of a subset $X$.
Solving this quadratic program for the $v_{j}$ we obtain the following optimal leaf values and objective:

$$ $\displaystyle v_{j}^{*}=\frac{P(X_{j})}{V(X_{j})}\,,$ $\displaystyle\text{ISE}^{*}\left(X_{1},\ldots,X_{J}\right)=-\sum_{j=1}^{J} \frac{P^{2}(X_{j})}{V_{f}(X_{j})}\,.$ (55) $$

One can therefore grow a tree by greedily choosing to split a parent leaf with support $X_{P}$ into two leaves with supports $X_{L}$ and $X_{R}$ so as to maximize the following criterion:

$$ $\max_{X_{L},X_{R}\in\mathcal{P}\left(X_{P}\right)}\frac{P^{2}(X_{L})}{V(X_{L}) }+\frac{P^{2}(X_{R})}{V(X_{R})}-\frac{P^{2}(X_{P})}{V(X_{P})}\,.$ (56) $$

### Maximum Likelihood

To maximize the likelihood,

$$ $\max_{q}\mathbb{E}_{\mathbf{x}\sim p}\log q(\mathbf{x})\,,$ (57) $$

rather than the ISE one can use the same approach. Here the optimization problem to solve is:

$$ $\displaystyle\max_{v_{1},\ldots,v_{J},X_{1},\ldots,X_{J}}$ $\displaystyle\sum_{j=1}^{J}P(X_{j})\log v_{j}$ (58) s.t. $\displaystyle\sum_{j=1}^{J}v_{j}V(X_{j})=1\,.$ $$

This is, again, easy to solve for $v_{j}$ since it is separable over $j$ after removing the constraint using Lagrange multipliers. The optimal leaf values and objective are in this case:

$$ $\displaystyle v_{j}^{*}=\frac{P(X_{j})}{V(X_{j})}\,,$ $\displaystyle L^{*}\left(X_{1},\ldots,X_{J}\right)=\sum_{j=1}^{J}P(X_{j})\log \frac{P(X_{j})}{V_{f}(X_{j})}\,.$ (59) $$

The only change is therefore to the splitting criterion which should become:

$$ $\max_{X_{L},X_{R}\in\mathcal{P}\left(X_{P}\right)}P(X_{L})\log\frac{P(X_{L})}{ V(X_{L})}+P(X_{R})\log\frac{P(X_{R})}{V(X_{R})}-P(X_{P})\log\frac{P(X_{P})}{V( X_{P})}\,.$ (60) $$

## Appendix C Greedy Tree-Based Multiplicative Boosting

In multiplicative generative boosting an unnormalized current density model, $\tilde{q}_{t-1}(\mathbf{x})$, is updated at each boosting round by multiplication with a new factor $\delta q_{t}^{\alpha_{t}}(\mathbf{x})$:

$$ $\tilde{q}_{t}(\mathbf{x})=\tilde{q}_{t-1}(\mathbf{x})\cdot\delta q_{t}^{\alpha _{t}}(\mathbf{x})\,.$ (61) $$

For our proposed NRGBoost, this factor is chosen in order to maximize a local quadratic approximation of the log likelihood around $q_{t-1}$ as a functional of the log density (see Section 3).
The motivation behind the greedy approach of Tu () or Grover & Ermon () is to instead make the update factor $\delta q_{t}(\mathbf{x})$ proportional to the likelihood ratio $r_{t}(\mathbf{x})\coloneqq\nicefrac{{p(\mathbf{x})}}{{q_{t-1}(\mathbf{x})}}$ directly,
which under ideal conditions would mean that the method converges immediately when choosing a step size $\alpha_{t}=1$.
In a more realistic setting, however, this method has been shown to converge under conditions on the performance of the individual $\delta q_{t}$ as discriminators between real and generated data (Tu, ; Grover & Ermon, ; Cranko & Nock, ).

While in principle this desired $r_{t}(\mathbf{x})$ could be derived from any binary classifier that is trained to predict a probability of a datapoint being generated (e.g., by training it to minimize a strictly proper loss)
and Tu () does not prescribe any particular choice, Grover & Ermon () propose relying on the following variational bound of an $f$-divergence to derive an estimator for this ratio:

$$ $D_{f}(P\|Q_{t-1})\geq\sup_{u\in\mathcal{U}_{t}}\left[\mathbb{E}_{\mathbf{x} \sim p}\;u(\mathbf{x})-\mathbb{E}_{\mathbf{x}\sim q_{t-1}}f^{*}(u(\mathbf{x})) \right]\,.$ (62) $$

Here $f^{*}$ denotes the convex conjugate of $f$. This bound is tight, with the optimum being achieved for $u_{t}^{*}(\mathbf{x})=f^{\prime}(\nicefrac{{p(\mathbf{x})}}{{q_{t-1}(\mathbf{x
})}})$,
if $\mathcal{U}_{t}$ is capable of representing this function.
$\left(f^{\prime}\right)^{-1}(u^{*}_{t}(\mathbf{x}))$ can thus be interpreted as an approximation of the desired $r_{t}(\mathbf{x})$.

Adapting this method to use trees as weak learners can be accomplished by considering $\mathcal{U}_{t}$ in Equation 62 to be defined by tree functions
$u=\nicefrac{{1}}{{J}}\sum_{j=1}^{J}w_{j}\mathbf{1}_{X_{j}}$ with leaf values $w_{j}$ and leaf supports $X_{j}$.
At each boosting iteration a new tree, $u^{*}_{t}$ can thus be grown to greedily optimize the lower bound in the r.h.s. of Equation 62 and setting
$\delta q_{t}(\mathbf{x})=\left(f^{\prime}\right)^{-1}(u^{*}_{t}(\mathbf{x}))$ which is thus also a tree with the same leaf supports and leaf values given by $v_{j}\coloneqq\left(f^{\prime}\right)^{-1}(w_{j})$.
This leads to the seaprable optimization problem:

$$ $\max_{w_{1},\ldots,w_{J},X_{1},\ldots,X_{J}}\sum_{j}^{J}\left[P(X_{j})w_{j}-Q( X_{j})f^{*}(w_{j})\right]\,.$ (63) $$

Note that we drop the iteration indices from this point onwards for brevity. Maximizing over $w_{j}$ with the $X_{j}$ fixed we have that $w_{j}^{*}=f^{\prime}\left(\nicefrac{{P(X_{j})}}{{Q(X_{j})}}\right)$ which yields the optimal value

$$ $J^{*}(X_{1},\ldots,X_{j})=\sum_{j}\left[P(X_{j})f^{\prime}\left(\frac{P(X_{j}) }{Q(X_{j})}\right)-Q(X_{j})(f^{*}\circ f^{\prime})\left(\frac{P(X_{j})}{Q(X_{j })}\right)\right]$ (64) $$

In turn, this determines the splitting criterion as a function of the choice of $f$. Finally, the optimal density values for the leaves are given by

$$ $v^{*}_{j}=\left(f^{\prime}\right)^{-1}(w^{*}_{j})=\frac{P(X_{j})}{Q(X_{j})}\,.$ (65) $$

It is interesting to note two particular choices of $f$-divergences.
For the KL divergence, $f(t)=t\log t$ and $f^{\prime}(t)=1+\log t=\left(f^{*}\right)^{-1}(t)$.
This leads to

$$ $J_{KL}(X_{1},\ldots,X_{j})=\sum_{j}P(X_{j})\log\frac{P(X_{j})}{Q(X_{j})}$ (66) $$

as the splitting criterion.
The Pearson $\chi^{2}$ divergence, with $f(t)=(t-1)^{2}$, leads to the same splitting criterion as NRGBoost.
Note however that for NRGBoost the leaf values for the multiplicative update of the density are proportional to $\exp\left(\nicefrac{{P(X_{j})}}{{Q(X_{j})}}\right)$ instead of the ratio directly.
Table 6 summarizes these results.

**Table 6: Comparison of splitting criterion and leaf weights for the different versions of boosting.**
|  | Splitting Criterion | Leaf Values (Density) |
| --- | --- | --- |
| Greedy (KL) | $P\log\left(\nicefrac{{P}}{{Q}}\right)$ | $\nicefrac{{P}}{{Q}}$ |
| Greedy ($\chi^{2}$) | $\nicefrac{{P^{2}}}{{Q}}$ | $\nicefrac{{P}}{{Q}}$ |
| NRGBoost | $\nicefrac{{P^{2}}}{{Q}}$ | $\exp\left(\nicefrac{{P}}{{Q}}-1\right)$ |

Another interesting observation is that a DET model can be interpreted as a single round of greedy multiplicative boosting
starting from a uniform initial model. The choice of the ISE as the criterion to optimize the DET corresponds to the choice of Pearson’s $\chi^{2}$ divergence and likelihood to the choice of KL divergence.

## Appendix D Reproducibility

### D.1 Datasets

We use 5 datasets from the UCI Machine Learning Repository (Dheeru & Karra Taniskidou, ): Abalone, Physicochemical Properties of Protein Tertiary Structure (from hereon referred to as Protein), Adult, MiniBooNE and Covertype.
We also use the California Housing dataset which was downloaded through the Scikit-Learn package Pedregosa et al. () and a downsampled version of the MNIST dataset Deng ().
Table 7 summarizes the main details of these datasets as well as the approximate number of samples used for train/validation/test for each cross-validation fold.

**Table 7: Dataset Information. We respect the original test sets of each dataset when provided, otherwise we set aside 20% of the original dataset as a test set. 20% of the remaining data is set aside as a validation set used for hyperparameter tuning.**
| Abbr | Name | Train + Val | Test | Num | Cat | Target | Cardinality |
| --- | --- | --- | --- | --- | --- | --- | --- |
| AB | Abalone | 3342 | 835 | 7 | 1 | Num | 29 |
| CH | California Housing | 16512 | 4128 | 8 | 0 | Num | Continuous |
| PR | Protein | 36584 | 9146 | 9 | 0 | Num | Continuous |
| AD | Adult | 32560 | 16280^∗ | 6 | 8 | Cat | 2 |
| MBNE | MiniBooNE | 104051 | 26013 | 50 | 0 | Cat | 2 |
| MNIST | MNIST (downsampled) | 60000 | 10000^∗ | 196 | 0 | Cat | 10 |
| CT | Covertype | 464810 | 116202 | 10 | 2 | Cat | 7 |

### D.2 DEF and NRGBoost Implementation Details

##### Discretization

In our practical implementation of tree based methods we first discretize the input space by binning continuous numerical variables by quantiles.
Furthermore we also bin discrete numerical variables in order to keep their cardinalities smaller than 256.
This can also be interpreted as establishing a priori a set of discrete values to consider when splitting on each numerical variable
and is done for computational efficiency, being inspired by LightGBM (Ke et al., ).

##### Categorical Splitting

For splitting on a categorical variable we once again take inspiration from LightGBM.
Rather than relying on one-vs-all splits we found it better to first order the possible categorical values at a leaf according to a pre-defined sorting function and then choose the optimal many-vs-many split as if the variable was numerical.
The function used to sort the values is the leaf value function.
E.g., for splitting on a categorical variable $x_{i}$ we order each possible categorical value $k$ by $\nicefrac{{\hat{P}(x_{i}=k,X_{-i})}}{{\hat{Q}(x_{i}=k,X_{-i})}}$ where $X_{-i}$ denotes the leaf support over the remaining variables.

##### Tree Growth Strategy

We always grow trees in best first order, i.e., we always split the current leaf node that yields the maximum gain in the chosen objective value.

##### Line Search

As mentioned in Section 3, we perform a line search to find the optimal step size after each round of boosting in order to maximize the likelihood gain in Equation 45.
Because evaluating multiple possible step sizes, $\alpha_{t}$, is inexpensive, we simply do a grid search over 101 different step sizes, split evenly in log space over $[10^{-3},10]$.

##### Code

Our implementation of the proposed tree-based methods is mostly Python code using the NumPy library (Harris et al., ). We implement the tree evaluation and Gibbs sampling in C, making use of the PCG library (O’Neill, ) for random number generation.
Our code is available in the supplementary material and we will make both the DEF and NRGBoost methods available in an open-source library.

### D.3 Hyperparameter Tuning

#### D.3.1 XGBoost

To tune the hyperparameters of XGBoost we use 100 trials of random search with the search space defined in Table 8.

**Table 8: XGBoost hyperparameter tuning search space. $\delta(0)$ denotes a point mass distribution at 0.**
| Parameter | Distribution or Value |
| --- | --- |
| learning_rate | $\textrm{LogUniform}\left(\left[10^{-3},1.0\right]\right)$ |
| max_leaves | $\textrm{Uniform}\left(\left\{16,32,64,128,256,512,1024\right\}\right)$ |
| min_child_weight | $\textrm{LogUniform}\left(\left[10^{-1},10^{3}\right]\right)$ |
| reg_lambda | $0.5\cdot\delta(0)+0.5\cdot\textrm{LogUniform}\left(\left[10^{-3},10\right]\right)$ |
| reg_alpha | $0.5\cdot\delta(0)+0.5\cdot\textrm{LogUniform}\left(\left[10^{-3},10\right]\right)$ |
| max_leaves | 0 (we already limit the number of leaves) |
| grow_policy | lossguide |
| tree_method | hist |

Each model was trained for 1000 boosting rounds on regression and binary classification tasks.
For multi-class classification tasks a maximum number of 200 rounds of boosting was used due to the larger size of the datasets and because a separate tree is built at every round for each class.
The best model was selected based on the validation set, together with the boosting round where the best performance was attained. The test metrics reported correspond to the performance of the selected model at that boosting round on the test set.

#### D.3.2 RFDE

We implement the RFDE method (Wen & Hang, ) after quantile discretization of the dataset and therefore split at the midpoint of the discretized dimension instead of the original one.
When a leaf support has odd cardinality over the splitting dimension a random choice is made over the two possible splitting values. Finally, the original paper does not mention how to split over categorical domains.
We therefore choose to randomly split the possible categorical values for a leaf evenly as we found that this yielded slightly better results than a random one vs all split.

For RFDE models we train a total of 1000 trees. The only hyperparameter that we tune is the maximum number of leaves per tree for which we test the values $[2^{6},2^{7},\ldots,2^{14}]$.
For the Adult dataset, due to limitations of our tree evaluation implementation we only test values up to $2^{13}$.

#### D.3.3 ARF

For ARF we used the official python implementation of the algorithm (Blesch & Wright, ).
We implement the single variable inference for this method ourselves since the python library does not provide density evaluation.
Due to memory (and time) concerns we train models with a maximum of 100 trees for up to 10 adversarial iterations.
We use a similar grid-search setup to DEF, tuning only the following two parameters:

- •
The maximum number of leaves per tree, for which we test the values $[2^{6},2^{7},\ldots,2^{12},\infty]$ ($\infty$ denotes unconstrained).
- •
The minimum number of examples falling on each leaf, for which we test the values $[3,5,10,30,50,100]$.
If one value of this parameter doesn’t improve upon the previous, we move on to the next value of the maximum number of leaves (i.e., we exit this loop early).

#### D.3.4 DEF

We train ensembles with 1000 DET models. Only three hyperparameters are tuned, using three nested loops, each loop running over the possible values of a single parameter in a pre-defined order. These are, in order of outermost to innermost:

- •
The maximum number of leaves per tree for which we test the values $[16384,4096,1024,256]$.
If the best model for one value of the maximum number of leaves (over the remaining inner-loop parameters) doesn’t improve upon the previous, we stop the hyperparameter tuning.
This is done because in most datasets, reducing the number of leaves always leads to a worse model for DEF.
- •
The fraction of features to consider when splitting a node. We test the values $[d^{-1/2},d^{-1/4},1]$ with $d$ being the dimension of the dataset.
If the best model for one value of this parameter doesn’t improve upon the best model for the previous value, we move on to the next value of the maximum number of leaves (i.e., we exit this loop early).
- •
The minimum number of data points that need to be left in each leaf for a split to be considered. We test the values $[0,1,3,10,30]$.
If one value of this parameter doesn’t improve upon the previous, we move on to the next value of fraction of features.

#### D.3.5 NRGBoost

We train NRGBoost models for a maximum of 200 rounds of boosting. We only tune two parameters for NRGBoost Models:

- •
The maximum number of leaves for which we try the values $[64,256,1024,4096]$ in order, stopping if performance fails to improve from one value to the next. For the CT dataset we also include 16384 in the values to test.
- •
The constant factor by which the optimal step size determined by the line search is shrunk at each round of boosting.
This is essentially the "learning rate" parameter.
To tune it we perform a Golden-section search for the log of its value using a total of 6 evaluations. The range we use is $[0.01,0.5]$.

For regularization we limit the ratio between empirical data density and model data density on each leaf to a maximum of 2.
We noticed that smaller values of this regularization parameter tend to work better but that it otherwise plays a similar role to the shrinkage factor so we opt to tune only the latter.

For sampling we use a sample pool of 80000 samples for all datasets except for Covertype where we use 320000.
These values are chosen so that the number of samples in the pool are at minimum similar to the training set size which we found to be a good rule of thumb.
At each round of boosting, samples are removed from the pool according to our rejection sampling approach (using Equation 12) after independently removing 10% at random.
These samples are replaced by samples from the current model using Gibbs sampling.

The starting point of each NRGBoost model was selected as a mixture model between a uniform distribution (10%) and the
product of training marginals (90%) on the discretized input space. We observed that this mixture coefficient does not have much impact on the results however.

#### D.3.6 TVAE

We use the implementation of TVAE from the SDV package.(^2^22[https://github.com/sdv-dev/SDV](https://github.com/sdv-dev/SDV))
To tune its hyperparameters we use 50 trials of random search with the search spaces defined in Table 9.

**Table 9: TVAE hyperparameter tuning search space. We set both compress_dims and decompress_dims to have the number of layers specified by num_layers, with hidden_dim hidden units in each layer. We use larger batch sizes and smaller number of epochs for the larger datasets (MBNE, MNIST, CT) since these can take signifcantly longer to run a single epoch.**
| Parameter | Datasets | Distribution or Value |
| --- | --- | --- |
| epochs | small | $\textrm{Uniform}\left(\left[100..500\right]\right)$ |
|  | large | $\textrm{Uniform}\left(\left[50..200\right]\right)$ |
| batch_size | small | $\textrm{Uniform}\left(\left\{100,200,\ldots,500\right\}\right)$ |
|  | large | $\textrm{Uniform}\left(\left\{500,1000,\ldots,2500\right\}\right)$ |
| embedding_dim | all | $\textrm{Uniform}\left(\left\{32,64,128,256,512\right\}\right)$ |
| hidden_dim | all | $\textrm{Uniform}\left(\left\{32,64,128,256,512\right\}\right)$ |
| num_layers | all | $\textrm{Uniform}\left(\left\{1,2,3\right\}\right)$ |
| compress_dims | all | (hidden_dim,) * num_layers |
| decompress_dims | all | (hidden_dim,) * num_layers |

#### D.3.7 TabDDPM

We used the official implementation of TabDDPM(^3^33[https://github.com/yandex-research/tab-ddpm](https://github.com/yandex-research/tab-ddpm)) adapted to use our datasets and validation setup.
To tune the hyperparameters of TabDDPM we use 50 trials of random search with the same search space that the original authors report in their paper (Kotelnikov et al., ).

#### D.3.8 Forest-Flow

We use the official implementation of Forest-Flow.(^4^44[https://github.com/SamsungSAILMontreal/ForestDiffusion](https://github.com/SamsungSAILMontreal/ForestDiffusion))
The only hyperparameters with impact on performance that the Forest-Flow documentation recommends tuning are the number of diffusion timesteps for which a recommended default value of 50 is offered and the number of noise values per sample for which the default value is 100.
Due to the fact that Forest-Flow models already take significantly longer to train than the other models, and that that increasing either value would lead to slower training we opted to use these default values.

### D.4 Evaluation Setup

##### Single variable inference

For the single variable inference evaluation, the best models are selected by their discriminative performance on a validation set.
The entire setup is repeated five times with different cross-validation folds and with different seeds for all sources of randomness.
For the Adult and MNIST datasets the test set is fixed but training and validation splits are still rotated.

##### Sampling

For the sampling evaluation we use a single train/validation/test split of the real data (corresponding to the first fold in the previous setup) for training the generative models.
The density models used are those previously selected based on their single variable inference performance on the validation set.
For TVAE and TabDDPM we directly evaluate their ML Efficiency using the validation data.

##### ML Efficiency

For each selected model we sample a train and validation sets with the same number of samples as those used on the original data.
For NRGBoost we generate these samples by running 64 chains in parallel with 100 steps of burn in and downsampling their outputs by 30 (for the smaller datasets) or 10 (for MBNE, MNIST and CT).
For every synthetic dataset, an XGBoost model is trained using the best hyperparameters found on the real data and using a synthetic validation set to select the best stopping round for XGBoost.
The setup is repeated 5 times with different datasets being generated for each method and dataset.

##### Discriminator Measure

We create the training, validation and test sets to train an XGBoost model to discriminate between real and generated data using the following process:

- •
The original validation set is used as the real part of the training set in order to avoid benefiting generative methods that overfit their training set.
- •
The original test set is split 20%/80%. The 20% portion is used as the real part of the validation set and the 80% portion as the real part of the test set.
- •
To form the synthetic part of the training, validation and test sets for the smaller datasets (AB, CH, PR, AD) we sample data according to the original number of samples in the train, validation and test splits on the real data respectively.
Note that this makes the ratio of real to synthetic data 1:4 in the training set. This is deliberate because for these smaller datasets the original validation has few samples and adding extra synthetic data helps the discriminator.
- •
For the larger datasets (MBNE, MNIST, CT) we generate the same number of synthetic samples as there are real samples on each split, therefore making every ratio 1:1 because the discriminator is typically already too powerful without adding additional synthetic data.

Because, in contrast to the previous metric, having a lower number of effective samples helps rather than hurts we take extra precautions to not generate correlated data with NRGBoost.
We draw each sample by running its own independent chain for 100 steps starting from an independent sample from the initial model which is a rather slow process.
The setup is repeated 5 times with 5 different sets of generated samples from each method.

##### Statistical Distance Measure

We follow a procedure similar to Jolicoeur-Martineau et al. ().
We min-max scale the numerical variables, and one-hot encode the categorical variables and scale them by 1/2.
We use a $L_{1}$ distance in the formulation of the optimal transport problem and use the python POT library to solve it. We sub-sample a maximum of 5000 samples from the original train or test sets and use an equal number of synthetic samples, repeating each experiment 5 times.

### D.5 Computational Resources

The experiments were run on a machine equipped with an AMD Ryzen 7 7700X 8 core CPU and 32 GB of RAM.
The comparisons with TVAE and TabDDPM further made use of a GeForce RTX 3060 GPU with 12 GB of VRAM.

## Appendix E Additional Results

### E.1 Computational Effort

In Figure 4 we show the training times for NRGBoost as well as the other methods for the best model selected by hyperparameter tuning.
We do not report the training times for the RFDE method because it is virtually free when compared to the other methods since the splitting process is random and depends only on the input domain. The data itself is only required for computing leaf probabilities which is inexpensive.

Note that the biggest computational cost for training a NRGBoost model is the Gibbs sampling (accounting for roughly 70% of the training time on average).
This could potentially be improved by leveraging higher parallelism than what we used in the experiments (16 virtual cores).

We note also that we believe there is still plenty of margin for optimizing the tree-fitting code used for DEF models. As such, the results presented are merely indicative.

The computational effort for fitting a tree-based model should scale linearly both in the number of samples and the number of features and this is the main source of variation for the training times between datasets.
However, we note that larger datasets can benefit more from larger models (e.g., with a larger number of leaves) which are also slower to train.

Figure: Figure 4: Wall time required to train the best model for each method.
Refer to caption: x4.png

### E.2 MCMC Chain Convergence on MNIST

In Figure 5 we show the convergence of a Gibbs sampler sampling from a NRGBoost model.
In only a few samples each chain appears to have converged to the data manifold after starting at a random sample from the initial model (a mixture between the product of training marginals and a uniform).
Note how consecutive samples are autocorrelated. In particular it can be rare for a chain to switch between two different modes of the distribution (e.g., switching digits) even though a few such transitions can be observed.

Figure: Figure 5: Downsampled MNIST samples generated by Gibbs sampling from a NRGBoost model. Each row corresponds to an independent chain initialized with a sample from the initial model $f_{0}$ (first column). Each column represents a consecutive sample from the chain.
Refer to caption: x5.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: I would like to increase my score according to the comments from other reviewers and the response.

(1) What about the inference time comparison between NRGBoost and other baselines?

(2) Why there is often a performance gap between discriminative GBDTs (i.e., XGBoost here) and NRGBoost? Is it possible to outperform these discriminative baselines from generative modeling paradigm?

(3) How about the discriminative performance of NRGBoost on large regression datasets (e.g., Microsoft, Yahoo dataset in [1]), large-feature-amount one (e.g., Epsilon in [1]) and large-class-number one (e.g., ALOI in [1])?

**Reference**

[1] Revisiting Deep Learning Models for Tabular Data, NeruIPS 2021.
soundness: 3
strengths: I am not expert in the fields of generative models or tree-based models, therefore I think I am not qualified to comment on the paper novelty, while there do exist certain aspects making me impressive.

**Solid formulation, proof & demonstration**: Based on the theory foundation of GBDT, the extension to its generative version NRGBoost is natural and its formulation derivation is clear and convincing. The proposed amortized sampling approach is reasonably designed utilizing the properties of the boosting algorithm. The experiments on tabular discriminative & generative tasks are well conducted with repeated trails, statistical tests and visualization analysis. Besides, complete theoretical and technical Appendix is given, including the computational time comparison for training each baseline.

**Wide application scenarios**: As a generative boosting algorithm, NRGBoost is able to be applied to both discriminative and generative tasks, which distinguishes it from other neural-based generative baselines.

**Reference-worthy exploration on tree-based models**: This paper explores the generative extensions in both Boosting- and Bagging-based tree models, giving the reference for the further research in the community.
summary: This paper explores generative extensions of tree-based models explicitly molding the data density for structured tabular data. Specifically, an energy-based generative boosting algorithm NRGBoost is proposed which is trained to maximize a local second order approximation to the likelihood at each boosting iteration. To save the major training budget caused by approximately estimating the event probability of input data with sampling, the authors further propose an amortized sampling approach by maintaining a fixed-size sample pool with rejection sampling at each round to reduce training time. Apart from designing NRGBoost, the authors also explore bagged ensembles of DET as a generative counterpart to Random Forest. Comprehensive experiments on tabular discriminative and generative (data synthetic efficiency & Discriminator Measure) tasks show NRGBoost can be comparable to the prevailing discriminative GBDTs on prediction tasks as well as has competitive generation quality compared to recent tabular generative neural models, all achieved in one gradient boosting algorithm.
weaknesses: Several weaknesses from application perspective.

(1) From Table 1, in discriminative tasks the performance gap between NRGBoost and traditional discriminative GBDTs (e.g., XGBoost) seems to be relatively more significant as the data scale increases, there is still room for further improvement before application in large-scale discriminative classification tasks.

(2) From the computational effort analysis during training, in Fig. 4, it seems NRGBoost is not computation-economical in the large datasets compared to other neural-network-based generative models, the data scalability of NRGBoost is not sufficiently discussed.

(3) From evaluated dataset information in Table 5, there lacks large regression datasets, which help us to further realize NRGBoost on large-scale regression tasks.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 4
questions: 1. How does the performance compare to DL-based generative models?
2. Can the tradeoffs of NRGBoost vs e.g. DET/DEF be quantified somehow?
3. Is there a concrete use case where NRGBoost is the only/best model?
soundness: 3
strengths: - (clarity) Very well-written paper that is easy to follow.
   I also appreciate the clear presentation and discussion of limitations of the proposed approach.
 - (originality) The proposed boosting framework seems to be novel and provides some refreshing insights.
 - (originality/quality) The overview over existing approaches to use decision trees for generative modelling is extensive and contextualises the work well.
 - (quality) The derivations are explained well and the appendix provides useful additional details.
 - (quality) The experiments seem to have been set up properly with error bars and competitive baselines.
summary: This manuscript presents a framework for using boosting to model unnormalised densities.
The proposed framework is used to build a model using decision trees as weak learners.
To accelerate sampling from the model (and thus training), a combination of Gibbs and rejection sampling is used.
Experiments on a down-sampled MNIST and UCI datasets demonstrate competitive performance to other tree-based generative models.
weaknesses: - (significance) The abstract mentions a comparison with neural network models, but I could not find these in the main paper.
   Furthermore, the authors claim that tree-based generative models are to be preferred over DL-based models
   because discriminative DL methods do not provide competitive performance on tabular data.
   However, no DL methods have been included as baselines in the experiments to confirm this.
   As a result there is no evidence that DL methods would not be able to outperform tree-based generative models.
   PS: DL methods based on normalizing flows have explicit (tractable) densities.
 - (significance) It is not entirely clear how the proposed model can/should be used.
   The discussion mentions the overhead due to sampling, but this overhead is never quantified.
   This makes it hard to get a feeling for whether the (sometimes minor) increase in performance is "worth it".
   Also, it is not entirely clear why shallower trees are preferrable, as implied in the discussion.
 - (significance) The authors put a lot of emphasis on the tractability of computing densities.
   However, there are no clear experiments that illustrate the importance of this feature.
   The conclusion hints at applications enabled by density models,
   but this seems to indicate that there are no real use cases for this model.
   Especially since the results seem to indicate that there are other competitive methods out there.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 2
questions: 1. The paper is hard to follow and haphazardly written. It uses some vague terms (mid-range consumer CPU, Line 88 and the best generative models, Line 89). Minor details Line 41: please avoid abbreviation like don't, can't.

2. The paper proposes generative model for tabular data in the abstract, but conducts experiments with vision data like MNIST. I am not sure how downsampling makes MNISt relevant to their setting. I did not understand what is going on in Figure 2. There is no description of the data generation process or distribution. 

3. The equations are incoherent and sometimes most probably wrong. For example, I do not see the purpose of Equation 2. It was never used later. The math is incoherent and written in such a way that it is hard for me to evaluate their correctness. 

4. Table 1: the authors say they bold the best performing algorithms, but they kept their algorithm bold always although it did not perform the best.

5. The paper has bad structure. Related works comes at Section 5 near the end of the paper. There is no flow diagram or pseudocode for their algorithm. 

Overall, I think the paper is far from being a coherent and legit study.
soundness: 3
strengths: The paper is confusing and hard to follow. I was unable to find any strength after investing a long time.
summary: The paper proposes a tree-based approach to model tabular data density which they claim outperforms traditional generative approaches.
weaknesses: The paper is poorly written, incoherent and uses inconsistent experiments.

---

confidence: 4
contribution: 4
initial_rating: 6
presentation: 2
questions: Questions

- The methods section and experiments devoted to DEF seems out of place, given that it is seemingly unrelated to NRGBoost and given that it performs worse. It is nice to tell more people about DETs and how they can be applied and improved, but I suggest perhaps omitting it entirely and submitting it as an ICLR Blogpost or similar venue.

- I realize that NGBoost does not support nonparametric probabilistic prediction, but I think seeing its results in Table 1 would still be helpful. 

Minor suggestions for improved readability:

Replace: We therefore estimate these quantities, with recourse to empirical data for P(X), and to samples approximately drawn from the model with MCMC.
With: We therefore estimate these quantities, with recourse to empirical data for P(X) and to samples approximately drawn from the model with MCMC.

Replace: Because even if the input space is not partially discrete, f is still discontinuous and constant almost everywhere we can’t use gradient based samplers and therefore rely on Gibbs sampling instead.
With: Because, even if the input space is not partially discrete, f is still discontinuous and constant almost everywhere, we can’t use gradient based samplers; we therefore rely on Gibbs sampling instead.
soundness: 2
strengths: 1. The method is original and well-motivated, and the approach for incorporating rejection sampling within the boosting process to improve efficiency was an excellent contribution towards making energy-based tabular modeling practical.

2. The paper is well-written, with the method and experiments clearly explained. Showing how the model progressively models the dataset over boosting iterations in Figure 2 gave helpful insight into how the method works.

3. The experiments on downsampled MNIST data were very nice to see, and convincing that NRGBoost goes beyond current SotA (ForestFlow) [2] in generation in certain settings, besides also offering density estimation.
summary: NRGBoost is an energy-based method for tabular data density estimation and generation using gradient boosted decision trees (GBDTs), in contrast to recent diffusion-based generative methods which cannot give density estimates.
weaknesses: 1. [UPDATE: fully addressed in rebuttal] For single-variable inference in Section 6.1, using R^2 as an evaluation metric seems strange, compared to using CRPS, RMSE, and MAE, as employed in [3]. Furthermore, the choice of baseline methods seems incomplete. For example, because Forest-Flow supports conditioning on covariates, one could also generate (say) 100 Forest-Flow conditional samples, then compute the mean and evaluate this.

2. [UPDATE: mostly addressed in rebuttal] The tabular datasets used for evaluating NRGBoost in Section 6.2 were nonstandard. Most recent papers either evaluate on the TabCSDI [1] benchmark of 6 datasets or the ForestFlow [2] benchmark of 27 datasets. Because using one's own choice of datasets does open the possibility for cherry-picking, one should either use a preexisting benchmark or justify one's selection of datasets. 

3. [UPDATE: fully addressed in rebuttal] Some relevant works are missing from related work (and perhaps the experiments). It would be good to discuss how this compares against a previous energy-based modeling method, TabEBM [4], albeit one which uses TabPFN instead of GBDTs. And there is UnmaskingTrees [5] which uses GBDTs to perform both generation and also density estimation via autoregression and recursive partitioning.
  
[1] Zheng, S., & Charoenphakdee, N. (2022). Diffusion models for missing value imputation in tabular data. arXiv preprint arXiv:2210.17128.

[2] Jolicoeur-Martineau, A., Fatras, K., & Kachman, T. (2023). Generating and Imputing Tabular Data via Diffusion and Flow-based Gradient-Boosted Trees. arXiv preprint arXiv:2309.09968.

[3] Beltran-Velez, N., Grande, A. A., Nazaret, A., Kucukelbir, A., & Blei, D. (2024). Treeffuser: Probabilistic Predictions via Conditional Diffusions with Gradient-Boosted Trees. arXiv preprint arXiv:2406.07658.

[4] Margeloiu, A., Jiang, X., Simidjievski, N., & Jamnik, M. (2024). TabEBM: A Tabular Data Augmentation Method with Distinct Class-Specific Energy-Based Models. arXiv preprint arXiv:2409.16118.

[5] McCarter, C. (2024). Unmasking Trees for Tabular Data. arXiv preprint arXiv:2407.05593.

Meta-review: The paper extends Gradient Boosted Decision Trees and Random Forests to generative modeling tasks. This paper went through significant changes during the review period, in a persistent effort from the authors to address comments from the reviewers, who eventually increased their original scores. My overall impression from the paper itself and the reviews and discussion is also positive. Modelling density is not a natural task for forests and this paper adds value to the small literature on this topic. It also may attract interest to the use of forests in generative modelling. I do concur with a reviewer that it is not clear what the exact use case is for NRGBoost, though.

Decision: Accept (Poster)