BUG Fix zero division error in GBDTs (scikit-learn#14024)

Author: NicolasHug

sklearn/ensemble/_hist_gradient_boosting/grower.py

@@ -16,6 +16,10 @@
 from .predictor import TreePredictor
 from .utils import sum_parallel
 from .types import PREDICTOR_RECORD_DTYPE
+from .types import Y_DTYPE
+
+
+EPS = np.finfo(Y_DTYPE).eps # to avoid zero division errors
 
 
 class TreeNode:
@@ -398,7 +402,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):

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)
+], ids=['binary_crossentropy', 'categorical_crossentropy'])
+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)

sklearn/utils/estimator_checks.py

@@ -2401,8 +2401,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()
+ # is 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))