Is there an implementation of xgboost for a single target variable but using multiple regression parameters

Question

Is it possible, using custom cost functions or otherwise, to run xgboost regression for a single target variable, but rather than outputting just the conditional mean (conditioned on the feature values) of the variable as a prediction, to output condition distribution parameters.

For example, the way to interpret Poisson regression is that if xgboost predicts 0.67 based on your feature values, that you would expect the true value of your target variable to be distributed according to a Poisson distribution with $\mu=0.67$. In this case, that is also the expected value which in some sense masks some of the richness of what is going on under the hood.

But what if instead, we wanted xgboost to take our features and output $(p(x),r(x))$ and we would then interpret this as our target variable's true value being negative binomially distributed according to $f(y;p(x),r(x))$

The built-in cost functions include gamma which allows the user to do "gamma-regression". But when you do a gamma regression an then hit reg.predict(x), you only get 1 value, rather than two as one might expect (see here for more details)

As a reminder (see docs), for every partitioning of the space, the (albeit approximate) optimal value of the objective function at iteration t implied by said partitioning is given by (where the sum runs over all terminal nodes):

$$obj_{t}^{*}=-\frac{1}{2}\sum_{j=1}^{T}\frac{G_{j}^{2}}{H_{j}+\lambda} + \gamma$$

where

$G_{j} = \sum_{i \in \mathbf{R}_{j}}\frac{\partial \ell}{\partial \hat{y}}|_{y_{i}, f^{t-1}(x_{i})} $

and

$H_{j} = \sum_{i \in \mathbf{R}_{j}}\frac{\partial ^{2}\ell}{\partial \hat{y}^{2}}|_{y_{i}, f^{t-1}(x_{i})}$

Once the optimal partitioning has been found, the update associated with each node in the partitioning must be calculated and it's given by

$$w_{j}^{*}=-\frac{G_{j}}{H_{j} + \lambda}$$

If instead, we wanted our boosted tree to produce two outputs (let's call these $f(x)$ and $g(x)$, and when we take derivatives wrt them, where before we took the derivative wrt $\hat{y}$ we now take derivates wrt $(\theta_{1},\theta_{2})$), a use-case for this would be heteroskedastic regression where we want for any given $x$ not just to output an expected value but also a variance, I have managed to derive very similar update equations.

For a given partitioning of the space, the optimal objective is given by

$$obj_{t}^{*}=\gamma T-\frac{1}{2}\sum_{j=1}^{Y}G_{j}^{T}\cdot H_{j}^{-1}\cdot G_{j}$$

in which

$G_{j}= \begin{pmatrix}\sum_{i\in \mathbf{R}_{j}}\frac{\partial \ell}{\partial \theta_{1}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})} \\ \sum_{i\in \mathbf{R}_{j}}\frac{\partial \ell}{\partial \theta_{2}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})}\end{pmatrix}$

and

$H_{j}=\begin{pmatrix}\lambda + \sum_{i\in \mathbf{R}_{j}}\frac{\partial^{2}\ell}{\partial \theta_{1}^{2}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})} & \sum_{i\in \mathbf{R}_{j}}\frac{\partial ^{2} \ell}{\partial \theta_{1}\partial \theta_{2}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})} \\ \sum_{i\in \mathbf{R}_{j}}\frac{\partial ^{2} \ell}{\partial \theta_{1}\partial \theta_{2}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})} & \lambda + \sum_{i\in \mathbf{R}_{j}}\frac{\partial^{2}\ell}{\partial \theta_{2}^{2}}|_{y_{i}, f^{t-1}(x_{i}), g^{t-1}(x_{i})}\end{pmatrix}$

And once this optimal partitioning has been found, the optimal outputs associated with each node are given by

$$\begin{pmatrix}\theta_{1j}^{*}\\\theta_{2j}^{*}\end{pmatrix} = -H_{j}^{-1}\cdot G_{j}$$

Now clearly this all simplifies if the Hessian is diagonal. It's instructive to think of the problem of using XGBoost to train a single tree to solve two tasks at once (multi-target classification), or multi-class xgboost, which both essentially boil down to this, the tree is outputting multiple numbers for a given node but the loss functions is a sum of terms each of which only depends on one of the outputs, thus no cross-derivatives and the problem simplifies.

So my question boils down to "is there a good reason this isn't implemented in xgboost?" ("this" being problems where are likelihood is parametrised by multiple parameters in a way that creates cross-derivatives). Is there a reason that when we go from one to two parameter boosting, despite the math theoretically working out nicely like above, it becomes computationally too expensive.

