NRGBoost: Energy-Based Generative Boosted Trees
Authors: João Bravo
Venue: ICLR 2025 Poster (2025)
Keywords: Energy-Based Models
Generative Models
Gradient Boosting
Tabular Data
Primary area: generative models
TL;DR: We explore extensions of Gradient Boosted Decision Trees and Random Forests to generative modeling tasks

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 libraries 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.
Code is available at https://github.com/ajoo/nrgboost.

Published as a conference paper at ICLR 2025
NRGBoOST: ENERGY-BASED GENERATIVE BOOSTED TREES
João Bravo
Feedzai joao.g www.ooout.oook.....oooooo com
ABSTRACT
Despite the rise 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 libraries like XGBoost. We show that, despite producing generative model capable of handling inference tasks over any input variable, proposed algorithm can achieve similar discriminative performance to GBDT on a number real world tabular datasets, outperforming alternative generative approaches. At the same time, show that it is also competitive with neural-network- based models for sampling. Code is available at https: http//gtttuueat com/ nrgboss....
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 a potential avenue for bypassing privacy concerns when sharing data. Estimating the data density, p(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 single target variable, density models can be used more flexibly for inference over different variables or for cut-of-diisributtinn detection. They also handle inference with missing data in principled way by marginalizing over unobserved variables.
The development of generative models for tabular data has mirrored its progression in computer vision, with many its deep learning approaches being adapted to the tabular domain (Jordon et al., 2018; Xu al., 2019; Engelmann Lessmann, 2021; Fan al., 2020; Zhao al., 2021; Kotelnikov et al., 2023). Unfortunately, these methods are only useful for sampling as they either do not model the density explicitly can not evaluate it efficiently due intractable marginalization over high- nigh-dimensional latent variable spaces. Furthermore, despite growing in popularity, deep learning 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 nuural-networkee methods many real world datasets (Grinsztajn et al., 2022; Borisov et al., 2022a).
While there have been recent efforts to extend the success of tree-based models to generative modeling (Correia et al., 2020; Wen & Hang. 2022; Nock & Guillame- Bert, 2022: Watson et al., 2023; Nock & Guillame-Ber Bert, 2023; SeriiccuurMaattnnen
al., 2024; McCarter, 2024), direct extensions of Random Forests (RF) and Gradient Boosted Decision Trees (GBDT) are still missing. We try to address this gap, 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:

---

Published as a conference paper at ICLR 2025
Sample Pool
Training Data
Current Model
2
Gibbs Sampling
q
P
Fit Tree
5f
Weak Learner
Drop Samples
Add
Figure I: Overview of an NRGBoost training iteration: 1) Draw new samples from an ensemble of trees representing a current energy function, f. and add them to the sample pool. 2) Fit a new weak learner (tree) representing a model update, 8f. 3) Update the model. 4) Use rejection sampling to discard samples from the sample pool that conform poorly to the new model.
We propose NRGBoost, a novel energy-based generative boosting model that is trained to maximize local second-order approximation to the ligeliiel at each boosting round.
We propose an amortized sampling approach that significantly reduces the training time of NRGBoost, which, like other energy-based models, is often dominated by the Markov chain Monte Carlo sampling required for training.
We explore the use of bagged ensembles of Density Estimation Trees (Ram & Gray, 2011) with feature subsampling generative counterpart to Random Forests.
The longstanding popularity 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, minutes on a typical consumer CPU and achieves similar discriminative performance to a standard GBDT model, while remaining competitive with other generative models for sampling.
implementation of NRGBoost can be trained in few
ENERGY BASED MODELS
An Energy-Base Based Model (EBM) parametrizes the logarithm of a probability density function directly (up to an unspecified normalizing constant):
(x) = exp (f(x))
Z[S]
Here, f(x) : X + is a real function over the input domain.' 1 We will avoid introducing any parametrization and instead treat the function E F(X), lying in appropriate function space over the input domain, as our model parameter directly. Zif] = ExEx exp (x(),, known as the partition function, then functional of giving the necessary normalizing constant.
This is the most flexible way one could represent a probability density function, making essentially compromises its structure. It can seamlessly handle input domains comprising mixture of continuous and categorical features, making it a natural choice for tabular data. The downside to this is that, for most interesting choices of F. computing the normalizing constant intractable which makes training these models difficult. Their unnormalized nature, however, does not prevent EBMs from being useful number of applications besides sampling, such as performing joint inference over small enough groups of variables. For example, by splitting the input domain into
We will assume that X is finite and discrete to simplify the notation and exposition but everything is applicable bounded continuous input spaces, replacing the with integrals as appropriate.
2

---

