NRGBoost: Energy-Based Generative Boosted Trees

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 09:21:05

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 effective energy-based generative boosting algorithm that is analogous to the
second order boosting algorithm implemented in popular packages like XGBoost.
We show that, despite producing a generative model capable of handling inference tasks over any input variable, 
our proposed algorithm can achieve similar discriminative performance to GBDT algorithms on a number of real world tabular datasets and outperform competing approaches for sampling.

Main paper: 
 introduction: 
 Generative models have achieved tremendous success in computer vision and natural language processing, where the ability to generate synthetic data guided by user prompts opens up many exciting possibilities. While generating synthetic table records does not necessarily enjoy the same wide appeal, this problem has still received considerable attention as a potential avenue for bypassing privacy concerns when sharing data. Estimating the data density, p(x), is another typical application of generative models which enables a host of different use cases that can be particularly interesting for tabular data. Unlike discriminative models which are trained to perform inference over a single target variable, density models can be used more flexibly for inference over different variables or for out of distribution detection. They can also handle inference with missing data in a principled way by marginalizing over unobserved variables.The development of generative models for tabular data has mirrored its progression in computer vision with many of its Deep Learning (DL) approaches being adapted to the tabular domain [Jordon et al., 2018 , Xu et al., 2019 , Engelmann and Lessmann, 2020 , Fan et al., 2020 , Zhao et al., 2021 , Kotelnikov et al., 2022] . Unfortunately, these methods are only useful for sampling as they either don't model the density explicitly or can't evaluate it due to untractable marginalization over high dimensional latent variable spaces. Furthermore, despite growing in popularity, DL has still failed to displace tree-based ensemble methods as the tool of choice for handling tabular discriminative tasks with gradient boosting still being found to outperform neural-network-based methods in many real world datasets [Grinsztajn et al., 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 and Hang, 2022 , Nock and Guillame-Bert, 2022 , Watson et al., 2023 , Nock and Guillame-Bert, 2023 , Jolicoeur-Martineau et al., 2023] , we find that direct extensions of Random Forests (RF) and Gradient Boosted Decision Tree (GBDT) are still missing. It is this gap that we try to address, seeking to keep the general algorithmic structure of these popular algorithms but replacing the optimization of their discriminative objective with a generative counterpart. Our main contributions in this regard are:• Proposing NRGBoost, a novel energy-based generative boosting model that, analogously to the boosting algorithms implemented in popular GBDT packages, is trained to maximize a local second order approximation to the likelihood at each boosting round.• Proposing an approximate sampling algorithm to speed up the training of any tree-based multiplicative generative boosting model.• Exploring the use of bagged ensembles of Density Estimation Trees (DET) [Ram and Gray, 2011] with feature subsampling as the generative counterpart to RF.The longstanding popularity of GBDT models in machine learning practice can, in part, be attributed to the strength of its empirical results and the efficiency of its existing implementations. We therefore focus on an experimental evaluation in real world datasets spanning a range of use cases, number of samples and features. We find that, on smaller datasets, our implementation of NRGBoost can be trained in a few minutes on a mid-range consumer CPU and achieve similar discriminative performance to a standard GBDT model while also being able to generate samples that are generally harder to distinguish from real data than state of the art neural-network-based models.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: Claims about proposal of novel methods are justified in Sections 3 and 4.Claims about empirical results are justified in Section 6 and Appendix G.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.• The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.• It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
energy based models: 
 An Energy-Based Model (EBM) parametrizes the logarithm of a probability density function directly (up to an unspecified normalizing constant):EQUATIONHere f (x) : X → R is a real function over the input domain. 1 We will avoid introducing any parametrization, instead treating the function f ∈ F(X ) lying in an appropriate function space over the input space as our model parameter directly. Z[f ] = x∈X exp (f (x)), known as the partition function, is then a functional of f giving us the necessary normalizing constant. This is the most flexible way one could represent a probability density function making essentially no compromises on its structure. The downside to this is that for most interesting choices of F, computing or estimating this normalizing constant is untractable which makes training these models difficult. Their unnormalized nature however does not prevent EBMs from being useful in a number of applications besides sampling. Performing inference over a small enough subset of variables requires only normalizing over the set of their possible values and for anomaly or out of distribution detection, knowledge of the normalizing constant is not necessary.One common way to train an energy-based model to approximate a data generating distribution, p(x), is to minimize the Kullback-Leibler divergence between p and q f , or equivalently, maximize the expected log likelihood functional:EQUATIONThis optimization is typically carried out by gradient descent over the parameters of f , but due to the untractability of the partition function, one must rely on Markov Chain Monte Carlo (MCMC)sampling to estimate the gradients [Song and Kingma, 2021] .Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We discuss limitations of our proposed method both in the section that introduces it (Section 3) as well as in the experiments (Section 6) and discussion (Section 7)sections.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
