np.linalg.norm is broken with float32 and arrays > 2 ** 28

[kohr-h]

I get this:

>>> a = np.ones(2 ** 28, dtype='float32')
>>> np.linalg.norm(a) == np.sqrt(a.size)
True
>>> a = np.ones(2 ** 28 + 2 ** 12 + 1024, dtype='float32')
>>> np.isclose(np.linalg.norm(a), np.sqrt(a.size))
True
>>> a = np.ones(2 ** 28 + 2 ** 12 + 2048, dtype='float32')
>>> np.isclose(np.linalg.norm(a), np.sqrt(a.size))
False

For float64, everything is fine even with size 2 ** 32, so 32-bit integer indices are not the immediate issue.

[eric-wieser]

Do you get the problem comparing a.sum() to a.size?

[kohr-h]

No, no problem there.

A potentially interesting finding: np.dot behaves like this:

>>> a = np.ones(2 ** 29, dtype='float32')
>>> int(np.dot(a, a))
268435456
>>> int(np.dot(a, a)) * 2 == a.size
True

The result is wrong by a factor of 2 exactly.

[kohr-h]

I get the same result with sdot from SciPy:

>>> sdot = scipy.linalg.blas.get_blas_funcs('dot', dtype='f')
>>> int(sdot(a, a, a.size))
268435456

I guess Numpy also uses BLAS dot in this case.

[eric-wieser]

On my machine it's wrong in a different way:

>>> int(np.dot(a, a))
67108864

[kohr-h]

On my machine it's wrong in a different way:

int(np.dot(a, a))
67108864

:-)

The issue only occurs with openblas it seems, just tested with MKL and results are correct.

[kohr-h]

Since this is likely platform-dependent I should add that I'm using 64 bit Linux.

[matthew-brett]

If it is a bug in OpenBLAS, we should raise an issue with them ...

[matthew-brett]

Here's me on numpy with Accelerate on Mac:

In [14]: a = np.ones(2 ** 29, dtype='float32')

In [15]: int(np.dot(a, a))
Out[15]: 268435456

[kohr-h]

If it is a bug in OpenBLAS, we should raise an issue with them ...

I'll do that.

[kohr-h]

I get a different but also wrong behavior of vanilla OpenBLAS with the develop branch. I'll test the release version that I currently use with Numpy next.

[kohr-h]

Repeating a comment on the OpenBLAS issue: How does Numpy use MKL for dot products? I'm asking since the MKL and OpenBLAS backends behave differently: the former does not suffer from this issue, the latter does.

[kohr-h]

From the upstream issue: there is also a dsdot function that takes float arrays and accumulates and returns in double precision, I think part of LAPACK. It's somewhat slower than sdot but gives correct results in the upstream issue case.

[ahmedahmedov]

I have a similar issue with np.dot.
np.dot(a,a) returns a negative number, which should never happen.
Related question with a working example here: (https://stackoverflow.com/questions/47206028/numpy-dot-product-of-two-same-vectors-returns-a-negative-number?noredirect=1#comment81364053_47206028)