roc_auc_score computation is wrong for large samples

[arogozhnikov]

Hi, we've recently found some strange inconsistency in the behavior of roc_auc_score. After long digging @tata-antares found out that it is wrongly computed for large data samples if weights passed are float32.

Minimal example

from sklearn.metrics import roc_auc_score
import numpy
numpy.random.seed(42)
​
n_samples = 4 * 10 ** 7
y = numpy.random.randint(2, size=n_samples)
prediction = numpy.random.normal(size=n_samples) + y * 0.01
trivial_weight = numpy.ones(n_samples)

roc_auc_score(y, prediction)
#prints 0.50273924526046887

roc_auc_score(y, prediction, sample_weight=trivial_weight.astype('float32'))
#**prints nan** and wrong warning
/moosefs/miniconda/envs/ipython_py2/lib/python2.7/site-packages/sklearn/metrics/ranking.py:530: UndefinedMetricWarning: No negative samples in y_true, false positive value should be meaningless
  UndefinedMetricWarning)

roc_auc_score(y, prediction, sample_weight=trivial_weight.astype('float64'))
#prints 0.50273924526046887

The value when dtype=float64 is correct, but in the case we investigated after passing float32 metric returned a regular number, which is very different from correct one (e.g 0.58 instead of 0.52). It is some occasion we detected this problem.

sklearn == '0.17.1'

[nelson-liu]

maybe related to issues # 6711 and 3864? Nevermind it's not.

[arogozhnikov]

#3864 is due to clipping, may be related to #6711, but seems no weights there.

[jnothman]

I can confirm getting "UndefinedMetricWarning: No negative samples in y_true, false positive value should be meaningless" but I'm not able to replicate the nan. Is there a reason you did not report the rest of your platform details?

[jnothman]

Scratch that. I do get nan. Just wasn't reading output correctly.

[jnothman]

This is what it comes down to:

In [25]: np.cumsum(trivial_weight.astype('float32'))[-1]
Out[25]: 16777216.0

In [26]: np.cumsum(trivial_weight.astype('float64'))[-1]
Out[26]: 40000000.0

In [27]: np.sum(trivial_weight.astype('float32'))
Out[27]: 40000000.0

In [28]: np.sum(trivial_weight.astype('float32')).dtype
Out[28]: dtype('float32')

Cumsum is losing precision, even when a float32 sum does not. This might be a numpy bug. If so, it lies in np.add.accumulate:

In [10]: np.add.accumulate(trivial_weight.astype('float32'))[-1]
Out[10]: 16777216.0

We can solve here by always passing dtype='float64' to cumsum, or doing so if the number of samples is high or something. There might also be a way to do cumsum in batches in a way that reduces the issue. Certainly we can identify the error by doing the cumsum and checking if its final value is equal to the sum.

[jnothman]

@arogozhnikov, please let me know if you report a bug to numpy. Feel free to submit a PR here.

and @nelson-liu, this is entirely independent of those other bugs.

[nelson-liu]

@jnothman ah ok, sorry I didn't look super closely into it; I just remembered seeing those other two threads and xref'd them here. thanks for letting me know.

[arogozhnikov]

@jnothman
it's not a bug, it is an expected behavior, since this code

x = trivial_weight
for i in range(len(x) - 1):
    x[i + 1] += x[i]

gives the same (wrong) result.

np.sum is a different story — one can change order of summation there in any possible way.
As for np.cumsum, this will drive to doubling number of computations (even for optimal implementation), so I don't think numpy team will do this (at least till the moment they need multithreading for np.cumsum).

[jnothman]

I don't find your claim to be true:

In [9]: x = trivial_weight.copy()

In [10]: for i in range(len(x) - 1):
   ....:         x[i + 1] += x[i]
   ....:
x[-1]

In [11]: x[-1]
Out[11]: 40000000.0

[jnothman]

whoops forgot astype. stupid.

[jnothman]

you're right.... though that's really quite nasty.

[jnothman]

Would we be happy with a fix that performs cumsum in float64?

[amueller]

+1 for cumsum in float64.

[arogozhnikov]

Would we be happy with a fix that performs cumsum in float64?

Yes, sounds good enough.

[jnothman]

Is there a way to test this? Would we better off with a (somewhat costly) helper to insure stability even in the float64 case? Like:

def stable_cumsum(arr):
    out = np.cumsum(arr, dtype=np.float64)
    expected = np.sum(arr, dtype=np.float64)
    if not np.allclose(out, expected):
        raise RuntimeError('cumsum was found to be unstable: its results do not correspond to sum')
    return out