Statistical Inference for Gradient Boosting Regression

We thank the reviewer for recognizing the strengths of the paper, and for offering very helpful suggestions. We respond below to your concerns. In particular, we offer a description on the advantages of each algorithm, and provide additional results comparing the confidence intervals with vanilla Boulevard. If our response sufficiently addresses your concerns, we would greatly appreciate it if you’d consider increasing your score.

Weaknesses

1. Methodological contribution past the CLT. We agree with the reviewer that constructing the baseline procedures for statistical inference are straightforward once the CLT is established. The intervals do indeed fall out of the CLT. The hypothesis test requires a little more insight, but also does so.

• Problem with the straightforward approach. Unfortunately, the straightforward procedures are not computationally tractable on large datasets commonly encountered when using gradient boosting! Although the CLT reduces the problem of inference to kernel ridge regression (with the right kernel), computing the intervals is an $O(n^3)$ problem. This is the chief difficulty that Zhou and Hooker (2022) face when constructing replication intervals for vanilla Boulevard.
  • They mention in Section 6.2 that: “Although the matrix inverse involved in the expression is too time consuming to compute for large $n$, it is possible to estimate it using Monte Carlo when $n$ and $d$ are small.”
• We do not have this problem. An important methodological contribution that we should have highlighted more is in making the intervals and tests practical and tractable via matrix sketching. We employ the Nystrom method and subsample $s$ landmarks to form a low-rank approximation of the kernel. This yields a solution that requires only $O(ns^2)$ time to precompute the kernel, and only $O(s^2)$ time for inference -- yielding practical statistical inference in near-linear time in the number of datapoints $n$.
  • With the Nystrom method of Musco and Musco (2017), we can choose $s$ to near-linear in the effective dimension $d_{\text{eff}}^\mu$ of a ridge regression problem with regularization parameter $\mu$ on the kernel matrix: $s = \widetilde{O}(d_{\text{eff}}^\mu)$.

2. Comparison to vanilla Boulevard. Certainly! We provide the requested comparison below, and will include it in the final version.

• Should I use Boulevard or Algorithm 1? Boulevard regularization is a special case of Algorithm 1 when no dropout is used $p=0$. Should we use dropout? There are good reasons to believe so.
  • Improvements in prediction error. In theory, as we discuss in Lines 59 and Corollary 4, Algorithm 1 achieves up to a fourfold improvement in terms of the asymptotic relative efficiency compared to vanilla Boulevard, when all other hyperparameters but the dropout rate are held equal.
    • In practice, this is similar to asking if one should use dropout with a neural network. There, it is almost always used, but when it is not, there are alternatives, such as weight decay, that also allow for regularization. With this analogy in mind, choosing $p$ can be a go-to knob for tuning the amount of regularization – though alternatives like tree depth and learning rate exist.
  • Improvements in statistical inference. From a coverage/statistical inference perspective, as we discuss in Line 339, there are drawbacks to using vanilla Boulevard instead of Algorithm 1. Setting $p \to 0$ empirically leads to wider intervals and less power in the hypothesis test for variable importance.
    • Why does this happen? We scale Algorithm 1 by $\frac{1+\lambda (1-p)}{\lambda}$ to correct the bias. As a result, the irreducible error will be scaled up as well -- but less so when $p$ is large. A higher dropout rate also introduces decorrelation between trees, suggesting a stabler variance estimate.