Published as a conference paper at ICLR 2025
set of observed variables, X and unobserved variables, X we have
exp (f (X xo))
(xxx)(xo) - Ex; exp (f(x,, Xa))
which only involves normalizing (i.e., computing a softmax) over all of the possible values of the unobserved variables. EBMs can also handle inference with missing data principled manner by marginalizing over unobserved variables. Furthermore, for anomaly out-of-distribution detection, knowledge of the normalizing constant is not necessary.
One common way train energy- model to approximate data generating distribution, p(x), to minimize the Kullback- divergence between and q,, or equivalently, maximize the expected log- likeeeildoood functional:
L Ex-p log (xx(xx)))x(((Xe (x) - log Z[
Ex-pp
(3)
This optimization is typically carried out by gradient descent over the parameters of f. However, due to the intractability of the partition function, one must rely on Markov chain Monte Carlo (MCMC) sampling to estimate the gradients (Song & Kingma, 2021).
NRGBoOST
Expanding the increase in log-likelihen liliiiiiihodd in Equation 3 due to a variation öf around an energy function f up to second order (see Appendix A), we have
Lf + o L[f] so Ex)p - Ex 8 of(x) Varx 8f(x) = [
(4)
The 5ʃ that maximizes this quadratic approximation should thus have a large positive difference between the expected value under the data and under ęs. while having low variance under df- Note that, just like the original log-likelihood, this Taylor expansion is invariant to adding a constant to 8f (i.e., L c| ر f] for constant function c). This means that, in maximizing Equation 4, we can restrict our attention to functions that have zero expectation under gy. in which case we can simplify g as
1 Exxqq o (x).
A = of
(5)
2
We can thus formulate a boosting algorithm that, at each boosting iteration, t, improves upon a current energy function, ft, by finding an optimal step, off, that maximizes A
ōʃ arg max AT,,155f1
SfER.
where H. is an appropriate space of functions (satisfying Ex-ay, fx( = if Equation is used). The solution to this problem can be interpreted Newton step in the space of energy functions. Because, for an EBM, the Fisher information matrix for the energy function and the Hessian of the expected logelikelihood are the same, (see Appendix A). This approach is analogous to the second-order step implemented in modern discriminative gradient boosting libraries such XGBoost (Chen & Guestrin, 2016) and LightGBM (Ke et al., 2017) and which can be traced back to Friedman et al. (2000).
also interpret this solution as a natural gradient step
In updating the current iterate, fr-1 ft + ar Öfř. we scale f. by an additional scalar step-size Or. 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 Oa via a line search based on the original log- likelihood. Common practice in discriminative boosting, however, to interpret this step-size size as regularization parameter and to select a fixed value in J0, 1]. with (more) smaller steps typically outperforming fewer larger ones when it comes to generalization. We choose to adopt a hybrid strategy, first selecting optimal step-size by line search and then shrinking it by 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, fur we can choose the logarithm of any probability distribution over X, as long as it is easy to evaluate. Sensible choices 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 toy 2D data density.
uniform distribution (i.e., = 0), the product of
3

---

Published as a conference paper at ICLR 2025
Model, t=
Model, L 3
Model, t= 10
Model, t= 100
Data
Figure 2: Density learned by NRGBoost different boosting iterations (1, 3, 10 and 100), starting from uniform distribution. The data distribution is depicted on the right (see Appendix D.1 for details). Weak learners are piecewise constant functions given by binary trees with 16 leaves.
3.1 WEAK LEARNERS
As a weak learner we will consider piccewise constant functions defined by binary trees over the input space. Letting -|=1+ X; = X' be the partitioning of the input space induced by the leaves of binary tree whose internal nodes represent a split along one dimension into two disjoint partitions, take H. as the set of functions such
J 8f(x) = Σ wy wwyxx (x).
5=1
Here, 1x denotes the indicator function of subset X and wy are values associated with each leaf j [1..J]. In standard decision tree, these values would typically encode an estimate of ply|x E Xj). with 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 the öf functions are thus parametrized by the values, Wj. as well as the structure of the tree and the variables and thresholds for the split at each node, which ultimately determine the Xj. We omit these dependencies for brevity.
By substituting the definition from Equation into our objective (Equation 5) we get the following optimization problem to find the optimal binary tree:
1 -(3)))==33((x))
w.
41's yours
(8)
s.t. Σ
[w, =0
where P(X,) and (X,) denote the probability of the event X, under the respective distribution and the constraint ensures that of has zero expectation under gs. With respect to the leaf weights, this is quadratic program whose optimal solution and objective values are respectively given by
P2 -1) 2 Qj(X,)
P(X,) w (X,,
1
Σ
1,
ALj (X 1. X,)
(9)
Because carrying out the maximization of this optimal value over the tree structure that determines the X, is hard, we approximate its solution by greedily growing tree that maximizes it when considering how to split each node individually. A parent leaf with support Xp is thus split into two child leaves, with disjoint support X₁. XR = -xp. so maximize over all possible partitionings along a single dimension, P (Xr), the following objective:
(2)) +
(²(XR)
P2(Xp)
max XxEXXE((X (x ()(xn) (x(Pp)
(10)
Note that, when using parametric weak learners, computing a second-order step would typically involve solving linear system with full Hessian. As see, this is not the when the
4

---

Published as a conference paper at ICLR 2025
weak learners are binary trees, where the optimal value assign to leaf does not depend any information from other leaves and, likewise, the optimal objective value a sum of terms, each depending only information from single leaf. This would not have been the case had tried to optimize the log-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 (X)),, which is infeasible due to the intractable 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 not use gradient based samplers. We therefore rely on Gibbs sampling, which only requires evaluating each f. along one dimension at a time, while keeping all others fixed. Importantly, this can be computed efficiently for a tree by traversing it only once. However, since our energy function at boosting iteration / is given by a sum of f trees, this computation scales linearly with the iteration number. As a result, the total computational cost of sampling, over the course of training model, scales quadratically in the number of boosting iterations, 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 current model distribution. It could, therefore, be helpful keep some of those samples, especially those that conform best the current model. Rejection sampling allows us to do just that. The boosting update in terms of the densities takes the following multiplicative form:
(11) Here, k2 is an unknown multiplicative constant and since of is given by tree, we can easily bound the exponential factor by finding the leaf with the largest value. To sample from 2+11 we can
(+)(x)) kt q: (x() exp (a))f((x))
therefore use the previous model, qt. as a proposal distribution for which we already have a set of samples. We retain each of these previous samples, with an acceptance probability of
Paccept (x) exp (os (5f((x) max (x)(
(12)
noting that knowledge of the normalizing constant, kir, is not necessary for its computation.
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 updated model using Gibbs sampling. Figure 1 depicts representation of the training loop with this amortized sampling approach. Note that qo is typically a simple model for which partition X) and efficiently draw exact samples. As such, no samples are required for the first round of boosting. For the second round, can initialize the sample pool by drawing exact samples from Y with rejection sampling, using 40 as a proposal distribution.
both directly evaluate the desired quantities (i.e., Qx(X) for a given
This amortized sampling approach works better when either the range of öf the step sizes oz are small, as this leads to larger acceptance probabilities. Note that, in practice, it can be helpful independently refresh fixed fraction of samples, Prefresh at each round of boosting in order to encourage a more diverse set of samples between rounds. This can be accomplished by keeping each sample with a probability Poccept (1) Prefresh) instead. Please refer to Algorithm in Appendix A.3 for a high level algorithmic description of the training process with amortized sampling.
3.3 REGULARIZATION
The simplest way to regularize a boosting model is to stop training when overfitting is detected by monitoring suitable performance metric validation set. For NRGBoost this could be the increase in likeliiel at each boosting round (see Appendix A.1) However, estimating this quantity would require drawing additional validation samples from the model. A viable alternative validation strategy which needs no additional samples is to simply monitor a discriminative performance metric. This amounts to monitoring the quality of o (3ra xx;) instead of the full gx().
5

---

Published as a conference paper at ICLR 2025
Besides early stopping, the binary trees themselves can be regularized by limiting the depth total number of leaves in each tree. Additionally, we can rely on other strategies, such as disregarding splits that would result in leaf with too little training data, P(X), model data, (x),,, volume, V(X). or too high of a ratio between training and model data, P(X)/Q,(x) We found the latter to he 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 training likelihood at each boosting round.
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 (2002); Tolstikhin et al. (2017). (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 dog-density) as additive sum of weak learners. Welling al. (2002), in particular, also approaches boosting from the point of view of functional optimization of the log-likelihood or the logistic loss of an energy-based model. However, it relies on first order local approximation of the objective given that it focuses on parametric weak learners, such as restricted Boltzmann machines, for which a second-order step would be impractical.
Another more direct multiplicative boosting framework was first proposed by Tu (2007). At each boosting round, a discriminative classifier trained distinguish between empirical data and data generated by the current model, by estimating the likelihood ratio p(x)/9;(x) This estimated ratio is used as direct multiplicative factor to update a current model, q.- In ideal conditions, this greedy procedure would converge in single iteration (when using step size of 1). While Tu (2007) does not prescribe a particular choice of classifier, Grover & Ermon (2018) 8) proposes a similar concept where the ratio is estimated based on variational bound for an f f-divergence and Cranko & Nock (2019) 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, cxp (az (())/ instead of (rx))//.(x))) directly. In Appendix C we explore the differences between NRGBoost and this approach, when it is adapted to use trees as weak learners.
Tabular Energy-Based Models Ma et al. (2024) proposes reinterpreting the logits of a TabPFN classifier (Hollmann et al., 2022) as an energy function over the input space. Similarly, Margeloiu et al. (2024) relies on auxiliary binary classification tasks, reinterpreting the logits of the classifiers as class-conditional energy functions. These approaches, while not as principled as maximum likelihood, avoid the need for MCMC sampling and were shown to work well empirically. However, their reliance on TabPFN can limit their applicability to larger datasets. Furthermore, both approaches rely on having a categorical target variable, which further limits their applicability and could introduce an inherent bias given that this variable is not treated in the same way as the others.
TreeeBased Density Models Density Estimation Trees (DET) were proposed by Ram & Gray (2011) as alternative to histograms and kernel density estimation. They model the density function as constant value over the support of each leaf in a binary tree, (x) Σ V(X;) 1x, (x), with V(X) denoting the volume of X. The partitioning of the imput space is determined directly by greedy optimization of a generative objective (ISE). Despite allowing for efficient exact sampling, DETs have received little attention as generative models. Appendix B, we explore extensions of this method to maximum likelihood estimation. In our experiments, also test the natural idea to ensemble DET models through bagging, leading to a generative algorithm analogous to Random Forests.
P(X,)
Other authors have proposed similar tree-based density models (Nock & Guillame-Bert, 2022) or additive mixtures of tree-based models (Correia et al., 2020; Wen & Hang. 2022; Watson et al., 2023). Distinguishing features of some of these alternative approaches are: (i) Not relying on density estimation goal to drive the partitioning of the input space. Correia et al. (2020) leverages a standard discriminative Random Forest, therefore giving special treatment to particular input variable whose

---

Published as a conference paper at ICLR 2025
conditional estimation drives the choice of partitions and Wen & Hang (2022) proposes using mid- point random tree partitioning. (ii) Leveraging more complex models for the data in the leaf of tree instead of a uniform density (Correia et al., 2020; Watson et al., 2023). This approach can allow for the use of trees that are more representative with a smaller number of leaves. (iii) Nock & Guillame-Ber (2022) and Watson et al. (2023) both propose generative adversarial frameworks where the generator and discriminator are both a tree or an ensemble of trees respectively. While this also leads an iterative approach, unlike with boosting or bagging the new model trained each round doesn't it add to the previous one but replaces it instead.
Nock Guillame- Bert (2023) proposes a different ensemble approach where trees do not have their own leaf values that are added or multiplied to produce a final density, but instead collectively define the partitioning of the input space. To train such models the authors propose a framework where, rather than adding new tree to the ensemble at every iteration, the model is initialized with fixed number of tree root nodes and each iteration adds a split to an existing leaf node. Jolicoeur- Marttnnea et al. (2024) proposes diffusion and conditional flow matching approaches where a tree-based model (e.g., GBDT) is used to learn the score function or vector field. However, these approaches do not support efficient density estimation. McCarter (2024) leverages the autoregressive framework of TabMT (Gulati & Roysdon, 2024), replacing transformers with XGBoost based predictors. While this approach is capable of density estimation, it requires training one predictor (consisting of one or more XGBoost models) per input variable. Furthermore, these models are trained on datasets created from multiple copies of the original dataset with different missing data patterns, which could limit their applicability to larger datasets.
5 EXPERIMENTS
We evaluate NRGBoost on five tabular datasets from the UCI Machine Learning Repository (Dheeru & Karra Taniskidou, 2017): 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 al., 201 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.I.
We split our experiments into two sections, the first to evaluate the quality of density models directly on a Lannle-variabll inference task and the second to compare the performance of our proposed model on sampling use cases against more specialized models.
5.1 SINGLE- VARIAA-A INFERENCE
We test the ability of generative model, trained to learn the density over all input variables, 9(xx,, to infer the value of a single one given the others. Specifically, we evaluate the quality of its estimate of q(xx/x=;) For this purpose we pick x; = y the original target of the dataset, noting that the models that we train do not treat this variable in any special way, except when selecting the best model in validation. As such, we would expect that a model's performance in inference over this particular variable is indicative of its strength on any other lingle-varibbe inference task and also indicative of the quality of the full x(x), from which the conditional probability estimate is derived.
We use XGBoost (Chen & Guestrin, 2016) as a baseline for what should be achievable by strong discriminative model. Note that this model is trained to maximize the discriminative likelihood, Ex-p log ((:; x)) not wasting model capacity in learning other aspects of the full data distribution. Because generative models provide an estimate of the full conditional distribution for the regression datasets, rather than single point estimate like XGBoost, we also include NGBoost (Duan et al., 2020) as another discriminative baseline that can estimate the full q(z, x) on these datasets. Note, however, that NGBoost relies a parametric assumption about this conditional distribution, that needs to hold for any x.;. We use the default assumption of Normal distribution for all datasets.
We compare NRGBoost against bagging ensemble of DET models (Ram & Gray, 2011), which we call Density Estimation Forests (DEF), trained for minimizing either ISE or KL divergence (see Appendix B). We also include two other tree-based generative baselines: RFDE (Wen & Hang, 2022)
7

---

Published as a conference paper at ICLR 2025
Table 1: Discriminative performance of different methods at inferring the value of a target variable. We use R² for regression tasks, AUC for binary classification and accuracy for multiclass classifi- cation. The reported values are means and standard errors over 5 cross-validatin generative method for each dataset is highlighted in bold and other methods that are not significantly worse (as determined by a paired (-test at 95% confidence level) are underlined.
folds. The best
R²² 1
AUC 1
Accuracy 1
AB
CH
PR
AD
MBNE
MNIST
CT
XGBoost
0.552 +0,05
0.849 =0.009
0.678 =0.004
0.927 +0,000
0.987 +0,000
0.976 =0002
0.971 =0.001
NGBoost
0.546 +0,40
0.829 +0.009
0.621 +0.005
RFDE
0.071 20.096
0.340
20.004
0,059 20.007
0.862 20.002
0.668 20.008
0.302 20010
0.679 20.002
ARF
0,531 20.032
0.758 20.009
0,591 +0.000
0,893 20,002
0,968 20.001
0,938 20.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 -00t1
0.790 -0.003
DEF (KI.)
0.482 +0.027
0.801 +0.008
0.639 +0.004
0.892 +0.001
0.939 +0.001
0.487 -00
0.852 +0.002
NRGBoost
0.547 20.036
0.850 s0.011
0.676
0.920 a0.001
0.974 s0.001
0.966 =0.001
0.948 20.001
Due to the discrete nature of the MNIST dataset and the fact that the ARF algorithm tries fit continuous distributions for numerical variables at the leaves, we could not obtain a reasonable density model for this dataset.
and ARF (Watson al., 2023). The former allows us to gauge the impact of the guided partitioning used in DEF models over random partitioning of the input space.
We use random search to tune the hyperparameters of XGBoost and NGBoost 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.
We find that NRGBoost outperforms the remaining generative models (see Table I), even achieving comparable performance to XGBoost on the smaller datasets and with a small gap on the three larger ones (MBNE, MNIST and CT). Furthermore, in Appendix E.3 we show that, for the same tasks but in the presence of missing covariate, NRGBoost can outperform XGBoost relying on popular imputation strategies due to its principled handling of missing data through marginalization.
5.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. In addition to the previous ARF and DEF (KI.) baselines, we compare against TVAE (Xu et al., 2019) and TabDDPM (Kotelnikov et al., 2023), two lere-learming-basee generative models, as well as Forest-Flow (Jolicoeur- Martineau et al., 2024), a tree-based conditional flow matching model.
Molicoeur-Martineau
Machine Learning Efficiency Machine learning (ML) efficiency has been a popular way to measure the quality of generative models for sampling (Xu et al., 2019; Kotelnikov et al., 2023; Borisov et al., 2022b). relies on using samples from the 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 5.1. In fact, the discriminative model is flexible enough, one would expect it to recover the generator's ' ((y(x),, and therefore its performance, in the limit where infinite generated data used to train it.
We use XGBoost as the discriminative model and train it on the same number of synthetic samples as the original training data. For the density models, we generate the samples from the best model found in the previous section. For TVAE and TabDDPM, by evaluating the ML efficiency in the real validation set (full details of the hyperparameter tuning provided in Appendix D.3). Note that this leaves these models at a potential advantage since the hyperparameter selection is based on the metric that is being evaluated.
select their hyperparameters
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 perf-pero model depending on the dataset (with two inconclusive cases), and NRGBoost
8

---

Published as a conference paper at ICLR 2025
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 for evaluating its performance. Reported values are the averages and standard errors over 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 l-test a 95% confidence level) are underlined.
R²² 1
1
Accuracy t
CH
PR
MBNE
MNIST
CT
Real Data
0.554
0.838
0.682
0.927
0.987
0.976
0.972
TVAE
0,483 20.006
0.758 20.005
0.365 20.005
0.898 20.001
0.975 20.000
0.770
0,750
TabDDPM
0.539 20.018
0.807 20.005
0.596 20.007
0,910 a0.001
0,984 10,000
0.818 20.001
Forest-Flow
0.418 20.019
0.716 20.003
0,412 20,009
0.879 a0.003
0.964 20.001
0.224 20.010
0,705 20.002
ARF
0.504 +0.020
0.739 +0.003
0.524 +0.00s
0,901 +0.001
0.971 +0.001
0.908 +0.002
0.848 +0.002
DEF (KL)
0.450 20.013
0.762 20.006
0.498 20.007
0.892 a0.001
0.943 #0.002
0.230 20.028
0,753 20.002
NRGBoost
0.528 20.016
0.801 20.001
0.573 20.008
0.914 +0.001
0.977 20.001
0.959 20.001
0,895 20.001
Training Data
NRGBoost
DEF (KL)
ARF
Laitude
Longitude
Longitude
Longitude
Longitude
TahDDPM
Forest-Flow
TVAE
fafifuce
Longitude
Longitude
Longitude
Figure 3: Joint histogram for the latitude and longitude for the California Housing dataset.
is never ranked lower than second on any dataset. We also note that, despite our best efforts with additional manual tuning, we could not obtain reasonable model with TabDDPM on MNIST.
Discriminator Measure Similar Borisov al. (2022b) test the capacity of 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 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 test data in Table 3, which shows that NRGBoost outperforms other methods except on the PR (where the result is inconclusive) and MBNE datasets.
In Appendix E. we report Wasserstein distance estimate computed using similar setup to Jolicoeur-Martieauu et al. (2024) as an alternative measure the statistical distance between real and model distributions. This evaluation also shows NRGBoost ranking the best average across datasets. Qualitatively speaking, its samples also look visually similar to the real data both the California and MNIST datasets (see Figures 3 and 4).

---

Published as a conference paper at ICLR 2025
Training Data
NR GRoos
Forest-Flow
ARF
TVAE 9 0
2
2
O
1
4 4
3
8
7
F4
5
H 0 4 8
4 8
7
6
3
8 8
3
2 4 7
9
3
6
3
8 4 3
4 3 3 7
4
5 4
9 5
5 7 6
8
3
5 1
4 9 5
5 4
&
7 1
37
0
7
7
4
$
3
3
3 6 9
0
0
4
4
8
2
$
/
3
7 A
7
6
9 3
Figure 4: Downsampled MNIST samples generated by the best generative models on this dataset. Despite being a simple dataset that would pose no challenges to image models, it is hard for tabular generative models due to the high dimensionality and complex structure of correlations between features. We find NRGBoost to be the only tabular model that is able to generate passable samples.
Table 3: AUC of an XGBoost model trained to distinguish real from generated data (lower is 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 a 95% confidence level) are underlined.
AB
CH
PR
AD
MBNE
MNIST
CT
TVAE
0,971 20.004
0.834 10.006
0.940 20.002
0.898 20.001
1.000 20,000
1.000 20,000
0.999
TabDDPM
0,818 20.015
0.667 20.005
0.628 20.004
0.604
10.002
0.789 10.002
0,915 20.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 20.005
0.973 20.004
0.795 20.008
0,992
40,000
0,998 20.000
1.000
20.000
0,989
20.001
DEF (KL)
0.823 20.013
0.751 20.008
0,877 20.002
0,956 20.002
1 1.000 20.000
1.000 20,000
0,999 20,000
NRGBoost
0.625 20.017
0.574 20.012
0.631 20.006
0.559 20.003
0.993 20.001
0.943 20.003
0.724 20.006
DISCUSSION
While additive tree ensemble models like DEF require sampling to train and are easy to sample from, we find that, in practice, they require very deep trees to model the data well and still offer subpar performance for sampling and density estimation.
In contrast, NRGBoost was able to model the data better, while using fewer and shallower trees. Its main drawback is that can only be sampled from approximately using more expensive MCMC and also requires sampling during the training process. Our fast Gibbs sampling implementation. coupled with our proposed amortized sampling approach, were able to mitigate the slow training, making these models more practical. Unfortunately. they remain cumbersome to use for sampling, as autocorrelation between samples from the same Markov chain makes them inefficient in scenarios requiring independent samples. We argue, however, that unlike in image or text generation, where fast sampling is necessary for an interactive user experience, this can less of concern when generating synthetic datasets, where the one-time cost of sampling can be less important than faithfully capturing the nata-generating distribution.
One advantage of the energy-based approach that we did not explore is that it allows for arbitrary conditional sampling. requiring only that one clamps the values of conditioning variables during Gibbs sampling. In contrast, other methods, such diffusion- based approaches, less flexible, since conditioning variables typically need to be predetermined at training time.
7 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, particularly in exploring other applications beyond sampling that are enabled by density models.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Kristin Blesch and Marvin N. Wright. arfpy: A python package for density estimation and generative modeling with adversarial random forests, 2023.
Vadim Borisov, Tobias Leemann, Kathrin SeBler, Johannes Haug. Martin Pawelczyk, and Gjergji Kasneci. Deep neural networks and tabular data: A survey. IEEE transactions on neural networks and learning systems, 2022a.
Vadim Borisov, Kathrin Seßler, Tobias Leemann, Martin Pawelczyk, and Gjergji Kasneci. Language models are realistic tabular data generators. arXiv preprint arXiv:22/0.066 0620000, 2022b.
Tianqi Chen and Carlos Guestrin. XGBoost: A Scalable Tree Boosting System. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 785-794, August 2016. doi: 10.1 00.1144/299977222339933 arXiv:1603.00 100.00274 [cs].
Alvaro Correia, Robert Peharz, and Cassio de Campos. Joints in random forests. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Curran Associates, Inc., 2020.
Systems, volume 33. pp. 11404-15 11415.
Zac Cranko and Richard Nock. Boosted Density Estimation Remastered. In Proceedings of the 36th International Conference on Machine Learning, pp. 1416-1425. PMLR, May 2019. ISSN: 2640-3498.
Li Deng. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine, 2966)1141-14,, 2012. 2012.
Dua Dheeru and Efi Karra Taniskidou. UCI machine learning repository, 2017.
Tony Duan, Avati Anand, Daisy Yi Ding, Khanh K Thai, Sanjay Basu, Andrew Ng, and Alejan- dro Schuler. Ngboost: Natural gradient boosting for probabilistic prediction. In International conference on machine learning, pp. 2690-2700. PMLR, 2020.
Justin Engelmann and Stefan Lessmann. Conditional wasserstein gan-based oversampling of tabular data for imbalanced learning. Expert Systems with Applications, 174: 14582, 2021.
Ju Fan, Junyou Chen, Tongyu Liu, Yuwei Shen, Guoliang Li, and Xiaoyong Du. Relational data synthesis using generative adversarial networks: design space exploration. Proceedings of the VLDB Endowment, 13(12):1962-199) August 2020. ISSN 2 2150-8097. doi: 10.14778/3407790. 3407802.
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Additive logistic regression: statistical view of boosting (With discussion and a rejoinder by the authors). The Annals of Statistics, 28(2): 337-407, April 2000. ISSN 0090-5364, 2168-8966 8966. doi:
Publisher:
Institute of Mathematical Statistics.
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-37,, 2007.
Yury Gorishniy, Ivan Rubachev, Valentin Khrulkov, and Artem Babenko. Revisiting deep learning models for tabular data. Advances in Neural Information Processing Systems, 34:18932-1 18943, 2021.
Léo Grinsztajn, Edouard Oyallon, and Gaël Varoquaux. Why do tree-based models still outperform deep learning on typical tabular data? Advances in neural information processing systems, 35: 507-520, 2022.
Aditya Grover and Stefano Ermon. Boosted generative models. In Proceedings of the AAAI Conference Artificial Intelligence, volume 32, 2018.
Manbir Gulati and Paul Roysdon. Tabmt: Generating tabular data with masked transformers.
Advances in Neural Information Processing Systems, 36, 2024.
11

---

Published as a conference paper at ICLR 2025
Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathamiel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gerrard-Marchan Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585(7825)357-362, September 2020. doi: 10. 10.10388/41186-000-22000
Noah Hollmann, Samuel Müller, Katharina Eggensperger, and Frank Hutter. 'Tabpfn: transformer that solves small tabular classification problems a second. arXiv preprint arXiv:2207.0/ 848, 2022.
Alexia Selicoeur-Martinee Kilian Fatras, and Tal Kachman. Generating and imputing tabular data via diffuusin and flow-based gradient-boosted trees. In International Conference on Artificial Intelligence and Statistics, pp. 1288-1296. PMLR, 2024.
James Jordon, Jinsung Yoon, and Mihaela Van Der Schaar. Pate-gan: Generating synthetic data with differential privacy guarantees, In International conference on learning representations, 2018.
Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie- Yan Liu. LightGBM: A Highly Efficient Gradient Boosting Decision Tree. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Akim Kotelnikov, Dmitry Baranchuk, Ivan Rubachev, and Artem Babenko. Tabddpm: Modelling tabular data with diffusion models. In International Conference on Machine Learning, pp. 17564- 17579. PMLR, 2023.
Junwei Ma, Apoory Dankar, George Stein, Guangwei Yu, and Anthony Caterini. Tabpfgen-taabl -tauuar data generation with tabpfn. arXiv preprint 2XXi:224060022666 2024.
Andrei Margeloiu, Xiangjian Jiang, Nikola Simidjievski, and Mateja Jamnik. Tabebm: A tabular data augmentation method with distinct class-specific energy-based models. arXiv preprins arXiv:2409.1611 2024.
Calvin McCarter. Unmasking trees for tabular data. arXiv preprint arXiv:2407.0552. 2024.
Richard Nock and Mathieu Guillame-Ber Bert. Generative trees: Adversarial and copycat. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Pro- ceedings of the 39th International Conference Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 16906-16951. PMLR, 17-23 Jul 2022.
Richard Nock and Mathieu Guillame-Bert. Generative forests. arXiv preprint arXiv:2308.03648. 2023.
Melissa E. O'eei... Peg: A family of simple fast space-efficient statistically good algorithms for random number generation. Technical Report -CCCCCS--014,005, Harvey Mudd College, Claremont, CA, September 2014.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830 2011.
Parikshit Ram and Alexander G. Gray. Density estimation trees. Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 627-635. San Diego California USA, August 2011. ACM. ISBN 978-1-4503-0081377 doi: 10.1 15.1145/200040. 2020507.
Saharon Rosset and Eran Segal. Boosting Density Estimation. In Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002.
Yang Song and Diederik P Kingma. How to train your energy-based based models. arXiv:2101.088 032888. 2021.
arXiv preprint
12

---

Published as a conference paper at ICLR 2025
Ilya O Tolstikhin, Sylvain Gelly, Olivier Bousquet, Carl-Johann Simon-Gabriel, and Bernhard Schölkopf. Adagan: Boosting generative models. Advances in neural information processing systems. 30, 2017.
Zhuowen Tu. Learning Generative Models via Discriminative Approaches. In 2007 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1-8, June 2007. doi: 10. 109/CVPR.2007. 383035. ISSN: 1063-6919.
David Watson, Kristin Blesch, Jan Kapar, and Marvin N. Wright. Adversarial random forests for density estimation and generative modeling. In Francisco Ruiz, Jennifer Dy, and Jan- lan-Willem van Meent (eds.), Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pp. 5357-5375. PMLR, 25-27 Apr 2023.
Max Welling, Richard Zemel, and Geoffrey E Hinton. Self Supervised Boosting. In Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002.
Hongwei Wen and Hanyuan Hang. Random forest density estimation. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 23701-23722. PMLR, 17-23 Jul 2022.
Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni. Modeling Tabular data using Conditional GAN. In Advances in Neural Information Processing Systems, volume 32.
Curran Associates, Inc., 2019.
Zilong Zhao, Aditya Kunar, Robert Birke, and Lydia Y Chen. Ctab-gan: Effective table data synthesizing. In Asian Conference on Machine Learning, pp. 97-112. PMLR, 2021.
13

---

Published as a conference paper at ICLR 2025
Appendix
Theory
14
A.1 Application to Piecewise Constant Functions
16
A.2 A Simple Convergence Result
17
A.3 NRGBoost Algorithm....
17
Density Estimation Trees and Density Estimation Forests
17
C Greedy Tree-Based Multiplicative Boosting
19
D Reproducibility
20
D.I Datasets
20
D.2 DEF and NRGBoost Implementation Details
21
D.3 Hyperparameter Tuning
21
D.4 Evaluation Setup
24
D.5 Computational Resources
25
Additional Results
25
E.1
Statistical Distance Measure
25
E,2
Uncertainty Estimation Regression Tasks
26
E.3 Inference with Missing Data
27
E,4 Other Challenging Tasks
27
E.5 MCMC Chain Convergence on MNIST
28
E,6
Computational Effort
28
A THEORY
The expected log-likelihood for energy-based model (EBM),
(x))) exp (f(x))
(13)
z(f]
given by
[[]] = Ex-p log yf (x))xxxx))x)))))>)) = (x) - loog Z[
(14)
The first variation of L can be computed as
gL[f:g] dLff + egl
Exxp g(x) 8 log Zif:gi Exxp g(x) Exway g(x). (15)
c=0
This is linear functional of its second argument, g. and can be regarded as directional derivative of I. at along variation g. The last equality comes from the following computation of the first
14

---

Published as a conference paper at ICLR 2025
variation of the log- partition function: ó log Zlf:gl
(16)
Z[f]
= Z[fl
1
Σ exp' (f(x)) 9(x)
(17)
exp (x())
g(x)
(18)
Zifl
= Ex Ex 9(x)
(19)
Analogous to a Hessian, we can differentiate Equation 15 again along a second independent variation of f. yielding symmetric bilinear functional which we will write as L g. h.. Note that the first term in Equation 14 is linear in f and thus has no curvature, so we only have to consider the log-partition function itself:
TT h) = L2 + eg + sh]
(20)
Ocde
(e.z)=0
8 - log Zif: g. h] = 88 8 (ólog Zif:g]] [f:h]
(21)
exp ((xx)))){} [s:h]
(22)
8Z1f; h]
Σ cxp (f(x)) g(xx - Zif] Σ exp' ((x))(x)
(23)
Z
87[[:h] h]
Ex g(x) z
Σ exp g())((x)
(24)
Z
=Exngy h(x) Bx g(x)
Ex h hxx)g(x)
(25)
-Co (g(x), h(x)).
(26)
Note that this functional is negative semi-definite for all f. In other words, we have 82 h] S 0 for any variation h, meaning that the log-likelihowe
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 + f + of up to second order öf:
1 = 8f] + 2
82 Lif:
(27)
1
0f Ex 6f(x) Var (x))
(28)
5J[ of a functional J implicitly by As an aside, defining the functional derivative, 3
Σ 8J1F] - {J[f:g).
Eafls)o(x o(x)
(29)
5f(x)
we can formally define, by analogy with the parametric case, the Fisher information "matrix" at as the following bilinear functional of two independent variations and h:
82 log gr (x)
Flf:g. g. h]
Σ
Ex 8f(y)8[(z)
g(y)h(z)
(30)
log Z[.]]
8f(y)8f(z) g(y)h(z)
(31)
82 log Zlf: g. h) -
(32)
To summarize, we have that the Fisher information is the same as the negative Hessian of the log- likelihood for these models. This is essentially due to the fact that the only term in log 9 (x) that non-linear in ff is the log-partition, which is not a function of x, and therefore it doesn't 't matter whether the expectation is taken under p or gf-
15

