Repository: scikit-learn/scikit-learn
Pull request: BUG Fix zero division error in GBDTs
Fixes #14018

Probably also #14014

Non-regression test fails on master.

Problem statement:
Zero division error in HistGradientBoosting
```python
from sklearn.datasets import fetch_openml
from sklearn.model_selection import cross_val_score
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

import numpy as np

# one hundred plants - margin
bunch = fetch_openml(data_id=1491)
X = bunch.data
y = bunch.target


res = cross_val_score(HistGradientBoostingClassifier(max_iter=100, min_samples_leaf=5), X, y)
np.mean(res)
```
NaN

This dataset is a bit weird in that it has 100 classes with 16 samples each. The default parameter don't work very well but we should fail more gacefully.

cc @NicolasHug

Hints from issue discussion:
I am just adding the traceback


```pytb
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-3-b0953fbb1d6e> in <module>
----> 1 clf.fit(X, y)

~/Documents/code/toolbox/scikit-learn/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py in fit(self, X, y)
    247                     min_samples_leaf=self.min_samples_leaf,
    248                     l2_regularization=self.l2_regularization,
--> 249                     shrinkage=self.learning_rate)
    250                 grower.grow()
    251 

~/Documents/code/toolbox/scikit-learn/sklearn/ensemble/_hist_gradient_boosting/grower.py in __init__(self, X_binned, gradients, hessians, max_leaf_nodes, max_depth, min_samples_leaf, min_gain_to_split, max_bins, actual_n_bins, l2_regularization, min_hessian_to_split, shrinkage)
    195         self.total_compute_hist_time = 0.  # time spent computing histograms
    196         self.total_apply_split_time = 0.  # time spent splitting nodes
--> 197         self._intilialize_root(gradients, hessians, hessians_are_constant)
    198         self.n_nodes = 1
    199 

~/Documents/code/toolbox/scikit-learn/sklearn/ensemble/_hist_gradient_boosting/grower.py in _intilialize_root(self, gradients, hessians, hessians_are_constant)
    260             return
    261         if sum_hessians < self.splitter.min_hessian_to_split:
--> 262             self._finalize_leaf(self.root)
    263             return
    264 

~/Documents/code/toolbox/scikit-learn/sklearn/ensemble/_hist_gradient_boosting/grower.py in _finalize_leaf(self, node)
    399         """
    400         node.value = -self.shrinkage * node.sum_gradients / (
--> 401             node.sum_hessians + self.splitter.l2_regularization)
    402         self.finalized_leaves.append(node)
    403 

ZeroDivisionError: float division by zero
```
At a glance, the softmax is bringing probabilities close enough to zero, which causes hessians to be zero for the cross entropy loss.
I think the right fix is just to add a small epsilon to the denominator to avoid the zero division.

As Thomas noted these cases happen when the trees are overly confident in their predictions and will always predict a probability of 1 or 0, leading to stationary gradients and zero hessians.

We hit this part of the code in `initialize_root()`:
```py
        if sum_hessians < self.splitter.min_hessian_to_split:
            self._finalize_leaf(self.root)
            return
```

and since hessians are zero `finalize_leaf()` fails.

That's because the learning rate is too high BTW.
Changing the learning rate to .05 I get a .67 accuracy (not bad over 100 classes).


Will submit a PR soon

Commit message: added comment

Patch to review:
diff --git a/sklearn/ensemble/_hist_gradient_boosting/grower.py b/sklearn/ensemble/_hist_gradient_boosting/grower.py
index 691d10f16ac61..a6fbf61772e22 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/grower.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/grower.py
@@ -18,6 +18,9 @@
 from .types import PREDICTOR_RECORD_DTYPE
 
 
+EPS = 1e-15  # to avoid zero division errors
+
+
 class TreeNode:
     """Tree Node class used in TreeGrower.
 
@@ -398,7 +401,7 @@ def _finalize_leaf(self, node):
         https://arxiv.org/abs/1603.02754
         """
         node.value = -self.shrinkage * node.sum_gradients / (
-            node.sum_hessians + self.splitter.l2_regularization)
+            node.sum_hessians + self.splitter.l2_regularization + EPS)
         self.finalized_leaves.append(node)
 
     def _finalize_splittable_nodes(self):
diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
index f293402d8d650..39a43c9bb2762 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
@@ -172,3 +172,20 @@ def test_binning_train_validation_are_separated():
                   int((1 - validation_fraction) * n_samples))
     assert np.all(mapper_training_data.actual_n_bins_ !=
                   mapper_whole_data.actual_n_bins_)
+
+
+@pytest.mark.parametrize('data', [
+    make_classification(random_state=0, n_classes=2),
+    make_classification(random_state=0, n_classes=3, n_informative=3)
+])
+def test_zero_division_hessians(data):
+    # non regression test for issue #14018
+    # make sure we avoid zero division errors when computing the leaves values.
+
+    # If the learning rate is too high, the raw predictions are bad and will
+    # saturate the softmax (or sigmoid in binary classif). This leads to
+    # probabilities being exactly 0 or 1, gradients being constant, and
+    # hessians being zero.
+    X, y = data
+    gb = HistGradientBoostingClassifier(learning_rate=100, max_iter=10)
+    gb.fit(X, y)
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index 3b09ec2874997..0ee4cfefbdcde 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -2352,8 +2352,11 @@ def check_decision_proba_consistency(name, estimator_orig):
             hasattr(estimator, "predict_proba")):
 
         estimator.fit(X, y)
-        a = estimator.predict_proba(X_test)[:, 1]
-        b = estimator.decision_function(X_test)
+        # Since the link function from decision_function() to predict_proba()
+        # sometimes not precise enough (typically expit), we round to the 10th
+        # decimal to avoid numerical issues.
+        a = estimator.predict_proba(X_test)[:, 1].round(decimals=10)
+        b = estimator.decision_function(X_test).round(decimals=10)
         assert_array_equal(rankdata(a), rankdata(b))
 
 


Code review comments:
diff_hunk: @@ -18,6 +18,9 @@
 from .types import PREDICTOR_RECORD_DTYPE
 
 
+EPS = 1e-15  # to avoid zero division errors
path: sklearn/ensemble/_hist_gradient_boosting/grower.py
text: @user1:
would it make sense to use numpy values?

```python
np.finfo(np.float64).eps
```