• Should I use Boulevard or Algorithm 2? Asking why one should use Algorithm 2 over Boulevard amounts to asking whether one wants to fit trees in parallel or not, and there are computational reasons (e.g. parallelism) to do so.
  • Improvements in prediction error. In theory, Algorithm 2 can achieve an unbounded improvement in relative efficiency over vanilla Boulevard as we take $\lambda \to 0$, as we discuss in Lines 295-296.
    • But this improvement is in terms of the number of boosting rounds, not ensemble size! This clarifies the intuition – boosting random forests (a strong learner) should lead to better prediction error than boosting a tree (a weak learner), with the same number of boosting rounds.
    • Our results in Figure 2 are somewhat unfair to Algorithm 2, as the x-axis is in the ensemble size and not the number of boosting rounds. The number of boosting rounds that Algorithm 2 encounters can be as few as 31 on Abalone and 50 on Wine Quality, in contrast to the 500 boosting rounds all other algorithms enjoy. Given the same number of boosting rounds, we would certainly expect an improvement.
  • Parallel computation. With modern multi-core processors and distributed computing, there is immense potential for exploiting the possibility of increased parallelism in boosting.

3. Comparison of confidence intervals. We provide the requested comparison below.
   • For a quantitative comparison, we provide a new figure that plots the width and coverage of our intervals against the dropout probability. We have lodged this with the AC. Unlike in Figure 3, we do not perform the conformal adjustment outlined in Lines 320 to 327, and we plot the coverage and width across all datapoints and repetitions for simplicity. Although not performing this leads to overcoverage in the PIs, we have to do so to compare the widths – otherwise any over-or-under-coverage would simply be adjusted to the nominal 0.95 level.
     • The mean CI and PI widths decrease as $p$ increases. The PI coverage remains close to 1 throughout (unlike in Figure 3, overcoverage occurs as we do not perform the conformal adjustment here), while the CI coverage decreases as $p$ increases -- as we see in the table below.
   • For a qualitative comparison, we also provide a version of Figure 1 that we have lodged with the AC. This uses the same hyperparameters, except for setting $p=0$. Empirically, it seems that the width of the intervals does indeed decrease in $p$, though this can come at a cost with respect to CI coverage (though achieving good CI coverage is difficult in nonparametric regression in general).

• We also examine test size and power as a function of dropout. Performance is strong for $p \in [0.2, 0.5]$, peaking at $p=0.5$. However, for $p \leq 0.1$, Type II error rises significantly (vanilla Boulevard has Type II error $0.43$), and for $p \geq 0.6$, Type I error becomes uncontrollable. We provide a detailed figure lodged with the AC.

Questions

1. The asymptotic variance is estimated. It is unclear whether the reviewer is referring to $\mathbf{r}_n$ or $\widehat{\sigma}^2$. We address both. The convergence in $\mathbf{r}_n$ in Theorem 2 is in the theorem statement, in the sense that this is the exact quantity used in the CLT. This is a happy result – showing that the estimation is baked into the result, and the empirical quantity is the quantity of interest. Regarding $\widehat{\sigma}^2$, we estimate this as the sample average of the residuals on a hold-out dataset, for which the convergence and consistency is clear.
2. The CIs are large. This is a known issue with nonparametric regression.

• In theory. Sadly, large intervals are expected due to the curse of dimensionality, as minimax optimal rates are exponentially dependent on dimension. Although Corollaries 1 and 2 show that our algorithm is suboptimal by a quadratic factor for Lipschitz functions, this is small compared to the exponential dependence in $d$ that no algorithm can escape without making further assumptions, e.g. on additive structure.
• In practice. Empirically, Figure 3 demonstrates our prediction intervals achieve comparable median widths to conformal baselines, while providing stronger conditional (rather than merely marginal) coverage guarantees. This should be viewed positively.

3. Size of Figure 3. We agree, thank you for pointing this out. We will enlarge it for the camera-ready version.

• Regarding MSE performance, there is no consistent winner. Our methods remain competitive: rounding to three decimals, GBT and LightGBM each win four times; SGBT wins three; Algorithm 1, random forests, and XGBoost each win twice; Algorithm 2 and Boulevard each win once.

We hope that this response has assuaged your concerns to your satisfaction. Again, if it has, we would greatly appreciate it if you’d consider increasing your score!

References

Musco and Musco (2016), Recursive Sampling for the Nystrom Method