nrgboost: 
 Expanding the increase in log-likelihood in equation 2 due to a variation δf around an energy function f up to second order we haveEQUATIONThe δf that maximizes this quadratic approximation should thus have a large positive difference between the expected value under the data and under q f while having low variance under q f . We note that just like the original log-likelihood, this Taylor expansion is invariant to adding an overall constant to δf . This means that, in maximizing equation 3 we can consider only functions that have zero expectation under q f in which case we can simplify ∆L f [δf ] asEQUATIONWe thus formulate our boosting algorithm as modelling the data density with an additive energy function. At each boosting iteration we improve upon the current energy function f t by finding an optimal step δf * t that maximizes ∆L ft [δf ]EQUATIONwhere H t is an appropriate space of functions (satisfyingE x∼q f t δf (x) = 0 if equation 4 is used).The solution to this problem can be interpreted as a Newton step in the space of energy functions.Because for an energy-based model, the Fisher Information matrix with respect to the energy function and the hessian of the expected log-likelihood are the same, we can also interpret the solution to equation 5 as a natural gradient step (see the Appendix A). This approach is essentially analogous to the second order step implemented in modern discriminative gradient boosting libraries such asXGBoost [Chen and Guestrin, 2016] and LightGBM [Ke et al., 2017] and which can be traced back to Friedman et al. [2000] .In updating the current iterate, f t+1 = f t + α t • δf * t , we scale δf * t by an additional scalar step-size α t . This can be interpreted as a globalization strategy to account for the fact that the quadratic approximation in equation 3 is not necessarily valid over large steps in function space. A common strategy in nonlinear optimization would be to select α t via a line search based on the original log-likelihood. Common practice in discriminative boosting however is to interpret this step size as a regularization parameter and to select a fixed value in ]0, 1] with (more) smaller steps typically outperforming fewer larger ones when it comes to generalization. We choose to adopt a hybrid strategy, first selecting an optimal step size by line search and then shrinking it by a fixed factor. We find that this typically accelerates convergence allowing the algorithm to take comparatively larger steps that increase the likelihood in the initial phase of boosting. For a starting point, f 0 , we can choose the logarithm of any probability distribution over X as long as it is easy to evaluate. Sensible choices are a uniform distribution (i.e., f ≡ 0), the product of marginals for the training set, or any mixture distribution between these two.As a weak learner we will consider functions defined by trees over the input space. I.e., letting J j=1 X j = X be the partitioning of the input space induced by the leaves of a binary tree whose internal nodes represent a split along one dimension into two disjoint partitions, we take as H the set of functions such asEQUATIONwhere 1 X denotes the indicator function of a subset X and w j are values associated with eachleaf j ∈ [1..J].In a standard decision tree these values would typically encode an estimate of p(y|x ∈ X j ), with y being a special target variable that is never considered for splitting. In our generative approach they encode unconditional densities (or more precisely energies) over each leaf's support and every variable can be used for splitting. Note that our functions δf are thus parametrized by the values w j as well the structure of the tree and the variables and values for the split at each node which ultimately determine the X j . We omit these dependencies for brevity.Replacing the definition in equation 6 in our objective (equation 4) we get the following optimization problem to find the optimal decision tree: max w1,...,w J ,X1,...,X J J j=1EQUATIONwhere P (X j ) and Q f (X j ) denote the probability of the event x ∈ X j under the respective distribution and the constraint ensures that δf has zero expectation under q f . With respect to the leaf weights this is a quadratic program whose optimal solution and objective values are respectively given byEQUATION)Because carrying out the maximization of this optimal value over the tree structure that determines the X j is hard, we approximate its solution by greedily growing a tree that maximizes it when considering how to split each node individually. A parent leaf with support X P is thus split into 2 child leaves, with disjoint support, X L ∪ X R = X P , so as to maximize over all possible partitionings along a single dimension, P (X P ), the following objective:EQUATIONNote that when using parametric weak learners, computing a second order step would typically involve solving a linear system with a full Hessian. As we can see, this is not the case when the weak learners are decision trees where the optimal value to assign to a leaf j does not depend on any information from other leaves and, likewise, the optimal objective value is a sum of terms, each depending only on information from a single leaf. This would have not been the case had we tried to optimize the likelihood functional in Equation 2 directly instead of its quadratic approximation.To compute the leaf values in equation 8 and the splitting criterion in equation 9 we would have to know P (X) and be able to compute Q f (X) which is infeasible due to the untractable normalization constant. We therefore estimate these quantities, with recourse to empirical data for P (X), and to samples approximately drawn from the model with MCMC. Because even if the input space is not partially discrete, f is still discontinuous and constant almost everywhere we can't use gradient based samplers and therefore rely on Gibbs sampling instead. This only requires evaluating each f t along one dimension at a time, while keeping all others fixed which can be computed efficiently for a tree by traversing it only once. However, since at boosting iteration t our energy function is a sum of t trees, this computation scales linearly with the iteration number. This makes the overall time spent sampling quadratic in the number of iterations and thus precludes us from training models with a large number of trees.In order to reduce the burden associated with this sampling, which can dominate the runtime of training the model, we propose a new sampling approach that leverages the cumulative nature of boosting. The intuition behind this approach is that the set of samples used in the previous boosting round are (approximately) drawn from a distribution that is already close to the new model distribution.It could therefore be helpful to keep some of those samples, especially those that conform the best to the new model. Rejection sampling allows us to do just that. The boosting update in terms of the densities takes the following multiplicative form:EQUATIONHere, k is an unknown multiplicative constant and since δf t is given by a tree, we can easily bound the exponential factor by finding the leaf with the largest value. We can therefore use the previous model, q t-1 (x), as a proposal distribution for which we already have a set of samples and keep each sample, x, with an acceptance probability of:EQUATIONWe note that knowledge of the constant k t is not necessary to compute this acceptance probability.After removing samples from the pool, we can use Gibbs sampling to draw a new set of samples in order to keep a fixed total number of samples per round of boosting. Note also that q 0 is typically a simple model for which we can both directly evaluate the desired quantities (i.e., Q 0 (X) for a given partition X) and cheaply draw exact samples from. As such, no sampling is required for the first iteration of boosting and for the second we can draw exact samples from q 1 with rejection sampling using q 0 as a proposal distribution.This approach works better when either the range of f t is small or when the step sizes α t are small as this leads to larger acceptance probabilities. Note that in practice it can be helpful to independently refresh a fixed fraction samples, p ref resh , at each round of boosting in order to encourage more diverse samples between rounds. This can be accomplished by keeping each sample with a probabilityp accept (x)(1 -p ref resh ) instead.The simplest way to regularize a boosting model is to stop training when overfitting is detected by monitoring a suitable performance metric on a validation set. For NRGBoost this could be the increase in log-likelihood at each boosting round. However, estimating this quantity would require drawing additional validation samples from the model (see Appendix A). An alternative viable validation strategy which needs no additional samples is to simply monitor a discriminative performance metric (over one or more variables). This essentially amounts to monitoring the quality of q f (x i |x -i ) instead of the full q f (x).Besides early stopping, the decision trees themselves can be regularized by limiting the depth or total number of leaves of each tree. Additionally we can rely on other strategies such as disregarding splits that would result in a leaf with too little training data, P (X), model data, Q f (X), volume V (X) or too high of a ratio between training and model data P (X) /Q f (X). We found the latter to be the most effective of these, not only yielding better generalization performance than other approaches, but also having the added benefit of allowing us to lower bound the acceptance probability of our rejection sampling scheme.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [No] Justification: Unfortunately we did not have time to clean up the code and document it so that it could be useful at the time of the paper deadline. But we intend to make our implementations of the proposed algorithms available as a python library as soon as possibleand will also open source the full experimental setup on Github.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
