[Python] Weird behavior in extremely-large complex128 NumPy arrays, imaginary-part only

I'm a physicist working on numerical simulations. I recently ran into a worrying issue with extremely large NumPy arrays (~ 26 GB). The system in question is a linux machine with 128GB of RAM. The arrays are of type complex128 (this is important).

The issue is as follows:

• Arrays are instantiated without errors (or if they fail, they do so quietly)
• Assigning array values proceeds similarly without complaint
• Memory appears to be properly allocated, or at least the correct amount is reserved by the process

At this point, I query the minimum and maximum values of the real and imaginary parts:

• The real parts are always correct
• The imaginary parts are often extremely incorrect

Here's a working example that demonstrates the issue, at least on my specific system:

import numpy as np
import numba


@numba.jit(cache=True)
def minmaxrealimag(x):
rmaximum = x[0].real
rminimum = x[0].real
imaximum = x[0].imag
iminimum = x[0].imag

for i in x[1:]:
if i.real > rmaximum:
rmaximum = i.real
elif i.real < rminimum:
rminimum = i.real
if i.imag > imaximum:
imaximum = i.imag
elif i.imag < iminimum:
iminimum = i.imag
return (rminimum, rmaximum, iminimum, imaximum)


testn = 2048
testx = 1024


field = np.empty((testx, testx), dtype=complex)
field.real[:] = np.random.rand(*field.shape)[:]
field.imag[:] = np.random.rand(*field.shape)[:]
profile = np.empty(testn, dtype=complex)
profile.real[:] = np.random.rand(*profile.shape)[:]
profile.imag[:] = np.random.rand(*profile.shape)[:]
correctslice = field[:, :] * profile[0]
(rmin, rmax, imin, imax) = minmaxrealimag(correctslice.flatten())

test1 = np.empty((testn, testx, testx), dtype=complex)
for itau in range(testn):
correctslice = field[:, :] * profile[itau]
(rmin2, rmax2, imin2, imax2) = minmaxrealimag(correctslice.flatten())
if rmin2 < rmin:
rmin = rmin2
elif rmax2 > rmax:
rmax = rmax2
if imin2 < imin:
imin = imin2
elif imax2 > imax:
imax = imax2

test1[itau] = correctslice[:, :]

print((rmin, rmax, imin, imax))
print(minmaxrealimag(test1.flatten()))



This seems over-complicated, but a much-simpler example (though somewhat slower and less informative) is:

import numpy as np

testn = 2048
testx = 1024


field = np.empty((testn, testx, testx), dtype=complex)
field.real[:] = np.random.rand(*field.shape)[:]
field.imag[:] = np.random.rand(*field.shape)[:]

print(np.max(field.imag))


Sometimes everything goes fine, but usually, the minimum and maximum of the imaginary parts are wildly incorrect. For these examples, the minimum is reasonable but wrong, and the maximum is either nan (very rare) or (much more likely) close to the double-precision floating point maximum, though a factor of two or four smaller. Sometimes, the maximum has an exponent that a factor of two or four smaller than than the float64 max. I've never seen it be Inf. From this, I assume that this weird maximum value in the imaginary part has the following properties:

• The sign bit is always 0
• The mantissa is usually all 1s
• Bit 11 is almost always 0
• Bit 12 is almost always 1

When I try to locate and explicitly assign a value (say, 1) to the offending pixel, the value remains unchanged.

This is really driving me absolutely insane, so I hope someone can tell me that I've just been stupid and the solution is obvious. The issue isn't always replicable in these little snippets, but it's consistently stalling my simulations that use multiple arrays of this size and frequently occupy >80% of system memory. The most frustrating part is the lack of an error, segfault, or exception, or else I'd think this was a memory-safety issue or something along those lines. Not even a warning, though! I have no way to know if a simulation is going to work, or when it's going off the rails.

The obvious answer is to move to a distributed-memory model of parallelization, and I will, but I first want to know that the issue isn't on my end, aside from asking too much of my poor computer.

PLEASE HELP!

PS I haven't tried this on a different computer system, as none of my other workhorses have as much memory. However, the computer in question is behaving normally in every other way.

Continue reading...