---

Published as a conference paper at ICLR 2025
A.1 APPLICATION TO PIECEWISE CONSTANT FUNCTIONS
Considering weak learner such
8ʃ(x) Σ wyyxx((x;
(33)
where the subsets Xj are disjoint and cover the entire input space, X, we have that
Ex
((x)
Σ wylx, (x)
(34)
XEX
J
g(x) = Σ Qx)))
(35)
j=1
xEX;
1=1
Similarly, making use of the fact that 1x, (x) (x) dig xx, (x), we can compute
of
E ((x)
Euj1x, (x))
Ew)e(x)).
(36)
xEX
In fact, we can extend this to any ordinary function of 8f:
Ex~g 9 (8f(x))
Σ (x)
Σ
1x, (8f(x))
(37)
XEX
j=1
Σ ((x)g(1;)
(38)
flxEx,
gg(w)))((,))
(39)
where we made of the fact that the 1x, constitute partition of unity:
Σ 1x, (x)
Vx
(40)
Finally, we can compute the increase in log-likelihood from a step f - f + a. of as
Llʃ + 8ʃ L Exxp (a 5[())] log Z ō.] + log Z
(41)
arl 0f - log Ex exp(co
(42)
wy P (X,) log Σ Qr ex
exp (aww)
(43)
where in Equation 42 made use of the equality:
log Z + a. 8ʃ - log Z[ʃ] = log Ex exp(f(x) + ad axf()))
log Σ (x) exp((6) exx()of(xx),
2[5]
(44)
and of the result in Equation 39 in the final step. This result can be used conduct line search over the step-size using training data and to estimate an increase in likeliioddoood at each round of boosting, for the purpose of early stopping, using validation data.
16

---

Published as a conference paper at ICLR 2025
A.2 A SIMPLE CONVERGENCE RESULT
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, wi ((x))/2((x,) - 1, in Equation 43 we get:
pa (Xz) Qx(X,)
Qr (X,) exp [ ( (((x))
1 (Q,(X,)
-1)
LV + a- -00 - -M -a
log Σ j=1
(45)
=xx(( log exp (aW):
(46)
where, with some abuse of notation, we denote by (P)()) the diverregee between two discrete distributions over [L..J] induced by the partitioning obtained in the current boosting round. This value is always non-negative, and can only be zero if P (X,) (X,)) for all leaves. It therefore corresponds to an increase in likelihood that as large as one can make the difference between the induced P and Q distributions by a judicious choice of partitioning, X,. Note that this quantity is precisely the objective that NRGBoost tries to greedily maximize (see Equation 9).
The second term in Equation 46 can be interpreted as the log of a moment generating function for a centered random variable W that takes values in {P(X))/Q, - }) (i.e., the set of the leaf values), with probabilities (Q)(())))E(1..)) With our proposed regularization approach, we limit the ratio P(X;)/Q)(X,) in choice of X5 maximum value R and therefore we have that W E [-1, R 11. Given that W bounded, by Hoeffding's Lemma (1963), sub-gaus uub-gaussian with variance proxy R²/ /4 and the log of its moment generating function is thus upper bounded by:
Ln a R
log E exp (aW)=
(47)
As result, we have the following lower bound on the likelhoodd increase:
1i22
L + a 8ʃ] L[Л] 2 a (P(())))) - 5
(48)
As long as, the current round of boosting, a partitioning of the input space can be found that yields a non-zero (P the likelihood is guaranteed to increase as long as we choose small enough step size, < R (PIIQ)) In particular, choosing step size = Rे (PI(Q)) produces an increase in likelihood of at least 2 [=)
We note that this result, while insightful, assumes that 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. We leave further analysis for future work.
A.3 NRGBOOST ALGORITHM
In Algorithm 1 we provide high-level overview of the training loop for NRGBoost.
B DENSITY ESTIMATION TREES AND DENSITY ESTIMATION FORESTS
Density Estimation Trees (DET), introduced by Ram & Gray (2011), model the density function as piecewise constant function,
q(x) = Σ By (x),,
(49)
where the X, are given by partitioning of the input space, X, induced by binary tree and the uy are the density values associated with each leaf that, for the time being, we will only require to be such that q(x) sums to one. Note that it is possible draw exact sample from this type of model
17

