Possible math error in Binomial deviance calculations

[ChrisBeaumont]

The calculations of the binomial deviance and its gradient in https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/ensemble/gradient_boosting.py are bothering me. I could be missing something, but:

The deviance is calculated as log(1 + exp(-2 * y * pred)) here. This matches equation 10.18 on p 346 of Elements of Statistical Learning. However, that derivation assumes that y is {-1, 1} valued, whereas in sklearn y is {0, 1}. Effectively, the calculation is insensitive to pred whenever y=0. The fix is to change the return line to np.sum(np.logaddexp(0.0, -2 * (2 * y - 1) * pred)) / y.shape[0].

The calculation of the gradient makes sense to me if the pred values map to class probabilities via P(y=1) = 1 / (1 + exp(-pred)). However, the loss function calculation above seems to follow the convention that P(y=1) = 1 / (1 + exp(-2 * pred)) (again, see the link above). One way to make the two equations consistent with each other is to remove the first 2 in the above equation:

np.sum(np.logaddexp(0.0, -(2 * y - 1) * pred)) / y.shape[0]

[pprett]

Hi Chris,

thanks for looking into this! You are correct - the deviance is indeed incorrect, however, as far as I can see the gradient is correct thus the bug might have only affected model selection based on deviance (as in example plot_gradient_boosting_regularization.py).
I used the derivative from G. Ridgeway, "Generalized Boosted Models, A guide to the gbm package", 2007 [1]

Actually, there is also a bug in the plot_gradient_boosting_regularization.py example. In https://github.com/scikit-learn/scikit-learn/blob/master/examples/ensemble/plot_gradient_boosting_regularization.py#L65 I pass y_test to BinomialDeviance.__call__ which contains labels in {-1, 1} instead of {0, 1}.

I'll submit a fix in a couple of minutes.

[1] http://cran.r-project.org/src/contrib/Archive/gbm/gbm_1.6-3.tar.gz

[pprett]

@ChrisBeaumont it would be great if you could review the above changes - thanks!

[amueller]

@pprett this was fixed in #1627, right?

[ChrisBeaumont]

Yes, I believe so

[pprett]

yes - we have an explicit test for that in test_gradient_boosting_loss_functions.py .

(still binomial deviance gives me a hard time lately... working on a PR )

[chuan92]

It seems the error not fixed in #1627, the binomialdeviance is still not consistent with definition of deviance from ESLII. In test_gradient_boosting_loss_functions.py, the test case only contains large positive or negative numbers. If you test

import numpy as np
from sklearn.ensemble.gradient_boosting import BinomialDeviance
from numpy.testing import assert_almost_equal

bd = BinomialDeviance(2)
datum = (np.array([1.0, 1.0, 1.0]), np.array([3.0, 3.0, 3.0]))
# check if same results as alternative definition of deviance (from ESLII)
alt_dev = lambda y, pred: np.mean(np.logaddexp(0.0, -2.0 * (2.0 * y - 1) * pred))
assert_almost_equal(bd(*datum), alt_dev(*datum))

, it will fail.

If we choose the definition of ESLII, then

class BinomialDeviance(ClassificationLossFunction):

def __call__(self, y, pred, sample_weight=None):
        """Compute the deviance (= 2 * negative log-likelihood). """
        # deviance = 2* neg log-likelihood ?  
        # logaddexp(0, v) == log(1.0 + exp(v))
        pred = pred.ravel()
        if sample_weight is None:
            return -np.mean(2*(y*pred)-np.logaddexp(0.0,2*pred))
            # not return -2.0 * np.mean((y * pred) - np.logaddexp(0.0, pred))

@pprett @amueller @ChrisBeaumont

[amueller]

@pprett can you have a look again?

[ColinSuess]

Repro'd. @chuan92 's test throws an assertion on 0.19.2 with Python 3.7.

[lorentzenchr]

After reading carefully, with the convention y = 0 or 1 (or even y in [0, 1]), I think the actual implementation of the binomial deviance, here is OK. It is missing a constant factor that guarantees deviance >= 0 and deviance=0 only if y=predict_proba (which is currently not possible in scikit-learn), but this is irrelevant for optimization, see Eq.(B.21) here.

On the other hand, the negative_gradient is missing a factor of 2, here.

More infos in the commits of PR #18250.