FIX _safe_divide should handle zero-division with numpy scalar

[glemaitre]

Fixing the main branch.

This PR should handle the case of zero-division with two numpy scalar.

ping @lorentzenchr @thomasjpfan @OmarManzoor @lesteve

Edit: The CI failure popped up after merging #26278 in the example examples/ensemble/plot_gradient_boosting_regularization.py

[github-actions]

✔️ Linting Passed

All linting checks passed. Your pull request is in excellent shape! ☀️

_{Generated for commit: ad04d77. Link to the linter CI: here}

[glemaitre]

Given the failure in the example, it seems that the provided fix is not enough. We can hand up returning a very large value because the denominator is extremely small

numerator=-0.2857142857142857
denominator=3.18528276243611e-217
numerator=-8.969824879715574e+215

The next iteration, the numerator and denominator will nan.

[glemaitre]

I am not sure what it the right fix here to handle the numerical instability: either clip in the safe division because we should not assign such a large value in raw_predictions or handle the nan values of the next iteration?

[glemaitre]

Apparently, once we get some nan, it seems that this is too late:

vs. previously

[glemaitre]

I also see that we had potentially another strategy in the past: https://github.com/scikit-learn/scikit-learn/pull/26278/files#diff-dac0eef4868535102a7b9aeb319e1501dfbd10f1256f91d532927e12d30bfb15L732

[glemaitre]

So this is not enough. The negative gradient vector has some nan.

In the previous loss it seems that np.nan_to_num was used. Here I tried but the algorithm is still diverging.

[lorentzenchr]

I‘ll also have a look, but I need some time. Fortunately, we‘ve plenty of time to fix it.

[lorentzenchr]

Which version of numpy do you use?
With numpy 1.24.1 on macos, I can't reproduce this failure. My blue curve is like the one at the bottom plot of #27312 (comment).
I was mistaken, now I see the error.

[OmarManzoor]

@lorentzenchr I have numpy 1.25.1 on mac M1 and I can reproduce the error when I use the main branch.

[lesteve]

As an aside (I don't know how much we care to be perfectly honest): code based on np.seterr + catching an exception is not going to work in Pyodide, i.e. an exception will never be raised. Floating point exceptions are not supported in WebAssembly and this is unlikely to change in the near future, see numpy/numpy#21895 (comment) for more details.

[OmarManzoor]

@lorentzenchr @glemaitre @lesteve

scikit-learn/sklearn/ensemble/_gb.py

Lines 237 to 244 in c634b8a

	for leaf in np.nonzero(tree.children_left == TREE_LEAF)[0]:
	indices = np.nonzero(terminal_regions == leaf)[0] # of terminal regions
	y_ = y.take(indices, axis=0)
	sw = None if sample_weight is None else sample_weight[indices]
	update = compute_update(y_, indices, neg_gradient, raw_prediction, k)
	# TODO: Multiply here by learning rate instead of everywhere else.
	tree.value[leaf, 0, 0] = update

I think we should have

indices = np.nonzero(masked_terminal_regions == leaf)[0]  # of terminal regions

if we compare with the original code.

[OmarManzoor]

Plot if I use on main

indices = np.nonzero(masked_terminal_regions == leaf)[0]

[glemaitre]

Yep this look much what we had. It makes sense also if we where computing gradient on data that we should not have :)

[glemaitre]

@lesteve @lorentzenchr Instead of the exception catching, we could use the previous trick that return 0.0 if smaller than smaller than a really small values. This would be OK for Pyodide?

[lesteve]

we could use the previous trick that return 0.0 if smaller than smaller than a really small values. This would be OK for Pyodide?

I think this would work with Pyodide.

I am not too sure to which extent we want to support Pyodide quirkiness, as I mentioned above. It feels like not getting some warnings in Pyodide would acceptable but having a gradient boosting algorithm that behave weirdly because the loss becomes NaN or inf is maybe not that great.

[glemaitre]

So I made the fix of @OmarManzoor and change the error catching as it was previously.

[glemaitre]

I also added a comment to remember why we are not using np.errstate and who knows, we could potentially change it in the future if Pyodide handles it.

[lesteve]

I also added a comment to remember why we are not using np.errstate and who knows, we could potentially change it in the future if Pyodide handles it.

Don't hold your breath too much though, the link above says there is no plan to support it right now, so that's at least 3 years away.

[lorentzenchr]

LGTM just some nits.

@glemaitre @OmarManzoor @lesteve Thanks for fixing my bugs.

[OmarManzoor]

LGTM. Thanks @glemaitre , @lesteve and @lorentzenchr