This forces the support of values to be in $[0, 1]$, but is this actually the appropriate implementation? Isn't the problem else where in the implementation?

I don't know. I was trying to understand what the lower and upper bound that are stored in the tree builder stack (and then passed to this function) are used for/how their value is computed but not sure.

The very first node that gets pushed to the stack uses -inf and +inf for the lower and upper bounds. Which would suggest it is possible for the values to be out side the range of zero to one?

You were right @jjerphan the first attempt was a quick and dirty fix. As suspected there was something else going on... Now with d74b012 we're closer to a real solution, although we still have to understand the remaining floating point precisions issues...

I think clipping this is a good hot fix for the support and we may first proceed with it.

But values are not supposed to be in [0, 1].

But values are not supposed to be in [0, 1].

Exactly ! I was not aware of this until your comment #27639 (comment) @glemaitre.

OK so the bug is here. We don't expect to have a sum at 1 but a weighted count.
So we should have something like the sum_total - dest[0] instead.
I don't know if we have access to sum_total thought at this stage.

Uhm actually this is not even sum_total because this is grouped by labels.

So with my limited knowledge of the wider implementation of the constraints, I would come with the following code:

        cdef float64_t total_weighted_count = dest[0] + dest[1]

        if dest[0] < lower_bound:
            dest[0] = lower_bound
        elif dest[0] > upper_bound:
            dest[0] = upper_bound

        dest[1] = total_weighted_count - dest[0]

Thanks a lot for looking into this @glemaitre !

We don't expect to have a sum at 1 but a weighted count

I was not aware of this obviously, so this is a big mistake in the monotonic constraints implementation indeed !

Uhm actually this is not even sum_total because this is grouped by labels.

Monotonic constraints are only supported with single label, so this should not be an issue.

Note that lower_bound and upper_bound are positive class frequencies so we would need to modify your solution to

        cdef float64_t total_weighted_count = dest[0] + dest[1]

        if dest[0] < lower_bound * total_weighted_count:
            dest[0] = lower_bound * total_weighted_count
        elif dest[0] > upper_bound * total_weighted_count:
            dest[0] = upper_bound * total_weighted_count

        dest[1] = total_weighted_count - dest[0]

The patch almost works !

test_monotonic_constraints_classifications currently fails due to floating point arithmetic fluctuations. The tests pass if the predicted probabilities are cast to np.float32 so the overall logic is valid, which is good news.

I'm looking into the float64 precision matter.

Nice. Now we only have to deal with the 25 CIs with some 32-bits arch ;)

If we have some casting in float32, we can probably play with the relative tolerance.

My understanding is that the numerical noise comes from the predictions on clipped leaves:

(lower_bound * total_weighted_count) / total_weighted_count

which do vary with total_weighted_count (in floating point arithmetic) close to lower_bound. So two leaves can have values (before clipping) trespassing the same constraints, say lower_bound, but with different total_weighted_count and their predictions fluctuate a bit and sometimes break the monotonicity.

I found a solution that works by loosing precision on lower_bound (and upper_bound) before performing the scaling:

        cdef float64_t total_weighted_count = dest[0] + dest[1]
        cdef float64_t scaled_lower_bound = (<float32_t>lower_bound) * total_weighted_count
        cdef float64_t scaled_upper_bound = (<float32_t>upper_bound) * total_weighted_count

        if dest[0] < scaled_lower_bound:
            dest[0] = scaled_lower_bound
        elif dest[0] > scaled_upper_bound:
            dest[0] = scaled_upper_bound

        dest[1] = total_weighted_count - dest[0]

However, it does have limitations:

• sample_weights that have a precision of float32 or more (tested float16 and this is OK)
• n_samples too close to 2**32, so larger than a few hundred of millions (not tested for n_samples beyond a few millions though...)

This is probably the best that we can do. I don't know if we should have a check on the sample weight then in this specific setting.

To void early loss of precision.
Why do we need to cast to float32 at all?

To void early loss of precision.
Why do we need to cast to float32 at all?

The precision loss for lower_bound and upper_bound (and not for their scaled counterpart) is by design.

