# Ordinal log-loss in a multiclass classification in XGBoost?

I have a multi-class problem that which classes are simultaneously mutually exclusive and have ordering. You can think of the classes as being some score: 0 (Low), 1 (Medium), 2 (High).

What I would like to do is make a custom loss function that will capture this characteristic. Following from Section 7.2.1.

XGBoost's documentation gives the loss function $$L(\theta)$$ for a binary classification, which is simply the log-likelihood as given by:

$$L(\theta) = \sum_{i}[y_{i}ln(1+e^{-\hat{y}}) + (1-y_{i})ln(1-e^{-\hat{y}})].\\$$

Although not documented, I would assume that for multiclass with softmax:

$$L(\theta) = \sum_{i}[y_{i,1}ln(\frac{e^{\hat{y_{i,1}}}}{\sum_j^{k}e^{\hat{y_{i,j}}}}) + y_{i,2}ln(\frac{e^{\hat{y_{i,2}}}}{\sum_j^{k}e^{\hat{y_{i,j}}}}) + ... y_{i,m}ln(\frac{e^{\hat{y_{i,m}}}}{\sum_j^{k}e^{\hat{y_{i,j}}}})]$$

From here, two questions follow:

1. Since XGBoost does not provide a ordinal loss, I can only adjust the loss function after the softmax transformation. This is problematic because in the proportional log-odds (ordinal) problem we have:

$$P(y=1) = P(y^*<\tau_{1}) = P(\alpha + \beta X + \epsilon < \tau_{1}) = P(\epsilon < \tau_{1} - \alpha - \beta X ) = \frac{1}{1+e^{(\tau_{1} - \alpha - \beta X)}}$$

which in XGBoost's case will be:

$$\frac{1}{1+e^{(\tau_{1} - \hat{y})}}$$

Since I cannot invert the softmax function without a constant to get $$\hat{y}$$, this is problematic .. If you could do this then I would assume that you can then simply just sum the log-probabilities to get the log-likelihood as your desired loss.

1. In XGBoost, when you are coding the custom loss (see, code snippet), predt and dtrain are provided to you; if you input a multiclass vector suitable for ordinal classification, dtrain will be, e.g., [0,0,1,3,2], and predt will be a matrix of probabilities, which is the softmax-transformed $$\hat{y_{i,j}}$$ for each observation $$i$$ and class $$j$$. However, the ordinal problem, following from above, seems to be one-dimensional (it resembles more a binary classification than multinominal), so I am not fully understanding how can you work with this. I have checked out the following post, Custom Objective and Evaluation Metric, which explains some differences in encoding, but this affects the dtrain if anything and not the probabilities. I have also looked at the statsmodels.miscmodels.ordinal_model.OrderedModel implementation and they do accept a vector of integers (e.g., [0, 0, 1, 3, 0, 2]) as their exog variable, so I am assuming this is possible to implement.

I've cross posted this question in: StackExchange - Data Science - Ordinal Log Loss in a Multiclass Classification in XGBoost

• The canonical loss for ordinal data is the Ranked Probability Score, would that be an alternative? Commented Nov 15, 2023 at 12:21
• @Ben, is this not allowed? Commented Nov 15, 2023 at 13:58
• How can I link them? Simply state in both of them that the same post can be found on the other stackexchange? Commented Nov 15, 2023 at 14:55
• Although you can use the log-likelihood, I agree with Stephan Kolassa that the Ranked Probability score probably makes more sense here. Unfortunately, I don't understand your questions. Where is the problem? For instance, what do you mean by " the ordinal problem [...] seems to be one-dimensional [...], so I am not fully understanding how can you work with this. "? Commented Nov 16, 2023 at 18:09
• My understanding of ordinal loss is that it is a 1-to-3 mapping, unlike softmax which is 3-to-3 mapping. Therefore, this 1-to-3 mapping, I think, cannot be implemented in XGboost because it doesn't allow you to input a vector of more than 2 classes without providing predt as a matrix in the custom loss implementation. Unless, I am misunderstanding something and the ordinal loss can simply be rewritten as a softmax instead of a sigmoid. Is it more clear what my struggle is? Commented Nov 17, 2023 at 12:31