Dear Reviewer AVWH,

Certainly. We appreciate you raising this concern in particular, and very much agree that highlighting this explanation within the paper would improve it. Thank you for your support!

We appreciate the referee's attention to the paper and their confirmation that our mathematical development is sound.

However, we do strongly feel that the reviewer's score is unreasonable, as is their framing of comparisons solely in terms of predictive performance. We think that inference should be regarded as an additional component -- you don't need it to predict well on aggregate, but it can be important for some applications.

On the importance of statistical inference.

To the question of "Why does one need statistical inference?" which we will interpret as meaning for these models specifically. Here we could provide a very considerable range of situations where uncertainty quantification would be relevant. By this we mean "if you collected a new data set and ran the model again, how different would your prediction be"?

• You might feel quite different about having knee reconstruction surgery if you were told that it had an 80% success rate empirically, than if you were told that the confidence interval for what this success rate might actually be ran from 40% to 95%.
• To be more concrete, if you had run XGBoost on geologic samples to give you a predicted yield for proposed rare earth mines, you might prioritise a mine with a smaller but more certain predicted yield than one with a larger prediction but an enormous confidence intervals.
• Conversely, a marketing campaign may want to avoid sending materials to someone who is very unlikely to respond in order not to antagonise them, but might want to send materials out to someone of whom we have very high uncertainty.

These are all immediately practical situations in which uncertainty about a predicted value is already relevant. This is without considering use in downstream tasks:

• Prioritizing data to collect to reduce uncertainty over all.
• Propagating uncertainty through to forecasts about the consequences of taking actions in a sequential decision-making poblem.
• Propagating uncertainty into explanation or interpretation methods which might also be relevant for.
• Model reduction, or testing whether expensive features are worth collecting given how they affect performance.

We note that in all these cases, inference is in addition to prediction. Our methods perform comparably with SoTA boosting, but also provide uncertainty quantification – we are not sacrificing one for the other.

The NeurIPS format does not allow space for an in-depth discussion of the uses of inferential quantities. We felt both that these would be broadly understood and that it was important to focus on the technical components of our results; but we will include some additional motivation in the introduction.

On the algorithmic changes to boosting and assumptions.

For the remainder of your comments; our results do require specific algorithmic changes to boosting. It is not clear how this can be done otherwise – for instance, boosting to a set number of rounds can lead to rampant overfitting and absolutely no hope for convergence or a CLT, while early stopping has its own issues, as we highlight in our response to Reviewer cjXq.

However, the assumptions about the data are not stronger than theoretical guarantees for most machine learning algorithms: we need the distributions of our features and response to behave reasonably, and we need the data we collect to be both independent and representative of the task. These assumptions are relevant for any ML modeling task. Regarding the assumptions about the trees and the learning algorithm itself, we address this in detail within Lines 131-137, 255-258, and 261-263. None of these assumptions are new – they are all inherited from prior work (Zhou and Hooker (2022), Athey et al. (2018)) on trees. As we allude to in Line 137, numerical experiments indicate that our boosting algorithms and inferential procedures work even when these assumptions are not explicitly verified.

We note that purely on an algorithmic level, our extension of Boulevard also provides an interpolation between boosting and bagging methods such as random forests, and while this was not our focus, it does allow us extra degrees of freedom to optimize between them.

At the end of the day, we do not think this review is fair. We strongly feel that both the reviewer's score, and the lion’s share of their objections, are unreasonable, and that this review does not provide a meaningful assessment of our work.

Given that you do acknowledge that our results appear to be correct, we hope that you will consider significantly revising your score; inadequacies in expressing motivation does not render a paper irredeemable.

Dear Reviewer oGZf,

We appreciate your support, and thank you for doing so!

Thank you for your kind review! We are heartened to see your kind comments on our paper being “very well written” (which is high praise!), that we offer both practical and theoretical advancements, and of you highlighting the utility of a formal justification for the use of dropout within boosting. We provide answers to your questions below (the second is addressed before the first, given that the second was also raised as a weakness).