---

Published as a conference paper at ICLR 2025
Algorithm 1 NRGBoost training
Input: Empirical training distribution p. initial model distribution go. size of the sample pool M,
fraction of samples to independently drop at each round Prefresh E [0, 1], number of boosting rounds T, maximum number of leaves in each tree J, maximum probability ratio in each leaf R > shrinkage parameter E ]0,
Output: Energy function f.
f < log qo
for t f- T do
if itt= = then
B No sampling needed since is normalized
g + go
else if 2 then
> Initialize sample pool
ÿ + (x.he[M) and exp
By rejection sampling using do as a proposal distribution
else
I Drop samples from 7 randomly with probability Prefresh
I < Keep each sample, x,, from ÿ with probability given by Equation 12 using öft-1
④ + Add M samples to q by Gibbs sampling from exp f
end if
f t FitTree J. J.
Greedily grow tree using Equation 10 as a splitting criterion Line search training likelihood (Equation 43)
04 f- arg max, AL 5ft
+ + 7 a: 6.e
end for
by randomly selecting leaf, j E [1..J], given probabilities vj. and then drawing a sample from uniform distribution over X,.
To fit DET, Ram & Gray (2011) 1) proposes optimizing the Integrated Squared Error (ISE) between data generating distribution, p(x) and model:
min Σ (p(x) - qxx))
(50)
XEX
Noting that U/=11 X, X' we can rewrite this as
Σ
(u) - 20yp(x))
min
XEX
j=1 xex,
(51)
ΣΣ 11
s.t.
x=
Since the first term in the objective does not depend on the model, this optimization problem can be further simplified as
min U1 2. X.
-2,
X
j=1
(52)
s.t.
e (X,) 11
where V(X) denotes the volume of subset X. Solving this quadratic program for the oy we obtain
the following optimal leaf values and objective:
P(X,) V(X,)*
p²(X,)
ISE (X.,..,,,))
Σ
(53)
F=1 V)(X))
One can therefore grow a tree by greedily choosing to split parent leaf, with support Xp, into two leaves, with supports X2 and XR. so to maximize the following criterion:
(2T) + P²(XR)
P²(Xp)
max (EP(() V(X.) V(XR) V(Xp)
(54)
18

---

Published as a conference paper at ICLR 2025
This leads to similar splitting criterion to Equation 10, but replacing the previous model's distribution with the volume measure V, which can be interpreted as the uniform distribution on X' (up to a multiplicative constant).
Maximum Likelihood Often, generative models are trained to maximize the log-likelihood of the observed data,
max Ex~p log q(x).
(55)
rather than the ISE. This was left for future work in Ram & Gray (2011) but, following similar approach to the above, the optimization problem to solve is:
max
P(X,) log bj
#1
(56)
s.t. Σ 1=1
This is, again, easy to solve for Uj since it is separable over after removing the constraint using Lagrange multipliers. The optimal leaf values and objective are, in this case:
P(X,) 0j V(X,)*
Σ P(X,) log P(X;) j=1
L' 1X X,)
(57)
VX())
The only change is, therefore, to the splitting criterion which should become:
P(Xx) log P(X1.) + P(Xn) log v(X2)
P(Xn) V(Xa)
P(X)) log P(Xp) V(Xp)
max XEE(P))(XP)
(58)
This choice of splitting criterion can be seen as analogous to the choice between Gini impurity and Shannon entropy in the computation the information gain in decision trees.
Bagging and Feature Subsampling Following the common approach in decision trees, Ram & Gray (2011) 1) suggests the use of pruning for regularization of DET models. Practice has, however, evolved prefer bagging as a form of regularization rather than relying on single decision trees. We employ the same principle for DETs, by fitting many trees on bootstrap samples of the data. We also adopt the common practice from Random Forests of randomly sampling subset of features to consider when splitting any leaf node, in order to encourage independence between the different trees the ensemble. The ensemble model, which we call Density Estimation Forests (DEF), is thus an additive mixture of DETs with uniform weights, therefore still allowing for normalized density computation and exact sampling.
C GREEDY TREE-BASHR MULTIPLICATIVE BOOSTING
In multiplicative generative boosting, an unnormalized current density model, (x), updated at each boosting round by multiplication with new factor dq6" (x):
d (x) = (x) 5q, (x).
(59)
For our proposed NRGBoost, this factor chosen in order to maximize a local quadratic approx- imation of the log-likelihood around 4t-1, 1. functional of the log-density (see Section 3). The motivation behind the greedy approach of Tu (2007) and Grover Ermon (2018)
instead make the update factor ōq. (x) proportional to the likelihood ratio rx( = p(x)/9. i(x) directly. Under ideal conditions, this would mean that the method converges immediately when choosing a step size a₂ = 1. In more realistic setting, however, this method has been shown to converge under conditions on the performance of the individual oq: as discriminators between real and generated data (Tu, 2007; Grover & Ermon, 2018; Cranko & Nock, 2019).
In principle, this desired r, (x) could be derived from any binary classifier that is trained to predict probability of datapoint being generated (e.g., by training it to minimize strietly proper loss).
19

---

Published as a conference paper at ICLR 2025
Table 4: Comparison splitting criterion and leaf weights for the different versions of boosting.
Splitting Criterion
Leaf Values (Density)
Greedy (KL)
Plog (P/2) pa Q
P/Q
Greedy (x2)
P/Q exp (P/Q - 1)
NRGBoost
1
Grover & Ermon (2018), however, proposes relying on the following variational bound of an f- divergence to derive an estimator for this ratio:
(y(PIQr11) 2 sup [Exxp u(x) ff
(60)
well.
where f" denotes the convex conjugate of f. This bound is tight, with the optimum being achieved for un (x) = (x(xxx//=-1x) 14, contains this function. (f)) (uf (x)) can thus be interpreted as an approximation of the desired rx(x).
Adapting this method 60 to defined by tree functions u = 1/1 Σj= wi Ix,, with leaf values wy and leaf supports Xj- At each boosting iteration, a new tree, uf, can be grown to greedily optimize the lower bound in the r.h.s. of Equation 60, and setting (x) = (f')) (ut (x)). This means that öq: is also given by a tree with the same leaf supports and leaf values ey = (f') (wy)). leading to the following separable optimization problem:
trees as weak learners can be accomplished by considering Ue in Equation
max
Σ ((Xx)),-Q((x)FFF)
(61)
Note that we drop the iteration indices from this point onward for brevity. Maximizing over wy with the X, fixed, we have that w = I' (P(X))/Qxx)). which yields the optimal value
J' X,) Σ [Px,)) (((X)) Q(x,)()) f')
((X,)
(62)
Q(1,)
Q(2,)
This, in turn, determines the splitting criterion as a function of the choice of f. Finally, the optimal density values for the leaves are given by
P(X,)
- = ( (w Q(x,)*
(63)
is interesting to note two particular choices of f-divergences. For the KL divergence, f(t) = flog and ft()) 1 + logt = (f*) (t), which leads to
P(X,) JKL(N X,) = = P(X,) log Q(X,)
(64)
as the splitting criterion. The Pearson x2 divergence, with f(t) = (t - 1)² 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 (P(X:)/Q(xx)) instead of the ratio directly. Table 4
summarizes these results.
Another interesting observation is that DET model can be interpreted as 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, x² divergence and logelikelihood to the choice of KL divergence.
D REPRODUCIBILITY
D.1 DATASETS
We use 5 datasets from the UCI Machine Learning Repository (Dheeru & Karra Taniskidou, 2017): Abalone, Physicochemical Properties of Protein Tertiary Structure (from hereon referred to as
20

---

Published as a conference paper at ICLR 2025
Protein), Adult, MiniBooNE and Covertype. We also use the California Housing dataset which was downloaded through the Scikit-Learn package Pedregosa et al. (2011) and a downsampled version of the MNIST dataset Deng (2012). Table 5 summarizes the main details of these datasets well the number samples used for train/validation/test for the first cross-validation fold.
Table 5: Dataset Information. We respect the original test sets of each dataset when provided, otherwise we set aside 20% of the original dataset as test set. 20% of the remaining data set aside validation used for hyperparameter tuning.
Abbr
Name
Train Val
Test
Num
Cat
Target
Cardinality
AB
Abalone
3342
835
7
Num
29
CH
California Housing
16512
4128
8
0
Num
Continuous
PR
Protein
36584
9146
9
Num
Continuous
AD
Adult
32560
16280*
6
8
Cat
2
MBNE
MiniBooNE
104051
26013
50
Cat
2
MNIST
MNIST (downsampled)
60000
10000*
196
Cat
10
CT
Covertype
464810
116202
10
Cat
7
Original test set was respected.
The toy data distribution used in Figure consists of a mixture of cight isotropic Gaussians with standard deviation of placed around circle with radius of 8. To obtain discrete input domain, we discretize the input space in 100 equally spaced bins between -1 and 11.
D.2 DEF AND NRGBOOST IMPLEMENTATION DETAILS
Discretization In practical implementation of tree based methods, first discretize the input space by binning continuous numerical variables by quantiles. Furthermore, we also bin discrete numerical variables in order keep their cardinalities smaller than 256. This also be interpreted as establishing a priori a set discrete values to consider when splitting on each numerical variable and is done for computational efficiency, being inspired by LightGBM (Ke al., 2017).
Categorical Splitting For splitting on 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 leaf according pre-defined sorting function and then choose the optimal many-vs-many split as if the variable numerical. The function used sort the values is the leaf value function. For splitting on categorical variable, Ti, we order each possible categorical value by P(a, k. (/ where X_; denotes the leaf support over the remaining variables.
Tree Growth Strategy We always grow trees in best-first order, meaning that we 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 43. Because evaluating multiple possible step sizes, oa, 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 al., 2020) and Numba. We implement the tree evaluation and Gibbs sampling in C, making use of the PCG library (O (O'Neill, 2014) for random number generation.
D.3 HYPERPARAMETER TUNING
D.3.1
XGBoOST
To tune the hyperparameters XGBoost we use 100 trials of random search with the search space defined in Table 6.
21

---

Published as a conference paper at ICLR 2025
Table 6: XGBoost hyperparameter tuning search space, 8(0) denotes point mass distribution O.
Parameter
Distribution or Value
learning_rate
LogUniform ([10 ..00)
max max_ Leaves
Uniform ((16, 32, 64, 128, 256, 512, 1024))
min_child
whiggh weight
LogUniform ([10"1,
reg Lan_da
0.5 8(0) + 0.5 LogUniform ([10-5 10])
reg_alpha
0.5 8()) + 0.5 LogUniform ([10 10])
max_leaves
0 (we already limit the number of leaves)
grew_policy
lossguide hist
tree_method
Each model is trained for 1000 boosting rounds on regression and binary classification tasks. For multi-class classification tasks maximum number of 200 rounds of boosting used due to the larger size the datasets and because a separate tree built at every round for each class. The best model was selected based 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 that boosting round on the test
D.3.2 NGBOOST
To tune the hyperparameters NGBoost we use 100 trials of random search with the search space defined in Table 7. We use maximum of 500 rounds of boosting and early stopping with patience of 50 rounds. For consistency with the hyperparameter tuning setup for the other models, we select the model that achieves the best R² score validation.
Table 7: NGBoost hyperparameter tuning search space.
Parameter
Distribution
learning rate
LogUniform ([10-9..0) 1.0))
dex_h min rrin_sampee____ leaf
Uniform ([2..00))
Uniform ((...0))
D.3.3 RFDE
We implement the RFDE method (Wen & Hang, 2022) after quantile discretization of the dataset and therefore split at the midpoint of the discretized dimension instead of the original one. When leaf support has odd cardinality over the splitting dimension, 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 random one-vs-all split.
For RFDE models train total of 1000 trees. The only hyperparameter that tune is the maximum number of leaves per tree, for which we test the values [2 27...., 214]. For the Adult dataset, due limitations of our tree evaluation implementation, only test values up to 213
D.3.4 DEF
We train ensembles with 1000 DET models. Only three hyperparameters are tuned, using three nested loops, with each loop running over the possible values of single parameter in dre-defined order. These are, in order of outermost to innermost:
The maximum number leaves per tree, for which we test the values 221 22, 20 22].
For the Adult dataset 213 instead of 214
22

---

Published as a conference paper at ICLR 2025
2. The fraction of features to consider when splitting node.
We test the values
[d-1/2, 10-1//,d11/,,1 1] with d being the dimension of the dataset.
3. The minimum number of data points per leaf. We test the values 30,
To speed up the hyperparameter tuning process, we exit each loop early when the best value found for the objective does not improve between consecutive iterations. For example, if for 2 leaves per tree, and feature fraction of d 1/2, the objective does not improve when going from minimum points per leaf to 10, we move on to the next value of the feature fraction. Similarly, if the best objective found for feature fraction d114 (when considering possible values of minimum number of points per leaf) does not improve over that for d 1/2, we move on to the next value of number of leaves.
D.3.5 ARF
For ARF we used the official python implementation of the algorithm (Blesch Wright, 2023). Due to memory (and time) concerns we train models with maximum of 100 trees for up to 10 adversarial iterations. We use a similar grid-search setup to DEF (with early stopping), tuning only the following two parameters:
1. The maximum number of leaves per tree, for which we test the values [2 2 ..... 2 oo] (oc denotes unconstrained).
2. The minimum number of examples per leaf, for which we test the values [3, 5, 10, 30, 50, 100).
D.3.6 NRGBOOST
We train NRGBoost models for maximum of 200 rounds of boosting. We tune only two parameters, using a similar early stopping setup in the outer loop:
1. The maximum number leaves, for which we try the values 22,,2,, 2,,200 2,0,212] in order. For the CT dataset also include 2 in the values to test.
2. The constant factor by which the optimal step size determined by the line search is shrunk at each round of boosting. This 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 for this search is [0.01, 0.5].
For regularization we limit the ratio between empirical data density and model data density on each leaf to maximum of 2. We observed that smaller values of this regularization parameter tend to perform better, but that it otherwise plays a similar role to the shrinkage factor. Therefore, we choose to tune only the latter.
For sampling, we use a sample pool of 80,000 samples for all datasets except for Covertype where use 320,000. These values are chosen so that the number of samples in the pool are, 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 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. However, we observed that this mixture coefficient does not have much impact on the results.
D.3.7 TVAE
We use the implementation of TVAE from the SDV package. 2 To tune its hyperparameters we use 50 trials of random search with the search spaces defined in Table 8.
https:://po thuithhub.ae com/ sdv dev/-dde/svv
23

---

Published as a conference paper at ICLR 2025
Table TVAE hyperparameter tuning search space. We set both compress_ __dims and decompress _erress__msssss to have the number of layers specified by num_ sam_layers, with hidden_doe diim 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 significantly longer to run a single epoch.
Parameter
Datasets
Distribution or Value
epochs
small
Uniform (100..500))
large
Uniform ((00.200))
batch_size
small
Uniform (f100, 200, ..0,,...:00000 500))
large Uniform (1500, 1000, .5000)
embedding_dim
all
Uniform (132, 64, 128, 256, 512))
hidden_ _im
all Uniform ((32, 64, 128, 256, 512))
num_layers
all Uniform 2,3))
compress _dims
all (hidden cim,)
layers
decompresasdiim
all (hidden_dim,) num_layers
D.3.8 TABDDPM
We used the official implementation of TabDDPM, adapted datasets and validation setup. To tune the hyperparameters TabDDPM we use 50 trials of random search with the same search space used by Kotelnikov al. (2023).
D.3.9 FOREST-FLOW
We use the official implementation of Forest-Flow." The only hyperparameterrs with impact performance that the Forest-Flow documentation recommends tuning are the number diffusion timesteps (default of 50) and the number of noise values per sample (default of 100). Due to the fact that Forest-Floo models already take significantly longer to train than the other models, and that increasing either value would lead slower training, opted these default values.
D.4 EVALUATION SETUP
Single-variable inference For the single-variate inference evaluation, the best models selected by their discriminative performance on validation set. The entire setup 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,
single train/validation/eest split of the real data (corresponding to the first fold in the previous setup) for training the generative models. The density models used those previously selected based on their single-variable inference performance on the validation set. For the sampling models (TVAE and TabDDPM) we directly evaluate their ML efficiency using the validation data.
ML Efficiency For each selected model sample train and validation sets with the same number of samples as those used in 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 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:
research http:::://gtttbb.oo/oo resaree/a-eeeerrccccaabb--ee ttchhttbbddmm 'https: //gi http:://ggghhhb..cossseaeeeeett www..////w/sssaaa..oo I
ffusion
24

---

Published as a conference paper at ICLR 2025
The original validation set is used 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 generated 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 respectivelly. 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 not generate correlated data with NRGBoost. We draw each sample by rumning its own independent chain for 100 steps, starting from an independent sample from the initial model, which is rather slow process. The setup is repeated 5 times with 5 different sets of generated samples from each method.
D.5 COMPUTATIONAL RESOURCES
The experiments were run on Linux machine equipped with an AMD Ryzen 7 7700X 8 core CPU and 32 GB of RAM. The comparisons with TVAE and TabDDPM additionally made use of a GeForce RTX 3060 GPU with 12 GB of VRAM.
E ADDITIONAL. RESULTS
E. 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 Maicteur-Martineau al. (2024). Namely, we min-max scale numerical variables, and one-hot encode categorical variables, scaling the latter by 1/2. We use a L1 distance in the formulation of the optimal transport problem, which we solve using the POT python library. Since finding the optimal solution scales cubically with the sample size in the worst case, we sub-sample at most 5,000 samples from the original train or test sets and 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). Our results for the distance between empirical training and test distribution and synthetic samples from each model are reported in Tables 9 and 10 respectively.
Table 9: 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 20.005
0.195
0.274 20.011
1.075 20/024
0.366 #0.055
20.59 20:01
0.808
TabDDPM
0.467 20,021
0.153 20.007
0.190 20.005
0.895 20,019
109.6 493.14
0.640 20:016
Forest-Flow
0.234 20,020
0.192 20,007
0,238 20.009
1.239 a0.024
4.214 23.683
20.73 20.574
0.896 20.022
ARF
0.199 10.015
0.233 10.009
0.232 20.007
1.363 10.018
0.698
10.684
18. 14 20.060
0.703 20.019
DEF (KL)
0.320 #0.006
0.199 #0,012
0,260 40.010
1.682 #0.024
18.21 +14.92
261.0 #15.60
0,933 #0.025
NRGBoost
0.509 40.047
0.170 20.012
0.230 40.011
1.028 #0,015
22.25 #20.09
15.60 #0.212 0.600 20.016
We found these results to be somewhat sensitive to the choice of normalization of the numericals and
categoricals, as well as 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.
25

