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:
- 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.
- 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
- 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?)