Numpy efficient big matrix multiplication

Question

To store big matrix on disk I use numpy.memmap.

Here is a sample code to test big matrix multiplication:

import numpy as np
import time

rows= 10000 # it can be large for example 1kk
cols= 1000

#create some data in memory
data = np.arange(rows*cols, dtype='float32') 
data.resize((rows,cols))

#create file on disk
fp0 = np.memmap('C:/data_0', dtype='float32', mode='w+', shape=(rows,cols))
fp1 = np.memmap('C:/data_1', dtype='float32', mode='w+', shape=(rows,cols))

fp0[:]=data[:]
fp1[:]=data[:]

#matrix transpose test
tr = np.memmap('C:/data_tr', dtype='float32', mode='w+', shape=(cols,rows))
tr= np.transpose(fp1)  #memory consumption?
print fp1.shape
print tr.shape

res = np.memmap('C:/data_res', dtype='float32', mode='w+', shape=(rows,rows))
t0 = time.time()
# redifinition ? res= np.dot(fp0,tr) #takes 342 seconds on my machine, if I multiplicate matrices in RAM it takes 345 seconds (I thinks it's a strange result)
res[:]= np.dot(fp0,tr) # assignment ?
print res.shape
print (time.time() - t0)

So my questions are :

1. How to restrict memory consumtion of aplication which is using this procedure to some value for example to 100Mb(or 1Gb or something else).Also I don't understand how to estimate memory consumtion of procedure (I think memory is only allocated when "data" variable is created, but how much memory used when we use memmap files?)
2. Maybe there is some optimal solution for multiplication of big matrices stored on disk? For example maybe data not optimally stored on disk or readed from disk, not properly chached, and also dot product use only one core.Maybe I should use something like PyTables?

Also I interested in algorithms solving linear system of equations (SVD and others) with restricted memory usage. Maybe this algorithms called out-of-core or iterative and I think there some analogy like hard drive<->ram, gpu ram<->cpu ram, cpu ram<->cpu cache.

Also here I found some info about matrix multiplication in PyTables.

Also I found this in R but I need it for Python or Matlab.

Answer 1

Dask.array provides a numpy interface to large on-disk arrays using blocked algorithms and task scheduling. It can easily do out-of-core matrix multiplies and other simple-ish numpy operations.

Blocked linear algebra is harder and you might want to check out some of the academic work on this topic. Dask does support QR and SVD factorizations on tall-and-skinny matrices.

Regardless for large arrays, you really want blocked algorithms, not naive traversals which will hit disk in unpleasant ways.

Answer 2

Consider using NumExpr for your processing: https://github.com/pydata/numexpr

... internally, NumExpr employs its own vectorized virtual machine that is designed around a chunked-read strategy, in order to efficiently operate on optimally-sized blocks of data in memory. It can handily beat naïve NumPy operations if tuned properly.

NumExpr may cover #2 in your breakdown of the issue. If you address #1 by using a streamable binary format, you can then the chunked-read approach when loading your data files – like so:

    with open('path/to/your-data.bin', 'rb') as binary:
        while True:
            chunk = binary.read(4096) # or what have you
            if not chunk:
                break

If that is too low-level for you, I would recommend you look at the HDF5 library and format: http://www.h5py.org – it’s the best solution for the binary serialization of NumPy-based structures that I know of. The h5py module supports compression, chunked reading, dtypes, metadata… you name it.

Good luck!