It's somewhat instructive to think about where the heavy lifting actually happens in xgboost. My understanding is that it's mostly in (in the one parameter case), evaluation to term for $obj_{t}^{*}$ for every possible partitioning of the space. If you know your $G_{j}$ and $H_{j}$ within said partition, you need to square the former and divide by the latter, so it's two floating point operations.

For two-parameter xgboost, if you did this naively you would be inverting a matrix. Now of course, you're only ever going to want to be doing this for relatively small matrices, there's unlikely going to be many use-cases for 7 parameter boosting. Sticking with the two parameter case, if you invert the matrix by hand and write out the equation for $obj_{t}^{*}$ in terms of the two first gradients and three second gradients (let's call these $(g_{1}, g_{2}, h_{1}, h_{2}, h_{3})$ for convenience), I derive an equation along the lines of (the exact details aren't really that important for now)

$\frac{1}{h_{1}h_{3} - h_{2}^{2}}\cdot\left(h_{3}f_{1}^{2} + h_{1}f_{2}^{3} - 2h_{2}f_{1}f_{2}\right)$

so we go from 2 floating point operations to 14 (I think?) to evaluate this term, somewhat expensive yes, but not totally prohibitive. And of course, the fair comparison is not to compare 2 to 14 but to compare 4 to 14, as it's pretty commonplace to train xgboost to simultaneously solve two classification problems (2-target xgb)

So to summarise, the differences between 1 and 2 parameter boosting are:

1. In 1-parameter boosting, you need to, at the start of every boosting round, calculate a vector of gradients and second gradients. In 2-parameter boosting, it's two vectors of first gradients, and 3 vectors of second gradients.
2. For every partitioning of data, in 1 parameter boosting you need to move the relevant points between two sums for two terms, in two parameter boosting it's for 5 terms
3. In single parameter boosting, the $obj_{t}^{*}$ equation which needs to be repeatedly evaluated requires two floating point operations, and it's significantly more (I think 14), but not prohibitively so

Clearly this will scale badly in the number of boosting parameters, but it seems like most of the utility would be in expanding xgb from single parameter only to 2-3 parameter boosting.

In conclusion, have I missed something? Are there other reasons why this might be computationally prohibitively expensive? Is it likely the case that somebody has worked through this before and decided it's going to be prohibitively slow, or is it possible that it's not implemented just because nobody has decided to work on it yet (or the people interested in these sorts of problems decided to go with ngboost instead?)

Answer 1

This is definitely desirable to have, but not trivial to get. The way XGBoost, LightGBM, CatBoost et al. are created is with a single output at a leaf node in mind. When there's extra distributional parameters (e.g. shape for accelerated failure time models), these are hyperparameters. There's a few options with the way the package is, but these all involve multiple model fits (e.g. quantile regression, but it's hard to enforce that the quantiles stay ordered if it's not from a single model, fitting for different hyperparameter values and then taking the distribution to be based on how well the different hyperparameters fit the data etc.).

In general, I think it's hard because of how boosting works: you fit a model, then fit another model to the residuals, again and again. While you tend to sample observations for each model you create, because of the way each model builds on the previous one, the left-out observations don't directly provide an estimate of out-of-bag error (like e.g. in random forest, which can more easily give prediction intervals).

There's been data science competitions that targeted such an output (from whatever model you wanted) such as this one (see also discussions like these). These offer sources of inspiration for how people approach this kind of problem.

Answer 2

I don't think this is immediately possible with XGBoost as you would have to write a multi-output / multi-parameter boosting variant of it. At best, XGBoost (and other usual boosting routines learners) are able to do multi-output predictions (for example estimating the parameters for a Gamma distribution) by having one model for each target and then putting meta-estimators on top.

That said, what you primarily want is NGboost and the concept of "natural gradients" (NGs). NGs allow us to take into account the geometry of the underlying probability space, they are implemented in the package ngboost. In a standard NGBoosts application, we use Maximum Likelihood Estimation as our scoring function to try and minimize the Kullback-Leibler (KL) divergence between the current distribution and the target distribution. NGBoost explores a divergence (scoring function) -specific statistical manifold where each point in the manifold corresponds to a probability distribution.

This is a significant change to standard point estimation and allows our booster to output a full probability distribution over the entire outcome space instead of a point estimate. We can then use the NGBoost provided pred_dist function such that we get estimated distributional parameters. With ngboost for example, for a Gamma distribution we estimate the $\alpha$ (shape) and $\beta$ (rate) parameters aiming to minimize the KL-divergence. In contrast, XGBoosts (and other boosters) if minimising a Gamma distributed target focus directly on the mean Gamma deviance (mean_gamma_deviance of the point estimates.)