# Loss function in for gamma objective function in regression in XGBoost?

Suppose I want to predict $$y$$ from a set of predictors $$x$$. $$y$$ is gamma distributed, so I want to use gamma regression with XGBoost. The help page of XGBoost specifies, for the objective parameter (loss function):

reg:gamma: gamma regression with log-link. Output is a mean of gamma distribution. It might be useful, e.g., for modeling insurance claims severity, or for any outcome that might be gamma-distributed.

What is the explicit formula for this loss function? I haven't been able to find it. Also, Why is the log-link required? As far as I understand, log is usually used as the link function when using Generalized Linear Models, but in gradient boosting it would make any sense...

I'm not actually sure whether "gamma regression" is officially defined (it doesn't appear to have a wikipedia page for example), but if I were to define it (and some googling around suggests I'm not alone here), I would define it as setting up my regression problem so that for a given input vector $$\underline{x}$$, I predict a value y. When I do this, I meant that I think that the true value of the target will be gamma distributed with mean y.

How does this differ from least squares? One setup which leads to OLS is that you say that for given $$\underline{x}$$, the target variable will be normally distributed and your prediction $$y(\underline{x})$$ is the mean of that distribution. Of course a normal distribution is parametrised through its mean and variance, but it turns out that you don't need to know the variance in order to calculate the cost function you need to optimise, and thus this parameter doesn't need to be passed to xgboost.

For the gamma distribution however, this is different. Let's go through the maths. For the gamma distribution parametrised as $$\frac{1}{\Gamma (k)\theta ^{k}}x^{k-1}e^{-\frac{x}{\theta}}$$, the mean is given by $$k\theta$$ and the variance by $$k \theta ^{2}$$

Thus let's reparametrise in terms of $$\mu$$ and $$\theta$$ so that the distribution is given by $$\frac{1}{\Gamma (\frac{\mu}{\theta})\theta ^{\frac{\mu}{\theta}}}x^{\frac{\mu}{\theta}-1}e^{-\frac{x}{\theta}}$$

So for a given dataset, if you predict a bunch of $$\hat{y}_{i}$$ for target values $$y_{i}$$, the likelihood is given by

$$\prod _{i=1}^{N} \frac{1}{\Gamma (\frac{\hat{y}_{i}}{\theta})\theta ^{\frac{\hat{y}_{i}}{\theta}}}y_{i}^{\frac{\hat{y}_{i}}{\theta} -1}e^{-\frac{y_{i}}{\theta}}$$

and thus the negative (xgboost assumes a cost function you're trying to minimise) log-likelihood is

$$\sum _{i=1}^{N} \ln \Gamma (\frac{\hat{y}_{i}}{\theta}) + \frac{\hat{y}_{i}}{\theta}\ln \theta - (\frac{\hat{y}_{i}}{\theta}-1) \ln y_{i} + \frac{y_{i}}{\theta}$$

Compare this to the Gaussian Regression case where the negative log likelihood is given by $$\frac{1}{\sigma ^{2}}\sum _{i=1}^{N} \left(y_{i} - \hat{y}_{i}\right)^{2}$$

The the latter case, the $$\frac{1}{\sigma ^{2}}$$ is a constant term out the front. If you were doing linear regression or even xgboost without regularisation, this would mean that no matter what value you changed $$\sigma$$ to, the linear regressor/xgboost you trained would turn out to be exactly the same, so "Gaussian regression with $$\sigma = 10$$ and Gaussian regression with $$\sigma = 1$$ lead to the same predictions". This is no longer true when you have a regulariser, but you can always suck the value of $$\sigma$$ into the definition of the regulariser to get around this, and this is why the OLS formula never includes a $$\sigma$$ in it.

In the gamma case however, because of the $$\theta$$ factor contained in the $$\Gamma$$ function and the $$\ln \theta$$, you can't just pull the factor of $$\theta$$ outside of the summation.

For xgboost, you now need to pass it the elementwise first and second derivatives of the cost function wrt $$\hat{y}_{i}$$. This is where basic calculus doesn't get you all the way, you'll likely need to look up that the derivative of the logarithm of the gamma function is given by the digramma function $$\psi (z)$$.

The (elementwise) first derivative of the loss will be given by (by the xgboost definition which is $$G_{i}=\frac{\partial L}{\partial \hat{y}_{i}}$$):

$$G_{i} = \frac{1}{\theta}\psi (\frac{\hat{y}_{i}}{\theta}) + \frac{1}{\theta}\ln \theta - \frac{1}{\theta}\ln y_{i}$$

The second derivative will require derivatives of the digamma function, I don't know much about this but some googling tells me you need the trigamma function $$\psi _{1}(z)$$ which is the derivative of the digamma function, thus

$$H_{i}=\frac{1}{\theta ^{2}}\psi_{1}(\frac{\hat{y}_{i}}{\theta})$$

Again, note that you still have to supply $$\theta$$ up front as a hyperparameter, pass this to xgboost and then train a new xgboost model every time you wish to investigate another $$\theta$$

Finally, it's worth noting that I just did this derivation myself today, I haven't lifted anything other than the definition of the gamma distribution from elsewhere, so there could easily be minor algebra error, I'd feel more comfortable if somebody else independently verified my workings.

Edit: Alternately, you could parametrise the other way around:

You could use k as your free parameter and $$\theta = \frac{\mu}{k}$$, thus your gamma distribution is $$\frac{1}{\Gamma(k)(\frac{\mu}{k})^{k}}x^{k-1}e^{-\frac{xk}{\mu}}$$

and thus your negative log-likelihood is given by

$$L = \sum_{i=1}^{N}\left[\ln \Gamma (k) + k\ln \hat{y}_{i} - k\ln k - (k-1)\ln y_{i} + \frac{y_{i}k}{\hat{y}_{i}}\right]$$

this parametrisation is easier to differentiate, you get

$$G_{i}=\frac{k}{\hat{y}_{i}} - \frac{y_{i}k}{\hat{y}_{i}^{2}}$$

and

$$H_{i}=-\frac{k}{\hat{y}_{i}^{2}}+ 2\frac{y_{i}k}{\hat{y}_{i}^{3}}$$

but you still need to pass k as a hyperparameter

• Thanks for your answer. The derivation of the negative likelihood seems right but, in XGBoost you can definitely use objective reg::gamma and no $\theta$ parameter has to be inputed. There is also a gamma-nloglik eval_metric available, but I have not been able to obtain a formula for the first and second derivatives of this loss function.
– D1X
Aug 27, 2020 at 9:41
• My c++ isn't good enough to fully understand the source code (github.com/dmlc/xgboost/blob/…) but it looks to me like it's hard-coded to be 1. I would suggest using the learning API and coding this up yourself, that will allow you to pass $\theta$ manually. Last I used it, the syntax (in python) was awkward and the cost function didn't allow for other kwargs so you need to make theta a global variable. Aug 27, 2020 at 9:56
• Okay, I may be wrong (please somebody correct me if I am), but this is how I understand it. In regression we are interested in predicting $E(y \ | \ x) = \mu$. For a gamma distributed value, this would be $\mu = k \theta$, for $k, \theta > 0$ shape and scale respectively. Thus, for predicting $\mu$ we can just fix $\theta = 1$ and then adjust just the $k$, as this can adjust for any parameter of $\mu$. This means the distribution of $y \ | \ x$ will probably not be very reliable, but still $\mu$ will be, which suffices for deterministic prediction.
– D1X
Aug 27, 2020 at 10:49
• (continuation) Setting $\theta = 1$ also simplifies the log likelihood, which would become $\sum _{i=1}^{N} \ln \Gamma (\hat{y}_{i}) - (\hat{y}_{i}-1) \ln y_{i} + y_{i}$.
– D1X
Aug 27, 2020 at 11:02
• I think you could parametrise the other way around, let k be your free parameter, then your cost function is: $\sum_{i=1}^{N}\left[\ln \Gamma(k) + k \ln \hat{y}_{i} -k\ln k - (k-1)\ln y_{i} +\frac{y_{i}k}{\hat{y}_{i}}\right]$, but you still have to arbitrarily choose k, or run for different k values and see which gives you the best training loss. I think this parametrisation looks more like the source code and also doesn't require di- and tr-gamma functions. Will update my original post in due course to include this parametrisation. Aug 27, 2020 at 11:56

Thanks everybody for the contributions! I am late to the party, but I wanted to add one more point regarding the log-link function, which was to me still unclear.

I take the formula for the Gamma distribution from the bottom of gazza89's answer: $$\frac{1}{\Gamma(k)(\frac{\mu}{k})^k}x^{k-1}e^{-\frac{xk}{\mu}}$$ Using the logarithm as link function amounts to substituting $$\mu$$ with $$e^\mu$$ in this formula: $$\frac{1}{\Gamma(k)(\frac{e^\mu}{k})^k}x^{k-1}e^{-\frac{xk}{e^\mu}}$$ We can then compute the negative log likelihood, which looks like: $$L = \sum_{i=1}^N\left[\log\Gamma(k) + k\hat y -k\log k-(k-1)\log y +yke^{-\hat y} \right]$$ The gradient is then: $$G_i = k(1 - ye^{-\hat y})$$ And the Hessian: $$H_i = yke^{-\hat y}$$ These are the formulas provided in the source code (link):

          _out_gpair[_idx] = GradientPair((1 - y / expf(p)) * w, y / expf(p) * w);


with the notation change $$k=w$$ and $$\hat y = p$$.

• I agree with your working, this is consistent with what I said previously. What confuses me, is how k (or w) comes into this. The docs are pretty light on this, all I can find is (xgboost.readthedocs.io/en/stable/python/examples/…), but nowhere does it seem to let you set k. My point is that this isn't some arbitrary constant whose value doesn't matter. If you assume it to be 1 or 10, you are making very different assumptions as to the distribution of your target variable around your prediction May 19 at 10:43
• I know this is old, but I think it’s kinda important and this is the first page that anyone finds when they google gamma loss xgboost… The k is set by the weights. Since the log likelihood is linear in it, only ratios matter for the purpose of fitting. The distribution is obviously going to be different, but that’s the case for all exponential dispersion families (like the number of trials in binomial regression or the variance in normal regression). If a probabilistic prediction is required, the dispersion parameter has to be determined after fitting. May 26 at 10:09
• So let’s assume you’re modelling insurance claims or conversion values or something like that. If you have 5 observations with the same features (let’s say same ad campaign, same click and conversion date, etc.) and the value in total was 5000, then your target would be 1000 and your weight would be 5, in agreement with the formula for the distribution of the sum of independent Ga-distributed RVs. Note that any dataset can have weights in XGBoost or LightGBM. May 26 at 10:14