density estimation trees and density estimation forests: 
 Density Estimation Trees (DET) were proposed by Ram and Gray [2011] as an alternative to histograms and kernel density estimation but have received little attention as generative models for sampling or other applications. They model the density function as a constant value over the support of each leaf in a binary tree, q = J j=1 P (Xj )V (Xj ) 1 Xj , with P (X) being an empirical estimate of probability of the event x ∈ X and V (X) denoting the volume of X. Note that it is possible to draw an exact sample from this type of model by randomly selecting a leaf, j ∈ [1..J], given probabilities P (X j ), and then drawing a sample from a uniform distribution over X j .To fit a DET, Ram and Gray [2011] propose optimizing the Integrated Squared Error (ISE) between the data and model distributions which, following a similar approach to Section 3.1, leads the following optimization problem when considering how to split a leaf node:EQUATION)For the ISE, D should be taken as the function D ISE (P, V ) = P 2 /V which leads to a similar splitting criterion to Equation 12 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 likelihood of the observed data. This was left for future work in Ram and Gray [2011] but, as we show in Appendix B, can be accomplished by replacing the D in Equation 12with D KL (P, V ) = P log ( P /V ).This choice of minimization criterion can be seen as analogous to the choice between Gini impurity and Shannon entropy in the computation of the information gain in decision trees.Bagging and Feature Subsampling Following the common approach in decision trees, Ram and Gray [2011] suggest the use of pruning for regularization of DET models. Practice has however evolved to prefer bagging as a form of regularization rather than relying on single decision trees. We employ same principle to DETs by fitting many trees on bootstrap samples of the data. We also adopt the common practice from Random Forests of randomly sampling a subset of features to consider when splitting any leaf node in order to encourage independence between the different trees in the ensemble. The ensemble model, which we call Density Estimation Forests (DEF) in the sequence, is thus an additive mixture of DETs with uniform weights, therefore still allowing for normalized density computation and exact sampling. 
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 and 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 (log density) as an additive sum of weak learners. Welling et al. [2002] in particular also approach boosting from the point of view of functional optimization of the likelihood or the logistic loss of an energy-based model. However, they rely on a first order local approximation of the objective since they focus on parametric weak learners such as restricted boltzman machines for which a second order step would be impractical.Greedy Multiplicative Boosting Another more direct multiplicative boosting framework was first proposed by Tu [2007] . At each boosting round a discriminative classifier is trained to distinguish between empirical data and data generated by the current model by estimating the likelihood ratio p(x) /qt(x). This estimated ratio is used as a direct multiplicative factor to update the current model q t (after being raised to an appropriate step size). In ideal conditions this greedy procedure would converge in a single iteration if a step size of 1 would be used. While Tu [2007] does not prescribe a particular choice of classifier to use, Grover and Ermon [2017] proposes a similar concept where the ratio is estimated based on an adversarial bound for an f -divergence and Cranko and Nock [2019] provides additional analysis on this method. In Appendix C we dive deeper into the differences between NRGBoost and this approach when it is adapted to use trees as weak learners. We note, however, that the main difference is that NRGBoost attempts to update the current density proportionally to an exponential of the ratio, exp (α t • p(x) /qt(x)), instead of the ratio directly.Tree-Based Density Modelling Other authors have proposed tree-based density models similar to DET [Nock and Guillame-Bert, 2022] or additive mixtures of tree-based models [Correia et al., 2020 , Wen and Hang, 2022 , Watson et al., 2023] but perhaps surprisingly, the natural idea of creating an ensemble of DET models through bagging has not been explored before as far as we are aware. Two distinguishing features of some of these alternative approaches are: (i) Unlike DETs, the partitioning of each tree is not driven directly by a density estimation goal. Correia et al. [2020] leverages a standard discriminative Random Forest, therefore giving special treatment to a particular input variable whose conditional estimation drives the choice of partitions and Wen and Hang [2022] proposes using a mid-point random tree partitioning. (ii) Besides modelling the density function as uniform at the leaf of each tree, other authors propose leveraging more complex models [Correia et al., 2020 , Watson et al., 2023] which can allow for the use of trees that are more representative with a smaller number of leaves. (iii) Nock and Guillame-Bert [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. Note that, unlike with boosting, in these approaches the new model doesn't add to the previous one but replaces it instead.  
experiments: 
 For our experiments we use 5 tabular datasets from the UCI Machine Learning Repository [Dheeru and 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) available through the Scikit-Learn package [Pedregosa et al., 2011] . 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 not possible with structured tabular data, and provides an example of the performance that can be achieved in an unstructured dataset with many features that are correlated among themselves. More details about these datasets are given in Appendix E.We split our experiments into two sections, the first to evaluate the quality of density models directly on a single variable inference task and the second to investigate the performance of our proposed models when used for sampling.In this section we test the ability of a generative model, trained to learn the density over all input variables, q(x), to infer the value of a single one. I.e., we wish to test how good is its estimate of q(x i |x -i ). For this purpose we pick x i = y as the original target of the dataset, noting that the models that we train do not treat this variable in any special way, except for the selection of the best model in validation. As such, we would expect that the model's performance in inference over this particular variable is indicative of its strength on any other single variable inference task and also indicative of the quality of the full q(x) from which the conditional probability estimate is derived.We use XGBoost [Chen and Guestrin, 2016] as a baseline for what should be achievable by a very strong discriminative model. Note that this model is trained to maximize the discriminative likelihood,E x∼p log q(x i |x -i ),directly, not wasting model capacity in learning other aspects of the full data distribution. As another generative baseline we use our own implementation of RFDE [Wen and Hang, 2022] which allows us to gauge the impact of the guided partitioning used in the DEF models over a random partitioning of the input space.We use random search to tune the hyperparameters of the XGBoost model and a grid search to tune the most important hyperparameters of the generative density models. We employ 5-fold cross-validation, repeating the hyperparameter tuning on each fold for all datasets except for the largest one (CT) for which we report results on a single fold. For the full details of the experimental protocol please refer to Appendix F. We find that NRGBoost performs better than the additive ensemble models (see Table 1 ) despite producing more compact ensembles. It often achieves comparable performance to XGBoost on the smaller datasets and with a small gap on the three larger ones. We note also that for the regression datasets the generative models provide an estimate of the full conditional distribution over the target variable rather than a point estimate like XGBoost. While there are other variants of discriminative boosting that also provide an estimate of the aleatoric uncertainty [Duan et al., 2020] , they rely on a parametric assumption about p(y|x) that needs to hold for any x.In this section, we compare the sampling performance of our proposed methods to neural-networkbased methods TVAE [Xu et al., 2019] and TabDDPM [Kotelnikov et al., 2022] on two metrics.The 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., 2022 , Borisov et al., 2022b] . It relies on using samples from the model to train a discriminative model which is then evaluated on the real data. Note that this is similar to the single variable inference performance from Section 6.1. In fact, if the density model's support covers that of the full data, one would expect the discriminative model to recover the generator's q(y|x), and therefore its performance, in the limit where infinite generated data is used to train it.We use an XGBoost model (with the hyperparameters tuned in real data) as the discriminative model and train it using a similar number of training and validation samples as in the original data. For the density models, we generate samples from the best model found in the previous section and for non-density models we select their hyperparameters by evaluating the ML Efficiency in the real validation set. Note that this leaves the sampling models at a potential advantage since the hyperparameter selection is based on the metric that is being evaluated rather than the direct inference performance of the previous section.Discriminator Measure Similar to Borisov et al. [2022b] we test the capacity of a discriminative model to distinguish between real and generated data. We use the original validation set as the real part of the training data in order to avoid benefiting generative methods that overfit their original training set. A new validation set is carved out of the original test set (20%) and used to tune the hyperparameters of an XGBoost model which we use as our choice of discriminator, evaluating its AUC on the remainder of the real test data.We repeat all experiments 5 times, with 5 different generated datatsets from each model. Results are reported in Tables 2 and 3 showing that (i) NRGBoost outperforms all other methods by substantial margins in the discriminator measure except for the PR and the MBNE datasets. (ii) On the ML Efficiency metric, TabDDPM outperforms NRGBoost by small margins on the small datasets which could in part be explained by the denser hyperparameter tuning favouring models that perform particularly well at inferring the target variable at the expense of the others. Nevertheless, NRGBoost still significantly outperforms all other models on MNIST and CT. Its samples also look visually similar to the real data in both the MNIST and California datasets (see Figures 1 and 2 ). Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: The experimental setup is described in as much detail as the space allows in Section 6. The full setup is described in Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
discussion: 
 While the additive tree models like DEF require no sampling to train and are easy to sample from, we find that in practice they require very deep trees to model the data well which, in turn, also requires using a large number of trees in the ensemble to regularize. In our experiments we found that their performance was often capped by the maximum number of leaves we allowed them to grow to (2 14 ).In contrast, we find that NRGBoost is able to model the data better while using shallower trees and in fewer number. Its main downside is that it can only be sampled from approximately using more expensive MCMC and also requires sampling during the training process. While our fast Gibbs sampling implementation coupled with our proposed sampling approach were able to mitigate the slow training, making these models much more usable in practice they are still cumbersome to use for sampling due to autocorrelation between samples from the same Markov Chain. We argue however that unlike in image or text generation where fast sampling is necessary for an interactive user experience, this can be less of a concern for the task of generating synthetic datasets where the one time cost of sampling is not as important as faithfully capturing the data generating distribution.We also find that tuning the hyperparameters of tree-based models is easier and less crucial than DL models for which many trials fail to produce a reasonable model. In particular we found NRGBoost to be rather robust, with different hyperparameters leading to small differences in performance.Finally, we note that like any other machine learning models, generative models are susceptible to overfitting and are thus liable to leak information about their training data when generating synthetic samples. In this respect, we believe that NRGBoost offers better tools to monitor and control overfitting than other alternatives (see Section 3.3) but, still, due consideration for this risk must be taken into account when sharing synthetic data.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: Sample standard deviations for all experiments are reported in Appendix G.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
conclusion: 
 In this work, 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 better overall sampling performance than alternatives. We hope that these results encourage further research into generative boosting approaches for tabular data, in particular exploring other applications besides sampling that are enabled by density models.Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: This information is provided in Appendix G.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: As far as we are aware there are no violations of the Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: We discuss the main potential misuse of our work in Section 7.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We do not believe our proposed models have a high risk of misuse but will nonetheless highlight the potential risks in the documentation when we release the code.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: As far as we are aware we cite all the sources of the data used in our experiments as well the main software packages used.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.• For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.• If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: We don't release any new assets at the time of submission. We plan to release the code later and will fully document it.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.• The paper should discuss whether and how consent was obtained from people whose asset is used.• At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: We don't conduct any experiments involving human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.• Justification: We don't conduct any experiments involving human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. 
 
 Appendix: The expected log-likelihood for an energy-based model (EBM),EQUATIONis given byEQUATIONThe first variation of L can be computed asEQUATION)This is a linear functional of its second argument, g, and can be regarded as a directional derivative of L at f along a variation g. The last equality comes from the following computation of the first variation of the log-partition function:EQUATIONAnalogous to a Hessian, we can differentiate Equation 15 again along a second independent variation h of f yielding a symmetric bilinear functional which we will write as δ 2 L[f ; g, h]. Note that the first term in equation 2 is linear in f and thus has no curvature, so we only have to consider the log partition function itself:EQUATION)Note that this functional is negative semi-definite for all f , i.e. δ 2 L[f ; h, h] ≤ 0, meaning that the log-likelihood is a concave functional of f .Using these results, we can now compute the Taylor expansion of the increment in log-likelihood L from a change f → f + δf up to second order in δf :EQUATIONAs an aside, defining the functional derivative, δJ[f ] δf (x) , of a functional J implicitly by:EQUATIONwe can formally define, by analogy with the parametric case, the Fisher Information "Matrix" (FIM) at f as the following bilinear functional of two independent variations g and h:EQUATIONThe only difference to the second-order variation of 2 computed in equation 20 would be that the expectation is taken under the model distribution, q f , instead of the data distribution p. However, because the only term in log q f (x) that is non-linear in f is the log-partition functional, which is not a function of x, this expectation plays no role in the computation and we get the result that the FIM is the same as the negative Hessian of the log-likelihood for these models.Considering a weak learner such asEQUATIONwhere the subsets X j are disjoint and cover the entire input space, X , we have thatEQUATIONSimilarly, making use of the fact that 1 Xi (x)1 Xj (x) = δ ij 1 Xi (x), we can computeEQUATIONIn fact, we can extend this to any ordinary function of δf :EQUATIONwhere we made use of the fact that the 1 Xj constitute a partition of unity:EQUATIONFinally, we can compute the increase in likelihood from a step f → f + α • δf asEQUATIONwhere in equation 42 we made use of the equality:log Z[f + α • δf ] -log Z[f ] = log x exp(f (x) + αδf (x)) Z[f ] = log x q f (x) exp(αδf (x)) ,(44) and of the result in equation 39 in the final step.This result can be used to conduct a line search over the step size using training data and to estimate an increase in likelihood at each round of boosting for the purpose of early stopping, using validation data.Density Estimation Trees (DET) [Ram and Gray, 2011] model the density function as a piecewise constant function,EQUATIONwhere X j are given by a partitioning of the input space X induced by a binary tree and the v j are the density values associated with each leaf that, for the time being, we will only require to be such that q(x) sums to one.Ram and Gray [2011] proposes fitting DET models to directly minimize a generative objective, the Integrated Squared Error (ISE) between the data generating distribution, p(x) and the model:EQUATION)Noting that q is a function as in Equation 45 and that J j=1 X j = X , we can rewrite this as min v1,...,v J ,X1,...,X J x∈XEQUATIONSince the first term in the objective does not depend on the model this optimization problem can be further simplified asEQUATIONwhere V (X) denotes the volume of a subset X. Solving this quadratic program for the v j we obtain the following optimal leaf values and objective:EQUATION)One can therefore grow a tree by greedily choosing to split a parent leaf with support X P into two leaves with supports X L and X R so as to maximize the following criterion:EQUATIONTo maximize the likelihood,EQUATIONrather than the ISE one can use the same approach. Here the optimization problem to solve is:EQUATIONThis is, again, easy to solve for v j since it is separable over j after removing the constraint using Lagrange multipliers. The optimal leaf values and objective are in this case:EQUATION)The only change is therefore to the splitting criterion which should become:EQUATIONIn multiplicative generative boosting an unnormalized current density model, qt-1 (x), is updated at each boosting round by multiplication with a new factor δq αt t (x):EQUATIONFor our proposed NRGBoost, this factor is chosen in order to maximize a local quadratic approximation of the log likelihood around q t-1 as a functional of the log density (see Section 3). The motivation behind the greedy approach of Tu [2007] or Grover and Ermon [2017] is to instead make the update factor δq t (x) proportional to the likelihood ratio r t (x) := p(x) /qt-1(x) directly, which under ideal conditions would mean that the method converges immediately when choosing a step size α t = 1. In more realistic setting, however, this method has been shown to converge under conditions on the performance of the individual δq t as discriminators between real and generated data [Tu, 2007 , Grover and Ermon, 2017 , Cranko and Nock, 2019] .While in principle this desired r t (x) could be derived from any binary classifier that is trained to predict a probability of a datapoint being generated (e.g., by training it to minimize a strictly proper loss) and Tu [2007] does not prescribe any particular choice, Grover and Ermon [2017] propose relying on the following variational bound of an f -divergence to derive an estimator for this ratio:EQUATION)Here f * denotes the convex conjugate of f . This bound is tight, with the optimum being achieved foru * t (x) = f ′ ( p(x) /qt-1(x)), if U t is capable of representing this function. (f ′ ) -1 (u * t (x)) can thus be interpreted as an approximation of r t (x).Adapting this method to use trees as weak learners can be accomplished by considering U t in Equation 56to be defined by tree functions u = 1 /J J j=1 w j 1 Xj with leaf values w j and leaf supports X j .At each boosting iteration a new tree, u * t can thus be grown to greedily optimize the lower bound in the r.h.s. of Equation 56 and setting δq t (x) = (f ′ ) -1 (u * t (x)) which is thus also a tree with the same leaf supports and leaf values given by v j := (f ′ ) -1 (w j ). This leads to the seaprable optimization problem: EQUATIONP /Q DiscBGM (χ 2 ) P 2 /Q P /Q NRGBoost P 2 /Q exp ( P /Q -1)Note that we drop the iteration indices from this point onward for brevity. Maximizing over w j with the X j fixed we have that w * j = f ′ ( P (Xj ) /Q(Xj)) which yields the optimal valueJ * (X 1 , . . . , X j ) = j P (X j )f ′ P (X j ) Q(X j ) -Q(X j )(f * • f ′ ) P (X j ) Q(X j ) (58)that in turn determines the splitting criterion as a function of the choice of f . Finally, the optimal density values for the leaves are given byEQUATION)It is interesting to note two particular choices of f -divergences. For the KL divergence, f (t) = t log t and f ′ (t) = 1 + log t = (f * ) -1 (t). This leads toEQUATIONas the splitting criterion. The Pearson χ 2 divergence, with f (t) = (t -1) 2 , leads to the same splitting criterion as NRGBoost. Note however that for NRGBoost the leaf values for the multiplicative update of the density are given by exp ( P (Xj ) /Q(Xj) -1) instead of the ratio directly. Table 4 summarizes these results.Another interesting observation is that a DET model can be interpreted as a single round of greedy multiplicative boosting starting from a uniform initial model. The choice of the ISE as the criterion to optimize the DET corresponds to the choice of Pearson's χ 2 divergence and likelihood to the choice of KL divergence.Discretization In our practical implementation of tree based methods we first discretize the input space by binning continuous numerical variables by quantiles. Furthermore we also bin discrete numerical variables in order to keep their cardinalities smaller than 256. This can also be interpreted as establishing a priori a set of discrete values to consider when splitting on each numerical variable and is done for computational efficiency, being inspired by LightGBM [Ke et al., 2017] .Categorical Splitting For splitting on a categorical variable we once again take inspiration fromLightGBM. Rather than relying on one-vs-all splits we found it better to first order the possible categorical values at a leaf according to a pre-defined sorting function and then choose the optimal many-vs-many split as if the variable was numerical. The function used to sort the values is the leaf value function. E.g., for splitting on a categorical variable x i we order each possible categorical value k by P (xi=k,X-i) / Q(xi=k,X-i) in the case of NRGBoost where X -i denotes the leaf support over the remaining variables.Tree Growth Strategy We always grow trees in best first order. I.e., we always split the current leaf node that yields the maximum gain in the chosen objective value.Line Search As mentioned in Section 3, we perform a line search to find the optimal step size after each round of boosting in order to maximize the likelihood gain in Equation 43. Because evaluating multiple possible step sizes, α t , is inexpensive, we simply do a grid search over 101 different step sizes in the range [10 -3 , 10] with their logarithm uniformly distributed. Random Forest Density Estimation (RFDE) We implement the RFDE method [Wen and Hang, 2022] after quantile discretization of the dataset and therefore split at the midpoint of the discretized dimension instead of the original one. When a leaf support has odd cardinality over the splitting dimension a random choice is made over the two possible splitting values. Finally, the original paper does not mention how to split over categorical domains. We therefore choose to randomly split the possible categorical values for a leaf evenly as we found that this yielded slightly better results than a random one vs all split.Code Our implementation of the proposed tree-based methods is mostly Python code using the NumPy library [Harris et al., 2020] . We implement the tree evaluation and Gibbs sampling in C, making use of the PCG library [O'Neill, 2014] for random number generation.We use 5 datasets from the UCI Machine Learning Repository [Dheeru and Karra Taniskidou, 2017] :Abalone, Physicochemical Properties of Protein Tertiary Structure (referred to as Protein in the sequence), 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 as well as the approximate number of samples used for train/validation/test for each cross-validation fold.To tune the hyperparameters of XGBoost we use 100 trials of random search with the search space defined in Table 6 . , 32, 64, 128, 256, 512 , 1024}) min_child_weight LogUniform 10 -1 , 10 3 reg_lambda 0.5 • δ(0) + 0.5 • LogUniform 10 -3 , 10 reg_alpha 0.5 • δ(0) + 0.5 • LogUniform 10 -3 , 10 max_leaves 0 (we already limit the number of leaves) grow_policy lossguide tree_method hist The fraction of features to consider when splitting a node as a function of the total number of features d[d -1/2 , d -1/4 , 1]min_data_in_leaf The minimum number of data points that need to be left in each leaf for a split to be considered , 3, 10, 30] product of training marginals (90%) on the discretized input space. We observed that this mixture coefficient does not have much impact on the results however.[0, 1We only tune two parameters for NRGBoost Models:• The maximum number of leaves for which we try the values [64, 256, 1024, 4096] in order, stopping if performance fails to improve from one value to the next. For the CT dataset we also include 16384 in the values to test.• The constant factor by which the optimal step size determined by the line search is shrunk at each round of boosting. This is essentially the "learning rate" parameter. To tune it we perform a Golden-section search for the log of its value using a total of 6 evaluations. The range we use is [0.01, 0.5].This means that at maximum we train only 24 NRGBoost models (30 for CT).All other relevant parameters are fixed and their values, along with a short description, is given in Table 9 . Single variable inference For the single variable inference evaluation, the best models are selected by their discriminative performance on a validation set. The entire setup is repeated five times with different cross-validation folds and with different seeds for all sources of randomness except on the CT dataset due to its large size. For the Adult and MNIST datasets the test set is fixed but training and validation splits are still rotated.Sampling For the sampling evaluation we use a single train/validation/test split of the real data (corresponding to the first fold in the previous setup) for training the generative models. The density models used are those previously selected based on their single variable inference performance on the validation set. For the sampling models (TVAE and TabDDPM) we directly evaluate their In Table 13 we report the best hyperparameters found for NRGBoost for the first cross-validation fold together with the time taken to train this best model.We will assume that X is finite and discrete to simplify the notation and exposition but everything is applicable to bounded continuous input spaces, replacing the sums with integrals as appropriate. 


Decision: The authors propose tree-based methods for generative modeling.
- **Reviews**: The reviews highlight several weaknesses, particularly criticizing the choice of datasets, fixed hyperparameters, and the absence of comparative methods (in particular, well-known generative alternatives for tabular data are missing).
- **Rebuttal**: The authors provided new results and included one of the suggested methods by reviewer `trn6` for comparison.
- **Reviewer-AC Discussion Phase**: During this phase, reviewers `UJVb` and `trn6` expressed concerns. While their scores are borderline (reject and accept, respectively), they indicate that the authors' rebuttal did not address all open issues and that their most significant concerns remain unresolved. Furthermore, both explicitly state that they believe the paper, in its current form, should not be accepted.
- **Decision**: While the paper is well-written, addresses an interesting problem, and has an above-average score, the latter is primarily due to the score by Reviewer `fQFZ` and, as `UJVb` noted, might not reflect the consensus among reviewers. Given that two out of four reviewers recommend against accepting the paper in its current form, I recommend rejecting the paper.

Reviews: review: paper_summary: paper_summary: The paper proposes to extend the success of tree-based methods in discriminative tasks to generative modelling, which is implemented via an energy-based generative boosting algorithm (NRGBoost). Specifically, NRGBoost directly extends the tree-based tabular models by replacing the discriminative objectives with a generative one, which seems novel to me.
strength_weakness: strength_weakness: 1. **[Important]** Some highly relevant benchmark methods are missing, including ARF (tree-based) [1], GOGGLE (diffusion) [2] and TabPFGen (energy-based) [3].

2. In Line 325, the authors claim that the proposed method “significantly” outperforms other methods, while the significance test seems missing.

3. I would suggest the authors add comparison results on the computation efficiency. Because NRGBoost basically employs the same architecture as traditional gradient boosting trees, the computation efficiency should be higher than most other network-based generative models.

4. There seem to be some typos throughout the main text: “I.e.” (Line 272)

5. **[Important]** Code is not provided. I remain conservative about the results claimed in the paper.

[1] Watson, David S., et al. "Adversarial random forests for density estimation and generative modeling." International Conference on Artificial Intelligence and Statistics. PMLR, 2023.

[2] Liu, Tennison, et al. "GOGGLE: Generative modelling for tabular data by learning relational structure." The Eleventh International Conference on Learning Representations. 2023.

[3] Ma, Junwei, et al. "TabPFGen–Tabular Data Generation with TabPFN." NeurIPS 2023 Second Table Representation Learning Workshop.
questions: questions: 1. How does NRGBoost perform single variable inference (Lines 271-277)?
limitations: limitations: The authors detail the limitations of NRGBoost in Section 7.
paper_summary: paper_summary: This paper proposed an energy based generative boosting algorithm analogous to XGBoost, which can be used as generative model as well as be applied to discriminative tasks.
strength_weakness: strength_weakness: My concerns regarding the proposed method and the experiments are detailed in the questions section.
questions: questions: 1. For Table 1:
Why is NRGBoost shown in bold when it's overperformed by XGBoost? In other words why is XGBoost used as a baseline?
Is there a particular reason for SE reported separately in the supplement? 

2. For the sampling experiments in section 6.2, is the authenticity of synthetic samples measured in some ways? (In an extreme case, if the synthetic samples are almost replicating training data, they will achieve high ML efficiency and good discriminator measures)

4. How is the scalability? 

5. Can the proposed method be used on (unordered) categorical data, which is very common in tabular data?
limitations: limitations: I do not identify significant limitations other than the ones discussed in the paper.
paper_summary: paper_summary: This paper proposes a boosted tree algorithm that performs distribution learning using an energy-based formulation. Inspired by methods like XGBoost, it is claimed to achieve high performance not only in generative ability but also in discriminative performance.
strength_weakness: strength_weakness: This study is not theoretical; therefore, the validity of the proposed method must be confirmed through experiments. However, there are some questions regarding the experimental settings. 

Also, from an algorithmic perspective, since methods unique to tree ensembles are incorporated, they could potentially offer advantages in terms of computational complexity compared to other methods. However, evaluations regarding efficiency have not been conducted. If there are advantages, it seems opportunity loss. 

Please refer to the Questions section.

Minor:
The evaluation of variance is presented separately in Appendix G, but it is difficult to compare. Therefore, I would like it to be summarized in one table.
questions: questions: **(1) Robustness**

In the discussion, it is stated that hyperparameters are not critical compared to deep learning, but it seems that there is no experimental evidence to support this. If you claim robustness to hyperparameters as a strength, please provide evidence for this.

**(2) Computational cost**

Please describe the computational cost and compare it with benchmarks. In the discussion, DEF is mentioned to not require sampling for training and that the processing for sampling is light, but please quantify how much of a difference there is.

**(3) Dataset**

How do you select datasets? It seems that only a few are picked from among the many datasets available, which may give the impression that you are cherry-picking the data that works well.

**(4) Fixed parameters**

It seems that the hyperparameters are fixed in advance as shown in Table 9, but how were these decided? It is mentioned that the hyperparameter search can be covered by only 24 combinations, but this also gives the impression that you are setting them in advance to work well and narrowing the apparent search space.

**(5) Parameter range**

It was also mentioned that DEF tends to have a large number of leaves, reaching 2^14, but isn't this because the search range is set that way? Even if you align this setting with NRGBoost, which has a maximum of 2^12, will it lead to the same conclusion?

**(6) Benchmark performance**

In Figure 1, the performance of the benchmark methods seems too poor. Is this figure created with the same settings as those measured in Table 2, for example? TabDDPM seems to show better performance than the proposed method in the ML Efficiency result, but it is hard to imagine getting similar results from the generated results in Figure 1.
limitations: limitations: Computational cost of sampling is larger than existing methods.
paper_summary: paper_summary: The authors propose a energy-based generative boosting method. They try to maximize the log-likelihood functional delta with a second-order expansion. This lead to a boosting algorithm with deltas as steps in the log-likelihood. Instead of scaling the steps with a fixed predefined value (as a anti-overfitting hyperparameter), they rely on linear search to get better performance. They initialize f0 as uniform and use trees as weak learner to learn the deltas. They derive the objective with trees, it is well explained and looks similar to other second-order method objectives like XGBoost. They use MCMC to sample from Q(x) as is typical from energy-based approaches. They use Gibbs so they only need to sample from one dimension while keeping the rest constant. However since there are t trees, this is quadratic in t. They use some form of accept/reject to sometimes accept previous iteration samples in order to not have to re-samples new samples given the quadratic cost. They include a probability of refresh in order to not just always accept old samples. They propose interesting ways of regularizing the approach. They provide good literature review of related methods. It seems like Section 4 could be integrated with the related work section.

Figure may appear unimpressive to a generative expert unfamiliar with trees. But, being able to generate good samples from MNIST using only decision is an impressive feat (even if MNIST is downsampled). This paves the way for trees being used on more complex data.

Table 2 shows nice prediction results, and they do extensive hyperparameter tuning. However, ML Efficiency is not the best metric, it only focus on prediction. A generative model could produce low-diversity fake data that lead to good classifiers. I would recommending adding some distribution metrics such as the Wasserstein distance (which is computable in low-dimension and is used for high-dim data such as images/videos in the form of the widely popular FID/FVD). There are also other useful metrics for quality and diversity, I recommend looking the extensive choice of tabular-data metrics described in https://proceedings.mlr.press/v238/jolicoeur-martineau24a.html. However, the only essential one in my opinion is having a distance in distribution. It can be the Wasserstein distance or something like the MMD distance.

Except for the missing distribution metric, everything else in the results is good and the methodology is sound with good hyperparameter tuning. I would encourage the authors to add a distribution metrics.

The method is sound, novel and quite interesting. The presentation is very well made. Being the first tree method achieving good image data samples is impressive in my opinion. This is a very high quality paper.
strength_weakness: strength_weakness: Missing a distribution metric, the ML efficiency is not a adequate metric on its own.
questions: questions: What is the memory requirements, especially when processing large datasets such as MNIST? Could you provide sample time and RAM requirements for the datasets? How hard is scaling with respect to N_samples and also scaling with respect to N_features?
limitations: limitations: If the authors had to downsample MNIST from 28x28 to lower, than there are scaling limitations that should be added.

Comments: title: Response to rebuttal
GOGGLE comparison
Thank you
Additional Experiment on Sample Memorization
Response to Reviewer Comment (2/2)
Response to Reviewer Comment (1/2)
Thank you!
response
Additional comments on MNIST results
Fixed Hyperparameters
comment: We thank the reviewer for their time and regret that our addition of the ARF comparison was not sufficient to address their concerns about a complete benchmark. 

We decided to focus our efforts on ARF because it is also a density estimator and thus, the more interesting baseline given that it could be used in the same types of applications as NRGBoost. Unfortunately, this required implementing single variable inference for this method ourselves which consumed a significant portion of our rebuttal week (together with running hyperparameter tuning on 5 cross-validation folds for every dataset). 

For the reasons we explained above, we don't believe that TabPFGen would be an appropriate comparison for our setup, even if we were able to run it on datasets of the size that we use.

Furthermore, assuming that we could resolve the performance issues we are currently experiencing with GOGGLE, due to the method being so slow to train we would have to severely limit the number of epochs relative to what was used in the original paper in order to bring training times in line with other methods. We will try to add a comparison along these lines in the future but we are not sure how valuable it will be given these constraints.

Best regards,

The authors
Thank you for your detailed response and the efforts you've put into addressing the concerns raised. After carefully considering your points and other reviewers' comments, I still believe that my initial evaluation remains valid: a fair and complete benchmark should be necessary before acceptance, so my score remains unchanged.
Thank you for aknowledging our response. We want to reiterate that we are open to discussing any remaining concerns you might have.
Thanks for your response. I have read your response.
We apologize for the extra comment but given the fast approaching end of the discussion period we thought it would be better to antecipate any questions the reviewer might have and provide an update on our ongoing efforts to add GOGGLE as a baseline as the reviewer originally requested.

We have managed to get a working version of the synthcity library with the GOGGLE plugin (and we can't stress enough that this was not an easy task). 
So far, we have only managed to train models with (mostly) the library's default hyperparameters on 5 of the datasets but, in the process, have faced significant challenges which we outline below.

# Slow training
The method is considerably slower than all the other methods we compare:
- As an extreme example, training a model for 500 epochs on `covertype` took **more than 12 hours** at the maximum batch size used in the paper (128).
Note that the paper (and the library default) uses 1000 epochs (with early stopping) which could lead to more than 1 full day to train a single model!
- We haven't been able to run GOGGLE on the `MiniBooNE` or `MNIST` datasets as it seems that it **scales poorly with the number of features in the dataset**. 
For context, **a single epoch on the `MiniBooNE` dataset takes 17 minutes to complete**.
- Even for the datasets where we are able to train models, the slow training times render us unable to do hyperparameter tuning in a reasonable timeframe like we do for all the other methods.

# Lackluster Performance
While we would like to explore the hyperparameter space for this model better, the results we have obtained so far with the default hyperparameters have yielded poor performance in general.

On all of our regression datasets we consistently get $R^2$ values below 0 for the ML efficiency metric.
While we are suspicious of these results and will need to investigate further we note that the paper appears to focus on binary classification tasks so we have no frame of reference.

On classification datasets, trained models tend to only generate data from majority classes:
- On the adult dataset, using our training setup, the trained model exclusively produces samples of the majority class.
- On the covertype dataset, it also rarely produces samples of any of the classes outside of the top two (out of 7 classes). This makes performance on this dataset poor as well.
- Even if we could conditionally generate based on the class, this goes against the idea of our setup which, as we explained above, is meant to produce models that are good for any inference task over any potential variable, not just a specific one. It would, furthermore, make it unfair to all other methods for which we don't use conditional generation (even though we could for some of them, including NRGBoost).

# Conclusion
We will continue to investigate our options and explore different choices of hyperparameters to try to achieve a setup that we believe is fair, both to GOGGLE and the other methods.
As it stands, given the long training times, we do not believe we will be able to provide a meaningful comparison with this baseline; at least not in any of the larger datasets.
However we would very much welcome any suggestions or recommendations from the reviewer on this subject.

We would also be very thankful if the reviewer could provide any feedback on whether their concerns have been adequately addressed and offer again to clarify any remaining questions that they might have.
We have made a significant effort to comply with every request of the reviewer to the best of what we believe is reasonably possible and we think that this has significantly improved our paper as a result. In light of this we would kindly ask the reviewer to consider re-evaluating their score.

Best regards,

The authors
We are glad we were able to address your concerns. Thank you for raising your score.
Thanks for the rebuttal and additional experiments. This addressed my concerns and I increase my score.
We apologize for the additional comment but the reviewer's question about memorizing the training data gave us an idea about an additional experiment that we believe supports our argumentation above and might interest the reviewer.

We introduce a **memorizing** baseline that simply memorizes a fraction of the training set.
When prompted for synthetic samples, it samples from this memorized set with replacement.
We try two settings for this baseline: 
- Memorizing 10% of the training set
- Memorizing the full training set

Finally we evaluate the same discriminator setup we use in the paper. 
This experiment allows us to gauge how well learning the training empirical distribution (with different sample sizes) performs.

Below we report the average discriminator measure results we obtained when repeating each experiment 5 times as in the paper. 
We also include the results for the two best methods (TabDDPM and NRGBoost) for ease of comparison.
Note that these results are ROC AUC for distinguishing synthetic from real (held out) samples and as such **lower is better**.


| Dataset         | AB | CH | PR |  AD  | MBNE | MNIST    | CT |
|-----------------|---------|------------|---------|---------|-----------|---------|-----------|
| TabDDPM         |   0.818 |      0.667 |**0.628**|  0.604  |     0.789 |  1.000  |     0.915 |
| NRGBoost        |**0.625**|   **0.574**|   0.631 |**0.559**|     0.993 |  0.943  |     0.724 |
| Memorize (10%)  |   0.996 |      1.000 |   0.995 |  1.000  |     0.959 |  0.960  |     0.958 |
| Memorize (100%) |   0.763 |      0.799 |   0.762 |  0.769  |  **0.607**|**0.611**|  **0.606**|

For the smaller datasets, even the best possible memorization baseline fails to beat the best method on each respective dataset.
For the larger datasets (MiniBooNE, MNIST and Covertype), while memorizing the full training set does beat all other methods, memorizing a more sensible fraction of it is not enough.
As we increase the memorized set, the empirical distribution becomes closer to the data generating distribution and therefore discriminator measure results improve but we hope that this shows that simply memorizing samples is often not enough to achieve best performance on this metric due to the failure to generalize.

We would like to kindly ask the reviewer if their concerns have been resolved and if not, what remaining concerns they have.
We would greatly appreciate this feedback in order to further improve our paper.

Best regards,

The authors
We agree with the reviewer that a fair benchmarking of all methods is important for any claim of superiority. 
In regards to tree-based models, we believe that we gave every possible advantage to our original baselines:
- We proposed DEF models in section 4 of our paper as an improvement over the original Density Estimation Trees (DET), introducing two algorithmic improvements (bagging and changing from ISE metric to KL divergence) in order to get them to perform well and to have a worthy baseline.
- Overall, DEF models took **longer** to run hyperparameter tuning and single models are usually slower to train. Despite this they still achieve lower performance consistently. 
- NRGBoost tends to scale much better with depth of the trees than DEF models and in many datasets we could reduce the number of leaves with only a small performance penalty. 
In contrast, we had to push DEF models to the limit of what our implementation would allow to even have reasonable performance. If we trained DEF models with trees of only 1024 leaves (typical for NRGBoost) performance would be abysmal in general.

Generally, the performance gap between NRGBoost and DEF models is so large that even if we were to further handicap NRGBoost by, e.g., reducing the number of trees and the maximum number of leaves it would still come out on top.
As a result we believe that our original statement that NRGBoost performs better than the additive tree-based ensemble models is accurate.


For neural network methods, because we are not experts, we chose to rely on guidance from the existing literature when it comes to hyperparameter tuning. 
But since these are such different methods, running on different types of hardware, and with very different hyperparameters and hyperparameter tuning requirements it can be particularly challenging to achieve an experimental setup that everyone will agree is fair. We try to find a setup that we believe is representative of what a practitioner may get when using these methods and, like for our tree-based baselines, generally try to err on the side of caution. As an example, we give neural networks methods the benefit of using ML Efficiency directly as the hyperparameter tuning metric which is one of the metrics they are later evaluated on.

Nonetheless we are commited to being as fair as possible when presenting our results and we value the feedback from the reviewer on this issue. We welcome any further suggestions or recommendations on this matter that could improve our paper.

Best regards,

The authors
We thank the reviewer for their time in reviewing our rebuttal and for being willing to discuss our paper.
Below we outline the computational complexities of all methods but note that these provide only asymptotic guidance.

## Tree-based methods
All tree based methods have to fit the ensemble of trees. The time complexity of this should be similar across all non-random methods (ARF, DEF, NRGBoost). Essentially, it scales as $O(NFTD)$ where $N$ is the number of data points, $F$ the number of features, $T$ the number of trees in the ensemble and $D$ the depth of the trees with the following caveats:
- For NRGBoost the factor $N$ should be replaced by $N + M$ where $M$ is the total number of samples in a pool of samples (see below). As mentioned previously in our comment on fixed hyperparameters, choosing $M$ similar $N$ works well in practice so the overall scaling **with the data characteristics** is still $O(NF)$.
- For ARF there is an additional multiplicative factor which is the number of times the adversarial loop runs. This is a complex function of the dataset as the algorithm stops when no improvement can be made. In practice we observe that it is somewhere between 2 and 6 iterations with smaller datasets generally stopping earlier.
### Sampling in NRGBoost
To fit a tree at each round of boosting, NRGBoost requires samples from the model. We keep a pool of samples of size $M$ and before each round of boosting we need to resample a portion of those which we will call $m$ and assume that it is roughly constant across rounds. Note that this sampling has nothing to do with implicit parallelization like the reviewers' statement seems to imply. It is merely an additional required step of our algorithm which we happen to be able to parallelize effectively. 

The total cost of sampling is on average $O(mFT^2 D)$. As can be seen this process also scales linearly with $F$ and $N$ (through $m$). The quadratic scaling with the $T$ hyperparameter is the reason we only train NRGBoost models with 200 trees vs the DEF models that are trained with 1000 trees.

**To summarize, the computational complexity of all tree-based methods scales linearly with both number of data points and number of features as far as the characteristics of the data are concerned.**

As for the variation across different datasets, this is mostly due to these characteristics, both directly and indirectly:
- **Directly:** Datasets with a larger $NF$ factor will be slower which is the main reason that MNIST is the slowest followed by either miniboone or covertype.
- **Indirectly:** The characteristics of the dataset also change what are the best hyperparameters. For tree-based models, the number of leaves (depth $D$) is perhaps the biggest factor where datasets that favor larger trees will have slower training times for the best model, but other hyperparameters also have an effect.

## Neural Networks
Computational complexity for neural network methods depends on the architecture of the network and type of method. While we are not experts on the matter we can make a few educated guesses.

As an example, just the matrix multiplication in an input layer for $T$ training steps with a batch size $B$ would be $O(TBFH)$ where $H$ is the number of hidden output units of this layer. 
This is similar to tree-based methods' scaling with $NF$ when we assume that a fixed number of epochs are used to train the model (i.e., $T$ is proportional to $N/B$). 
Depending on the method there would be other multiplicative factors such as the number of diffusion timesteps and there are also the other layers of course.
These other layers may be responsible for most of the training time in practice and their input and output sizes are hyperparameters that may be chosen not to scale with $F$.
These models can also leverage massive parallelism from the use of GPUs which makes these comparisons based on computational complexity not very helpful in our opinion.


# Summary

We believe all methods scale similarly with the characteristics of the data under reasonable assumptions. 
However this does not tell the full picture as they also scale differently with the very different hyperparameters of each method.
Controlling for computational time is difficult because there is hyperparameter tuning in the loop which might favor hyperparameters that lead to slower or faster models depending on the dataset. 

We will follow up with an additional comment justifying our benchmarking approach and our statements due to space limitations.
Thank you for your response. I would like to discuss complexity.

I believe that differences due to variations in implementation (e.g., programming languages) are not important. If all benchmarks are not implemented consistently during evaluation, the comparison becomes unequal, and no claim of superiority can be made. Additionally, it seems that implicit parallelization is also being done through techniques like sampling, so it is difficult to objectively understand the actual complexity by only looking at the computation time.

Since your paper is empirical, I believe a fair assessment can only be made by properly evaluating the trade-off between performance and cost.

Could you please describe the complexity of the proposed method and the benchmarks using Big-O notation and compare them? This might better demonstrate the validity of the proposed method rather than just computational speed.
According to the Rebuttal PDF, there are significant variations across different datasets. What is the cause of this?
Thank you once again for your time in reviewing our rebuttal and for increasing your score.

We are glad we were able to address your concerns.
I agree that the Discriminator score can be linked to a IPM metric loosely, but having a true metric like Wasserstein is better. I agree that there is no perfect way to normalize/scale the data when you have continuous and categorical data, this is a well known problem. The Gower distance specifically tries to tackle this problem and this the distance that you used for the Wasserstein distance (by min-max continuous and /2 one-hot), which at least makes the choice rational.

Thank you for addressing my requests for a distribution metric and memory concerns. Given this, I am increasing my score by 1.
We thought we should provide some additional comments on the performance of benchmark methods on MNIST which we could not fit in the rebuttal due to space limitations.

NRGBoost and the other methods we compare to are meant for **tabular data**. As such, they do not benefit from the helpful inductive biases of convolutional neural networks and have no awareness of the spatial relation between pixels.
They are in fact invariant to a permutation of the pixels and need to learn the complex relationships between every pixel and all other pixels. 
The dimensionality of this dataset is also higher than other typical benchmarks for tabular data.
All these properties make this problem very challenging for tabular models and (as reviewer **fQFZ** also seems to agree) we consider it a big achievement for NRGBoost to do as well as it does. 
We do not believe that other methods underperforming on this dataset is unexpected as it's not really the type of data that tabular models are designed for. Still, we thought it would be interesting to include it as it allows us to visually assess actual samples from the model.

Still, regarding the performance of competing methods, we note that downsampled MNIST was also used as a benchmark in [[Xu et al., 2019]](https://arxiv.org/abs/1907.00503) which introduces TVAE and CTGAN. However, samples from neither of the methods are ever shown in the paper.
On the request of reviewer **trn6**, we have also included a comparison to a tree-based generative model (ARF) which seems to perform better than our original baselines (please refer to the PDF attached to the author's rebuttal).

Finally, ML Efficiency can be very deceptive about the ability of the model to accurately capture the input distribution, $p(x)$. 
The most important property to achieve a high score in ML efficiency is that the model captures well the relationship between the target $y$ variable and the remaining input variables $x$ (i.e., has a good internal model for $p(y \vert x)$ since this is what the downstream discriminative model is trying to estimate). 
Modeling well the remaining $p(x)$ (which is what is responsible for the images themselves) is secondary because a wrong $p(x)$ only means that the samples are inefficiently placed for optimal learning of the discriminative function. 
Reviewer **fQFZ** makes a similar point in their review that models can produce low diversity data that still lead to good classifiers.

We note also that it is surprisingly easy to score well in **classification** on MNIST. Input features (i.e. pixels) are very redundant for predicting $y$. Even just modeling the relationship between a small subset of pixels and the target variable y can be enough to achieve a much better than random accuracy. As an example, using only the **8 pixels** in a central vertical line XGBoost is already able to achieve an accuracy of 0.649. Therefore a model presumably doesn't even need to model well the relationship between $y$ and all pixels, just a small subset in order to achieve a good ML efficiency score.

We hope this clarifies some of the reviewer's doubts about the results on this dataset.
Below we explain our rationale for the choice of the fixed hyperparameters for NRGBoost:
- `num_rounds`: In earlier versions of the method, before we introduced importance sampling and rewrote our Gibbs sampling implementation in C, going above 200 rounds was prohibitevely expensive. While increasing this will tend to always improve results as long as the shrinkage parameter is also being tuned (and using early stopping), we decided to keep 200 rounds as the overall cost of sampling during training still scales quadratically with this number. Furthermore, this value is in line with values that already work well for discriminative boosting algorithms like XGBoost or LightGBM (which has a default of 100).
- `splitter`: Inspired by LightGBM's default approach to always split the leaf that leads to the best improvement in the objective. This allows more expressive models with fewer leaves than alternative approaches like depthwise splitting (XGBoost default) which wastes splits on leafs that don't really lead to much improvement.
- `line_search`: As explained in the paper this typically allows the model to converge faster by taking larger steps in the beginning of training where the quadratic approximation we optimize is farther from the likelihood. Incidently we find that the step size tends to settle closer to 1 as training progresses which is the optimal step size predicted by the quadratic approximation we try to optimize.
- `max_ratio_leaf`: We found early on with toy data that limiting this to small values tends to work better for regularization. Since this also limits the growth of the energy function at each round and therefore reduces the number of samples rejected from the pool we found no reason to increase it. In the end, this plays a similar role to the shrinkage parameter which we already tune so it was a bit redundant to tune both.

As for the sampling parameters, we have always used a similar number of samples to the original training set and this has always worked as a reasonable rule of thumb. Going over that rarely improved the results by much. But since it also never hurts to have more samples and the smaller datasets trained fast enough already, we kept the number 80000 fixed as a rough upper bound to the training set size for all datasets except for covertype (for which we used 4 times more to roughly match its #samples in training).

We chose to use 16 MCMC chains since these run in parallel and that is the number of virtual cores in the CPU used to run the experiments. For covertype we also multiplied the number of chains by 4 to keep the #samples/chain similar.

For the `burn_in` and `p_refresh` parameters we just picked seemingly reasonable values as they don't seem to have a large impact (other than `p_refresh` increasing the sampling time which is why we chose to set it at 10%).

We hope that this clarifies the choices made. Please let us know if you have any further questions.