---

Published as a conference paper at ICLR 2025
Table 10: 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.025
0.425 20,032
20.61 +0.042
0.807 +0.023
TabDDPM
0,500 +0,023
0.169 20,009
0.195 s0.00
0.938
s0.029
109.5 193.13
0,646 +0.022
Forest-Flow
0.248 10.020
0.205 s0.009
0.231 s0.008
1.243 s0.029
4.295 x3,697
20.73 10:574
0.893 s0.021
ARF
0.218 20.015
0.240 20.009
0.238 20.006
1.378 #0.020
0.732 #0,664
18.14 20.060
0.709 20.018
DEF (KL)
0.337 10.006
0.210 20.012
0.264 20.011
1.696 a0.023
18.29 114.90
261.0 101.0.m +35.01
0,935 20.024
NRGBoost
0.502 10.045
0.178 20.011
0.225 20.010
1.051 10.013
22.33 120.11
15.70 10.188
0.602 20.015
Table 11: (C)RPS and MAE results. The reported values are means and standard errors over cross-validation folds. The best method for each dataset is highlighted in bold and other methods that are not significantly worse (as determined by a paired t-test 95% confidence level) are underlined.
AB
CH
PR
RPS
MAE
CRPS +
MAE +
CRPS J
MAE J
XGBoost
N/A
1.524 50032
N/A
0.290 20,009
N/A
2.282 20.045
NGBoost
1.094 20.000
1.525 20012
0.238 20.013
0.314 20.01
1.976 20018
2.665 =0057
NRGBoost
1.075 s0.01
1.569 20:025
0.201 20.008
0.276 20.010
1.490 20.023
2.063 =0,036
E.2 UNCERTAINTY ESTIMATION IN REGRESSION TASKS
One advantage of density models over discriminative regression models like XGBoost is that they estimate the full conditional distribution, p(y)x x)) rather than provide a point estimate. Section 5.1, we focused on a metric that evaluates the point estimate E[y|x]. To further evaluate the quality of the models' probabilistic predictions, in Table 11 we report two other metrics:
We compute the average Mean Absolute Error given a point estimate of y. For NGBoost and NRGBoost, we use the median of the estimated q(ylx) since this is the point estimate that minimizes MAE. Note that, while it is possible to train XGBoost to directly optimize this metric instead of MSE, models estimating (z)(x) are flexible choose different point estimate post-training.
MAE:
The Continuous Ranked Probability Score a proper scoring rule commonly used to evaluate a probabilistic forecast (Gneiting & Raftery, 2007). For both NGBoost and NRGBoost, we evaluate this score analytically for the forecasted q xx) of every test point. Because the Abalone dataset has discrete numerical targets, compute RPS for this dataset instead. For this purpose, the continuous distribution obtained with NGBoost first discretized by rounding the nearest discrete value within the domain of the target.
(C)RPS:
Note that we omit the results for ARF from this analysis because the Python implementation does not provide neither conditional nor unconditional density evaluation at the time of writing. For the results reported in Table 1, we implemented our own evaluation of Efy] x for these models. Extending this implementation to compute the median or the CRPS analytically is, however, not trivial due to ARF fitting more complex continuous distributions at the leaves than the other tree- models.
On MAE, we find that NRGBoost outperforms XGBoost convincingly in the California Housing and Protein datasets, in part due to the use of a better point estimate for this metric. The Abalone dataset is the only dataset where predicting the mean of the forecasting distribution would outperform the median for MAE. We also find that NRGBoost achieves a lower (C)SPS than NGBoost on all datasets, which attribute the parametric choice made for NGBoost not necessarily being the most well-suited for these datasets. This highlights a key drawback of parametric approaches compared to maramereit ones: they require domain expertise to avoid suboptimal results.
26