Weaknesses and Questions

1. Practical guidance for choosing hyperparameters.

• Learning rate $\lambda$. We found that a learning rate close to $\lambda=1$ seemed to work well empirically. At no point did we choose a learning rate smaller than $0.1$. As a strategy, we would try $\lambda=1$ first, and decrease this in increments of $0.2$ based on cross-validated performance or its performance in a hold-out validation set. $\lambda$ should not be set to more than $1$, as specified in the setup of the paper in Line 99. Doing so invalidates the theory, and leads to instability in practice.
• Dropout probability $p$. Regarding the choice of dropout, we found that intermediate values of $p=0.3, 0.6, 0.9$ seemed to work well. The extremes lead to interesting behavior – $p=0$ (yielding vanilla Boulevard) did not yield stable prediction intervals, and $p=1$ yielding a random forest. In Lines 339-341, we highlight that the variance of the estimator, and therefore the CIs, grows as $\lambda$ and $p$ decrease – but we certainly agree that more discussion is warranted, and will make this amendment in the final paper.
• Number of boosting rounds $K$. We would set this to the highest possible $K$ that does not lead to a significant slowdown in the time taken for one (parallel) boosting round. This is likely to be a multiple of the number of processors one has available. Why? Having a larger $K$ essentially allows us to boost a mini-random forest. Boosting random forests (a strong learner) should lead to better prediction error than boosting a tree (a weak learner), if one keeps the number of boosting rounds the same. Algorithm 2 can in theory achieve an unbounded improvement in asymptotic relative efficiency compared to vanilla Boulevard. However, this improvement is in terms of the number of boosting rounds, not ensemble size!
  • Our results in Figure 2 regarding the MSE of Algorithm 2 are somewhat unfair to Algorithm 2 – as the x-axis is in the ensemble size, not the number of boosting rounds. The number of boosting rounds that Algorithm 2 encounters can be as few as 31 on Abalone and 50 on Wine Quality, in contrast to the 500 boosting rounds all other algorithms enjoy! Given the same number of boosting rounds, we would certainly expect an improvement.

2. Visual diagnostics on empirical normality. Recent response guidelines do not allow us to provide new figures. However, we have created the requested QQ plot from a simulation using Algorithm 1 and have lodged it with the AC.

• Qualitatively it gives us a good fit to the proposed CLT – this is almost a linear fit, with a slope that is just smaller than one. There is some indication of heavy tails in the extreme quantiles that involve only six out of 100 simulations, and the worst-behaved datapoint had a theoretical quantile of around $-2.75$ and a sample quantile of $-6$. This is not entirely unexpected – trees are notoriously bad at the very edge of the support of the data.
• In case this is inaccessible to the reviewer, we also conducted normality tests (Shapiro-Wilk and Jarque-Bera tests) based on the simulation. Results suggest that aside from outliers (that can be blamed on the poor behavior of trees near the edge of the data support), overall the predictions are pretty normal:

Overall, the results are encouraging for verifying empirical normality, and we thank the reviewer for this suggestion. We will include this QQ plot in the revised manuscript.

We hope that if our response and additional figures that we have generated in response to your review and that of Reviewers cjXq and AVWH have improved the paper in your eyes, that you might consider increasing your confidence in your assessment of the paper. Thank you for your kind review again!

Dear Reviewer kHmv,

Of course, checking for normality to see if the theory indeed holds in practice is probably something we should have done anyway in the original paper. Thank you for your kind review again!

We thank the reviewer for their kind comments on the strengths of the paper! We are very pleased to hear that you think that this is an important contribution towards understanding boosting algorithms. We are also glad that you think that the paper is structured clearly – we put a lot of thought on how this can be presented in 8 pages. We also thank you for your kind comments on how you think that our theory is involved and impressive.

