2013-10-14 222 views
6

要在磁盤上存儲大矩陣,我使用numpy.memmap。Numpy高效矩陣乘法

這裏是一個示例代碼來測試大矩陣乘法:

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) 

所以我的問題是:

  1. 如何限制aplication的內存單耗這是使用此過程,例如一些價值到100Mb(或1Gb或別的東西)。我也不知道如何估計過程的內存消耗(我認爲內存僅在創建「數據」變量時分配,但在使用memmap文件時使用了多少內存?)
  2. 也許存在一些存儲在磁盤上的大矩陣乘法的最佳解決方案?例如,可能數據沒有最優地存儲在磁盤上或從磁盤上獲取,沒有正確分配,並且dot產品只使用一個內核。也許我應該使用PyTables之類的東西?

另外我對解決線性方程組(SVD和其他)有限的內存使用情況的算法感興趣。 也許這種算法叫核心或迭代,我覺得有一些比喻像硬盤驅動器< - > ram,gpu ram < - > cpu ram,cpu ram < - > cpu cache。

另外here我在PyTables中找到了關於矩陣乘法的一些信息。

此外,我發現this在R,但我需要它的Python或Matlab。

+0

「如何限制正在使用此過程的應用程序的內存消耗達到某個值,例如100Mb」您的意思是,如果應用程序嘗試使用更多的內存,它應該會失敗?使用'psutil.set_rlimit'很容易,但AFAIK只能在linux上運行。 – Bakuriu

+0

不,我的意思是應用程序必須正常工作,但使用少於聲明的內存(一般來說,它會更慢,更少的內存,但它有用,當我們想限制應用程序內存使用或如果我們沒有足夠的內存適合整個矩陣)。我在Windows上工作。 – mrgloom

+0

你的'res'行沒有意義(而res是最大的數組......)。重讀'np.dot'文檔字符串,你會發現一些有用的... – seberg

回答

1

考慮使用NumExpr爲您處理:https://github.com/pydata/numexpr

...內部,NumExpr採用了圍繞分塊讀取策略設計了自己的量化的虛擬機,以最佳尺寸塊有效地運作上內存中的數據。如果調整得當,它可以輕鬆擊敗幼稚的NumPy操作。

NumExpr可能會覆蓋問題明細中的#2。如果您使用流式傳輸的二進制格式地址#1,你就可以分塊讀取加載數據文件時的方法 - 就像這樣:

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

如果太低級了你,我會建議你看看HDF5庫和格式:http://www.h5py.org - 這是我所知道的基於NumPy的二進制序列化結構的最佳解決方案。 h5py模塊支持壓縮,分塊閱讀,dtypes,元數據......你的名字。

祝你好運!

+0

我不太明白在矩陣乘法中可以使用numexpr,可以擴展答案並提供使用示例嗎?是的,我知道hdf5和pytables,但numpy.memmap更方便,因爲它可以像平常一樣使用numpy數組。 – mrgloom

+0

'numpy.memmap'是一個全或無的操作,是的?...我建議hdf5,因爲它支持分塊和流I/O:http://www.hdfgroup.org/HDF5/doc/Advanced/Chunking/ index.html - 雖然簡單地使用HDF5加載相同的數據可能比內存映射文件慢,但對基於塊的處理循環進行分塊讀取可能會更有效。 – fish2000

+0

關於NumExpr,我不知道它是否提供點積算子,但我確實知道它會爲您提供更快的視圖操作 - 如果您要使用NumExpr重做您的數組操作,則可能不需要例如耗費內存的換位。 – fish2000

3

Dask.array爲使用阻塞算法和任務調度的大型磁盤陣列提供了一個numpy接口。它可以輕鬆地進行非核心矩陣乘法和其他簡單的numpy操作。

阻塞線性代數很難,你可能想看看關於這個話題的一些學術研究。 Dask確實支持高維矩陣上的QR和SVD因式分解。

不管大數組,你真的想要阻塞的算法,而不是天真的遍歷,它會以不愉快的方式擊中磁盤。