---

Published as a conference paper at ICLR 2025
Table 12: Results for inference with missing covariate. The original results of each method when having access to the full data are also reported for reference. The best approach for dealing with missing data in each scenario is highlighted in bold.
Model
Imputation
CH (A² 1)
AD (AUC1)
CT (Accuracy 1)
XGBoost
Full Data
0.849 10.009
0.927 20.000
0.971 20.001
Mean
-0.283
N/A
0.610 #0.004
Median/Mode
-0.117 10.107
0.914 200003
0.621 40.002
KNN (K=5)
0.150 20.107
0.910 20003
0.883 20.001
NRGBoost
Full Data
0.850 10.011
0.920 s0.001
0.948 +0.001
Marginalization
0.773 10.010
0.920 .0.001
0.923 a0.001
E.3 INFERENCE WITH MISSING DATA
Generative density models can be used more flexibly for inference than their discriminative counter- parts. In this section we present a case study on a problem of inference in the presence of missing covariate that highlights this advantage.
For each type of task, we choose dataset and remove an important feature from the respective test set (Longitude for the California dataset, education for the Adult dataset and Elevation for the Covertype dataset). We then attempt to predict the value of the target variable using XGBoost with the help of several imputation techniques that are popular with data science practitioners:
Mean and median imputation (for a missing numerical feature).
Mode imputation (for missing categorical feature).
Imputation with the mean or mode of the 5 nearest neighbors on the training data. Note that this imputation requires one to have access to the entire training set at test time.
NRGBoost can handle this type of problem in principled manner by marginalizing over the possible values of the missing feature, rather than relying on ad-hoc imputation. Our results, in Table 12, show that this approach is able to outperform imputation, even on datasets where performance gap between NRGBoost and XGBoost existed on the full data.
E.4 OTHER CHALLENGING TASKS
In order to test the limits of NRGBoost, we compare its performance to XGBoost on the following two challenging tasks taken from Gorishniy et al. (2021):
ALOI:
This is an image classification dataset with 128 features and 100k examples dis- tributed evenly over 1000 classes. Because implementation of NRGBoost is limited to variables with a cardinality of 255, we consider only data for the first 250 classes. This still constitutes 25 times more classes than MNIST, presenting challenge to both discriminative and generative models.
This regression task contains both a large number of samples (1.2M) and features (136). making it one order of magnitude larger than MNIST and Covertype in terms of the product of the two. Because of this, we found that our usual hyperparameter tuning protocol for NRGBoost to be impractical and had resort to manual tuning instead. The best model we found had 256 leaves and was trained for 500 boosting iterations with a maximum ratio in each leaf of 4 and a shrinkage factor of 0.
Microsoft:
We use the same train, validation and test splits used Gorishniy et al. (2021) and report the same metrics in Table 13. We find that NRGBoost outperforms XGBoost in the ALOI dataset but that there still a significant gap between the best NRGBoost model that we were able to find on the Microsoft dataset and XGBoost.
27