Below are responses and answers to the weaknesses and questions you raised. These include more numerical results that pertain to whether or not early stopping affects our algorithms. We hope that these assuage your concerns.

Weaknesses and Questions

1. Are there other methods available for frequentist inference for gradient boosting? This is a good question. We hope that there are! Alternatives to Boulevard regularization are very much needed. We thank the reviewer for their comment – we had read Luo et al. (2016), but it did not come to mind that statistical inference could be possible with their approach. There are two things that you might mean, and we address them below.

• Do you mean that one can extract valid intervals out of their post-L2 boosting estimator, given that it is a valid post-selection inference algorithm? This is likely! Though we aren’t sure that this would provide inference for boosting in the sense that we had in mind. Their method performs boosting on component-wise linear models, and not trees.
  • Furthermore, their post-L2 boosting method refits a linear model on the predictors selected by their boosting procedure. This does not provide statistical inference for the model that was fitted by boosting.
  • This observation does raise a question. Viewing trees as OLS on adaptively selected rectangular indicator variables, can we fit a kernel regression on the kernel given by gradient boosting and perform inference on that model? Yes. This would constitute a similar algorithm to Friedberg et al. (2018), except applied to boosting and not random forests. We are working on it right now for another related paper. Please don’t scoop us!
• Do you mean that one may be able to obtain valid intervals in general given rates of convergence? Luo et al. (2016) show rates of convergence, but neither provide a CLT, nor do they use the convergence rates established to provide methods for statistical inference.
  • A CLT alone is not sufficient for statistical inference – one needs also an estimate of the asymptotic variance. It is unclear how this can be done in the framework of Luo et al. (2016). An analysis of the convergence rates is also not sufficient for statistical inference – constants are untracked, and the limiting distribution, or even whether asymptotic normality holds, is unknown.
  • Their paper makes no mention of confidence etc. intervals, and it is unclear to us how their results can be used for that. We think that this may be possible if one can either estimate the complexity of the early-stopped function class on the fly, or if one can track down the exact constants within their analysis -– but this extends far beyond the scope of their paper.
  • Having said that, we do agree that we should have included this within the literature review – we did read the paper, and found it to yield spiritually similar conclusions to Zhang and Yu (2005), which should also have been included. We thank the reviewer for pointing this out, apologize for the omission, and will fix this in the final version of the paper.

2. Stopping times. On questions about stopping criteria: one of the features of the Boulevard framework is that we do not need an explicit stopping criterion to avoid overfitting. Because it works on a running average rather than a sum of trees, it is akin to random forests where building a larger ensemble never hurts you, but leads to diminishing returns.

• In that sense, we do not need explicit stopping rules in our theory – in fact the case we analyse is based on the limit of infinite boosting rounds – and in practice would use similar tools as in random forests; monitoring test set performance and stopping at the point that we no longer see improvements.
• Having said that, the reviewer alludes to an important point: gradient boosting with early stopping is a promising direction for statistical inference. Certainly, if rates of convergence are available (such as with the work of Luo et al. (2016) and Zhang and Yu (2005)), perhaps something can be pulled out from them. Asymptotic normality may not hold, but non-normal limits and nonasymptotic results are still possible even then. We are interested in addressing this in future work.
• How does early stopping affect the empirical results? Last minute guidance restricts us from providing figures in this response, but we created a new figure that plots the test set performance against the number of boosting rounds which verifies this. We have lodged it with the AC. Under a moderate noise-to-signal ratio DGP with the ground truth being a sinusoid: $y = \sin(4\pi x) + \epsilon, \epsilon \sim \mathcal{N}(0,4)$, we fit Algorithm 1 with $p=0.5$ and Algorithm 2 with $K=2$, boost to 500 rounds, and keep all other hyperparameters the same.
  • The test MSE decreases rapidly for all models. For gradient boosting, and stochastic gradient boosting, it starts to increase dramatically after 10 to 20 rounds. XGBoost and LightGBM achieve the minimum at around the same mark, but the test MSE increases to a lesser degree (the former more than the latter) after this point. We attach a table below, depicting the results from highest final MSE to lowest final MSE at the end of 500 boosting rounds.
  • Algorithm 1, Algorithm 2, and vanilla Boulevard achieve approximately the same minimum, but the test set MSE remains flat thereafter – unlike the previously described boosting algorithms, it does not increase in the number of boosting rounds past this point. Early stopping is not necessary for our algorithms, both in theory and in practice!