For a start, you these packages (whether we're talking xgboost or lightgbm or whatever) are almost certainly working with non-softmaxed logits. You can see that by e.g. getting them for a trained multi-class model (here for the python lightgbm package) via model.predict(mydata, raw_score=True). That also illustrates for you what happens under the hood: there's one tree per category that outputs a number that acts as the non-softmaxed logit for a category. You should be able to access those and apply whatever loss function you want to these, it's just a matter of finding the right bit of documentation.
However, if you want to work with an ordinal loss-function, you also have another option you work with a single continuous output $$\hat{\theta}$$ (i.e. just one tree per iteration instead of one-per-category, which is the default for unordered classes). You'd assume an underlying $$N(\hat{\theta}, 1)$$ variable (for which you predict the mean) and define cut-points $$c_\text{Low-Medium}$$ and $$c_\text{Medium-High}$$ that indicate the boundaries between the categories. You then define that $$P(\text{Low}|\hat{\theta}) = P(N(\hat{\theta}, 1), $$P(\text{Medium}|\hat{\theta}) = P(N(\hat{\theta}, 1) \in [c_\text{Low-Medium}, c_\text{Medium-High})$$, and $$P(\text{High}|\hat{\theta}) = P(N(\hat{\theta}, 1)\geq c_\text{Medium-High})$$. You can even fix one of the boundaries, e.g. $$c_\text{Low-Medium}:=0$$, while the other one would be a hyperparameter to tune. You may or may not find this easier to implement (perhaps at the cost of having additional hyperparameters to tune).
• It would be a vector of integers the class of a record. The loss from a record $i$ with class $j_i$ is $-\log \hat{p}_{j_i}$. We get $\hat{p}_{j_i}$ from the continuous model output $\hat{\theta}_i$ as $\hat{p}_{j_i} = \Phi(c_\text{upper}(j_i)| \hat{\theta}_i, 1) - \Phi(c_\text{lower}(j_i) | \hat{\theta}_i, 1)$, where $\Phi(x|\mu, \sigma)$ is the CDF of a normal distribution with mean $\mu$ and SD $\sigma$, and where these limits are as above (and the lower limit is defined as $-\infty$ for the first class, while the upper limit is defined as $\infty$ for the last class). Commented Nov 28, 2023 at 11:50
• I am unsure if I am following. In regular softmax, i.e., loss from record $i$, with class $j_{i}$ is $L(\hat{y}) = -log\hat{p}_{j_{i}}$. Since$\hat{p}_{j_{i}}$ is defined as $\frac{e^{\hat{y}_{j}}}{\sum_{k=1}^{K}e^{\hat{y}_{k}}}$, $y$, directly appears in the function ($\hat{y}$ being the model output, so, you can take the first derivative w.r.t to the model output $\frac{\partial L}{\partial\hat{y}}$. But, you have $\theta$, which is latent variable I believe (so, not $\hat{y}$ which is given by XGboost), and moreover, you will have to take the derivative w.r.t. a difference in CDFs? Commented Nov 29, 2023 at 10:22
• By $\hat{\theta}$ I did mean what you call $\hat{y}$ (i.e. directly the xgboost ouput), but didn't want to call it that bc $y$ often refers to the observable data (here class label) when $\theta$ is a latent variable. Yes, you'd ideally provide xgboost with analytical gradient + Hessian & that it may not be so trivial in this case (haven't tried it, you can also try other CDF defined on all reals - e.g. logistic dist). Alternatively, you could try automatic differentiation via e.g. the autograd library, but I would expect this to be (perhaps a lot) slower than an analytic solution. Commented Nov 29, 2023 at 10:56