The subtlety here is that for leaves affected by clipping (say constrained by lower_bound), the predicted probabilities (without casting) would be (lower_bound * total_weighted_count) / total_weighted_count, which varies with total_weighted_count close to lower_bound.
So if two leaves end up having the same lower_bound desired prediction, but two different total_weighted_count, their predictions might end up being (slightly) different and then (slightly) breaking the monotonicity constraint.

By casting lower_bound to float32 we make ((<float32_t>lower_bound) * total_weighted_count) / total_weighted_count effectively constant for a pretty large set of total_weighted_count values, most interestingly integers < 2**30. The reason is that:

• the multiplication by total_weighted_count in this subset of values can be done without rounding
• the division by total_weighted_count (the same number) will be exact thus requiring no rounding either and equal to the initial <float32_t>lower_bound.

Please let me know if any of this is not clear enough !

Thanks for the explanation. It would be good to have a short summary of this as code comment.

Then, about the solution: By casting to float32, we loose 7 significant digits, while the loss of accuracy by a multiplication and a division is roughly 1 ulp. I would therefore prefer to have something like

float64_t scaled_lower_bound = nextafter(lower_bound, 1) * total_weighted_count
float64_t scaled_upper_bound = nextafter(upper_bound, 0) * total_weighted_count

Thanks for the suggestions. Adding an explanatory comment sounds like a good idea, not sure how short I can make it though...

Regarding the precision loss on lower_bound and upper_bound, we want the sequence of multiplication and division by the same number to be exactly neutral, which requires more loss of precision than a single bit.
As you can see

def check_nextafter(a, n): return (nextafter(a, 1) * n) / n == nextafter(a, 1)
check_nextafter(0.8793444460978181, 700)
>>> False

whereas

def check_castfloat32(a, n): return (np.float32(a) * n) / n == np.float32(a)

requires n > 2**29 for returning False, which seems good enough.

That's not a fair comparison and not the one we are interested in. We want scaled_lower_bound / total_weighted_count >= lower_bound, but as close as possible (that's where float32 sucks, in the end, it might not influence the estimator in a noticable way).

from math import nextafter
import numpy as np

def check_nextafter(a, n):
     """a between 0 and 1, n any number"""
     return (nextafter(a, 1) * n) / n >= a

sum([not check_nextafter(0.8793444460978181, n) for n in range(1, 1000_000)])  # 0
sum([not check_nextafter(0.8793444460978181, n) for n in np.logspace(0, 100,  1_000_000)])  # 0

Which assert in which test does fail?
I don't now of a "multi_nextafter", but we could just use lower_bound * (1 + x * np.finfo(float).eps) and lower_bound * (1 - x * np.finfo(float).epsneg) for some x > 1.

for some x > 1

If this is for the steps, the nextafter from python as a steps parameter.

Which assert in which test does fail?

In test_monotonic_constraints_classifications:

assert np.all(est.predict_proba(X_test_0incr)[:, 1] >= proba_test[:, 1])

the nextafter from python as a steps parameter.

Indeed, this would be useful, although it was only added in Python 3.12...

Arff my bad. I did not check the full documentation. Python 3.12 is quite ahead.

Do we need to protect dest[1] to be non-negative?

lower_bound and upper_bound are computed as positive class frequencies (or averages of positive class frequencies), so they should always be in the [0, 1] range, and thus dest[1] should always be non-negative.

If we decide to protect as you suggest though, note that we should make sure the invariant dest[0] + dest[1] = total_weighted_count holds.

Thanks for the explanation. It would be good to have a short summary of this as code comment.

Then, about the solution: By casting to float32, we loose 7 significant digits, while the loss of accuracy by a multiplication and a division is roughly 1 ulp. I would therefore prefer to have something like

float64_t scaled_lower_bound = nextafter(lower_bound, 1) * total_weighted_count
float64_t scaled_upper_bound = nextafter(upper_bound, 0) * total_weighted_count

Do those tests fail without this PR? (They should!)

    assert np.logical_and(
        proba_test >= 0.0, proba_test <= 1.0
    ).all(), "Probability should always be in [0, 1] range."

definitely fails.

    assert_allclose(proba_test.sum(axis=1), 1.0)

doesn't though as we made sure the sum is 1..

This is indeed much cleaner.