---

Published as a conference paper at ICLR 2025
Table 13: Comparison of XGBoost and NRGBoost on ALOI and Microsoft datasets.
ALOI (Accuracy f)
Microsoft (RMSE
XGBoost NRGBoost
0.943
0.744
0.962
0.787
7
7
1
G
G
3
T
2
5
5
&
8
8
3
3
3
3
3
3
3
6
6
6
2
2
2.
2
2
2
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
2
3
7
3
3
7
Figure 5: Downsampled MNIST samples generated by Gibbs sampling from an NRGBoost model. Each row corresponds to an independent chain initialized with a sample from the initial model fo (first column). Each column represents a consecutive sample from the chain.
E.5 MCMC CHAIN CONVERGENCE ON MNIST
Figure 5 shows the convergence of Gibbs sampling from an NRGBoost model. In only few samples each chain appears have converged to the data manifold after starting at 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 a chain to switch between two different modes of the distribution (e.g., switching digits) even though a few such transitions can be observed.
E.6 COMPUTATIONAL EFFORT
E. TRAINING
In Figure 6 we report the training times for NRGBoost, 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 the other methods given that the splitting process is random and does not depend on the data, only on the input domain.
Note that the biggest computational cost for training 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 (16 virtual cores).
We note also that we believe that there is still plenty of margin for optimizing the tree-fitee soe-fitting code used for DEF models. As such, the results presented here are merely indicative.
The computational effort required for fitting a tree-based model should scale linearly with both the number of samples and the number of features. This is the main source of variation for the training times across datasets. However, we note that larger datasets can benefit more from larger models (c.g., with a larger number of leaves) which also take a longer time train.
E.6.2 SAMPLING
In Table 14, we report the time required to draw 10,000 independent samples from each model on each dataset. We do not report sampling times for DEF models because current implementation
28

---

Published as a conference paper at ICLR 2025
Smaller Datasets
Larger Datasets
12
120
TVAE
10
TahDDPM
100
ARF
DEP (ISE)
Wall Time
DEP (KL)
NRGBoost
Wall
abalone
california
adhll
miniboone
mnist dataset
covertype
dataset
Figure 6: Wall time required to train the best model for each method.
Table 14: Time required draw 10000 independent samples in seconds. The last columns shows how much faster method is compared to NRGBoost in the median case. We find that TabDDPM, the only method that realistically competes with NRGBoost in sample quality is, in the median case, 10x faster for sampling.
AB
CH
PR
AD
MBNE
MNIST
CT
Faster (Median)
TVAE
0.06
0.05
0.04
0.05
0.22
1.33
0.07
984.0
TabDDPM
42.3
7.36
6.82
1.68
1.13
9.69
9.67
ARF
0.38
0.22
0.40
2.62
3.84
4.47
55.7
NRGBoost
23.1
49.2
54.6
65.0
194.2
622.4
109.9
of this procedure is very suboptimal, taking longer than sampling from NRGBoost. We helieve, however, that more optimized approach would, at most, be somewhat slower than sampling from ARF models, since the DEF ensembles larger.
For NRGBoost, independent samples are generated by starting from an independent sample from go and then running fixed number, of Gibbs sampling iterations, taking only the last iterate sample. However, this process only produces an approximate sample from the energy model, as it remains biased by the initial sample, which is drawn from an incorrect distribution. The larger B is. compared to the Markov chain's 's mixing time, the more independent the final sample is from the initial one and the smaller this bias. This parameter, therefore, trades-off correctness for computational time, which scales linearly with B.
our experiments we always set = 100, this may excessive for some datasets. We find that B = 30 already leads to similar ML efficiency results on most datasets, while cutting down sampling times by a factor of three. Note also that the process of drawing samples is embarrassingly parallel. The numbers we were to leverage a CPU with higher core count, we could similarly cut down on sampling times. In fact, the process of sampling from EBM is not dissimilar the process of sampling from diffusion model like TabDDPM, but the latter can leverage the massive parallelism afforded by GPUs.
report correspond running 16 chains in parallel on a CPU with 16 virtual cores. If
29


Reviews:
[{"review_id": "DIily5lQM5", "reviewer": "Reviewer_cBTV", "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.", "questions": "I would like to increase my score according to the comments from other reviewers and the response.\n\n(1) What about the inference time comparison between NRGBoost and other baselines?\n\n(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?\n\n(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])?\n\n**Reference**\n\n[1] Revisiting Deep Learning Models for Tabular Data, NeruIPS 2021.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 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.\n\n**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.\n\n**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.\n\n**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.", "weaknesses": "Several weaknesses from application perspective.\n\n(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.\n\n(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.\n\n(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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "85jwmNOJKZ", "reviewer": "Reviewer_pEfp", "summary": "This manuscript presents a framework for using boosting to model unnormalised densities.\nThe proposed framework is used to build a model using decision trees as weak learners.\nTo accelerate sampling from the model (and thus training), a combination of Gibbs and rejection sampling is used.\nExperiments on a down-sampled MNIST and UCI datasets demonstrate competitive performance to other tree-based generative models.", "questions": "1. How does the performance compare to DL-based generative models?\n2. Can the tradeoffs of NRGBoost vs e.g. DET/DEF be quantified somehow?\n3. Is there a concrete use case where NRGBoost is the only/best model?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- (clarity) Very well-written paper that is easy to follow.\n   I also appreciate the clear presentation and discussion of limitations of the proposed approach.\n - (originality) The proposed boosting framework seems to be novel and provides some refreshing insights.\n - (originality/quality) The overview over existing approaches to use decision trees for generative modelling is extensive and contextualises the work well.\n - (quality) The derivations are explained well and the appendix provides useful additional details.\n - (quality) The experiments seem to have been set up properly with error bars and competitive baselines.", "weaknesses": "- (significance) The abstract mentions a comparison with neural network models, but I could not find these in the main paper.\n   Furthermore, the authors claim that tree-based generative models are to be preferred over DL-based models\n   because discriminative DL methods do not provide competitive performance on tabular data.\n   However, no DL methods have been included as baselines in the experiments to confirm this.\n   As a result there is no evidence that DL methods would not be able to outperform tree-based generative models.\n   PS: DL methods based on normalizing flows have explicit (tractable) densities.\n - (significance) It is not entirely clear how the proposed model can/should be used.\n   The discussion mentions the overhead due to sampling, but this overhead is never quantified.\n   This makes it hard to get a feeling for whether the (sometimes minor) increase in performance is \"worth it\".\n   Also, it is not entirely clear why shallower trees are preferrable, as implied in the discussion.\n - (significance) The authors put a lot of emphasis on the tractability of computing densities.\n   However, there are no clear experiments that illustrate the importance of this feature.\n   The conclusion hints at applications enabled by density models,\n   but this seems to indicate that there are no real use cases for this model.\n   Especially since the results seem to indicate that there are other competitive methods out there.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "RkTCu2pXyB", "reviewer": "Reviewer_a4Ct", "summary": "The paper proposes a tree-based approach to model tabular data density which they claim outperforms traditional generative approaches.", "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.\n\n2. 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. \n\n3. 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. \n\n4. 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.\n\n5. 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. \n\nOverall, I think the paper is far from being a coherent and legit study.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper is confusing and hard to follow. I was unable to find any strength after investing a long time.", "weaknesses": "The paper is poorly written, incoherent and uses inconsistent experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EHN1F1s6PW", "reviewer": "Reviewer_NED1", "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.", "questions": "Questions\n\n- 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.\n\n- I realize that NGBoost does not support nonparametric probabilistic prediction, but I think seeing its results in Table 1 would still be helpful. \n\nMinor suggestions for improved readability:\n\nReplace: We therefore estimate these quantities, with recourse to empirical data for P(X), and to samples approximately drawn from the model with MCMC.\nWith: We therefore estimate these quantities, with recourse to empirical data for P(X) and to samples approximately drawn from the model with MCMC.\n\nReplace: 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.\nWith: 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.", "limitations": "", "rating": 10, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 4, "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.\n\n2. 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.\n\n3. 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.", "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.\n\n2. [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. \n\n3. [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.\n  \n[1] Zheng, S., & Charoenphakdee, N. (2022). Diffusion models for missing value imputation in tabular data. arXiv preprint arXiv:2210.17128.\n\n[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.\n\n[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.\n\n[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.\n\n[5] McCarter, C. (2024). Unmasking Trees for Tabular Data. arXiv preprint arXiv:2407.05593.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)