Python: numpy.corrcoef Memory Error

Question

I was trying to calculate the correlation between a large set of data read from a text. For extremely large data set the program give a memory error. Can anyone please tell me how to correct this problem. Thanks

The following is my code:

enter code here

import numpy
from numpy  import *
from array import *
from decimal import *
import sys

Threshold = 0.8;
TopMostData = 10;

FileName = sys.argv[1]

File = open(FileName,'r')

SignalData = numpy.empty((1, 128));
SignalData[:][:] = 0;

for line in File:

    TempLine = line.split();
    TempInt = [float(i) for i in TempLine]
    SignalData = vstack((SignalData,TempInt))

del TempLine;
del TempInt;

File.close();

TempData = SignalData;
SignalData = SignalData[1:,:]
SignalData = SignalData[:,65:128]

print "File Read | Data Stored" + " | Total Lines: " + str(len(SignalData))

CorrelationData = numpy.corrcoef(SignalData)

The following is the error:

Traceback (most recent call last):
  File "Corelation.py", line 36, in <module>
    CorrelationData = numpy.corrcoef(SignalData)
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1824, in corrcoef
    return c/sqrt(multiply.outer(d, d))
MemoryError

Answer 1

You run out of memory as the comments show. If that happens because you are using 32-bit Python, even the method below will fail. But for the 64-bit Python and not-so-much-RAM situation we can do a lot as calculating the correlations is easily done piecewise, as you only need two lines in the memory simultaneously.

So, you may split your input into, say, 1000 row chunks, and then the resulting 1000 x 1000 matrices are easy to keep in memory. Then you can assemble your result into the big output matrix which is not necessarily in the RAM. I recommend this approach even if you have a lot of RAM, because this is much more memory-friendly. Correlation coefficient calculation is not an operation where fast random accesses would help a lot if the input can be kept in RAM.

Unfortunately, the numpy.corrcoef does not do this automatically, and we'll have to roll our own correlation coefficient calculation. Fortunately, that is not as hard as it sounds.

Something along these lines:

import numpy as np

# number of rows in one chunk
SPLITROWS = 1000

# the big table, which is usually bigger
bigdata = numpy.random.random((27000, 128))

numrows = bigdata.shape[0]

# subtract means form the input data
bigdata -= np.mean(bigdata, axis=1)[:,None]

# normalize the data
bigdata /= np.sqrt(np.sum(bigdata*bigdata, axis=1))[:,None]

# reserve the resulting table onto HDD
res = np.memmap("/tmp/mydata.dat", 'float64', mode='w+', shape=(numrows, numrows))

for r in range(0, numrows, SPLITROWS):
    for c in range(0, numrows, SPLITROWS):
        r1 = r + SPLITROWS
        c1 = c + SPLITROWS
        chunk1 = bigdata[r:r1]
        chunk2 = bigdata[c:c1]
        res[r:r1, c:c1] = np.dot(chunk1, chunk2.T)

Some notes:

• the code above is tested above np.corrcoef(bigdata)
• if you have complex values, you'll need to create a complex output array res and take the complex conjugate of chunk2.T
• the code garbles bigdata to maintain performance and minimize memory use; if you need to preserve it, make a copy

The above code takes about 85 s to run on my machine, but the data will mostly fit in RAM, and I have a SSD disk. The algorithm is coded in such order to avoid too random access into the HDD, i.e. the access is reasonably sequential. In comparison, the non-memmapped standard version is not significantly faster even if you have a lot of memory. (Actually, it took a lot more time in my case, but I suspect I ran out of my 16 GiB and then there was a lot of swapping going on.)

You can make the actual calculations faster by omitting half of the matrix, because res.T == res. In practice, you can omit all blocks where c > r and then mirror them later on. On the other hand, the performance is most likely limited by the HDD preformance, so other optimizations do not necessarily bring much more speed.

Of course, this approach is easy to make parallel, as the chunk calculations are completely independent. Also the memmapped array can be shared between threads rather easily.