3. Random forests. A final weakness was suggested to be the similarity structure between random forest and our parallelization structure.

• The referee could be referring to the contrast between the parallel structure of parallelized boosting (which starts to look more like a random forest) and the sequential processing in gradient boosting. This applies only to Algorithm 2 – as Algorithm 1 is a bona-fide sequentially processed gradient boosting algorithm for the L2 loss.
• In either case, this does not mean that Algorithms 1 and 2 are a reduction of boosting to random forests. While trees are averaged instead of summed, these trees are still grown in sequential order, one after the other. This leads to a convergence to a kernel ridge regression $\widehat{\mathbf{y}}_n^{(b)} \xrightarrow[b \rightarrow \infty]{a . s .}\left(\lambda^{-1}\mathbf{I}+(1-p) \mathbb{E}\left[\mathbf{S}_n\right]\right)^{-1} \mathbb{E}\left[\mathbf{S}_n\right] \mathbf{y}_n$, rather than convergence to a random forest $\mathbb{E}\left[\mathbf{S}_n\right] \mathbf{y}_n$. The analysis in terms of kernel ridge regression is in fact quite different, and these lead to different estimators in practice as well.

Again, we thank the referee for your kind comments and insightful questions. We hope that the response above is helpful and has addressed your concerns. If you think our response has adequately addressed your concerns, we will be glad if you choose to increase your score!

References

1. Tong Zhang and Bin Yu (2005), Boosting with early stopping: Convergence and consistency.
2. Friedberg et al. (2018), Local linear forests
3. Luo et al. (2016), High-dimensional L2 boosting: Rate of convergence

Dear Reviewer cjXq,

Thank you for the compliment! We appreciate the remark. Is your objection solely due to the fact that we are happy to boost for infinite rounds and do not need early stopping? If so, then the following points may be helpful as a point of clarification:

1. This is something inherited from Zhou and Hooker (2022), who construct the original Boulevard procedure that shares the same willingness to boost for infinite rounds.
2. Within Algorithms 1 and 2, trees are grown dependent on the output of previous trees, where we fit new trees on the residuals of ensembles of old trees. This is faithful to the boosting procedure, where the steps are dependent on each other. The only caveat is that the procedure does admittedly require eventual non-adaptivity, but this only needs to happen in the limit.
3. Stopping times are only one of the ways that we can ensure that boosting algorithms do not overfit. Alternatives include restricting tree depth, reducing the learning rate, employing L2 regularization, or even including dropout (e.g. DART). Boulevard regularization just so happens to be one way this can be done that is especially amenable to statistical inference. With this, everything remains the same about the original boosting procedure -- we simply update $f_{b+1} = \frac{b-1}{b}f_b + \frac{\lambda}{b} t_{b+1}$ instead of $f_{b+1} = f_b + \lambda t_{b+1}$. This can be thought of as a decreasing learning rate and ensemble shrinkage.

If your objection lies in the fact that this does require a non-standard modification to the boosting algorithm from what is commonly used in practice, that is a fair objection. We aren't too happy with it either, and are currently exploring alternatives in ongoing projects that amount to future work.

Again, if you think our response has adequately addressed your concerns, we will be glad if you choose to increase your score! Otherwise, please let us know if you do still have objections -- your feedback will be very helpful in improving the paper and how we present our results.