在Numpy中创建笛卡尔积时的内存错误

Question

我有3个numpy数组，需要在它们之间形成笛卡尔积。数组的维数不是固定的，因此它们可以采用不同的值，例如A=（10000，50），B=（40，50），C=（10000,50）。在

然后，我执行一些处理（比如a+b-c），下面是我用于产品的函数。在

def cartesian_2d(arrays, out=None):

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.shape[0] for x in arrays])
    if out is None:
        out = np.empty([n, len(arrays), arrays[0].shape[1]], dtype=dtype)

    m = n // arrays[0].shape[0]
    out[:, 0] = np.repeat(arrays[0], m, axis=0)
    if arrays[1:]:
        cartesian_2d(arrays[1:], out=out[0:m, 1:, :])
        for j in range(1, arrays[0].shape[0]):
            out[j * m:(j + 1) * m, 1:] = out[0:m, 1:]
    return out

a = [[ 0, -0.02], [1, -0.15]]
b = [[0, 0.03]]

result = cartesian_2d([a,b,a])

// array([[[ 0.  , -0.02],
    [ 0.  ,  0.03],
    [ 0.  , -0.02]],

   [[ 0.  , -0.02],
    [ 0.  ,  0.03],
    [ 1.  , -0.15]],

   [[ 1.  , -0.15],
    [ 0.  ,  0.03],
    [ 0.  , -0.02]],

   [[ 1.  , -0.15],
    [ 0.  ,  0.03],  
    [ 1.  , -0.15]]])

输出与itertools.product相同。但是，我使用我的自定义函数来利用numpy矢量化操作，与itertools.product对我来说。在

在这之后，我会的

^{pr2}$

所以这是最终的预期结果。在

只要数组适合内存，函数就可以正常工作。但是我的用例要求我处理大量的数据，我在np.empty()行得到一个MemoryError，因为它无法分配所需的内存。 我目前正在处理大约20GB的数据，将来可能会增加。在

这些数组表示向量，必须存储在float中，因此我不能使用int。而且，它们是密集数组，所以使用sparse不是一个选项。在

我将使用这些数组进行进一步的处理，理想情况下我不想在这个阶段将它们存储在文件中。因此memmap/h5py格式可能没有帮助，尽管我不确定这一点。在

如果有其他方法来形成这种产品，那也没问题。在

因为我确信有些应用程序的数据集比这个大得多，我希望有人以前遇到过这样的问题，并想知道如何处理这个问题。请帮忙。在

Answer 1

在磁盘上高效地写入结果

起初，一些人会考虑结果数据的大小。在

结果数据的大小

size_in_GB = A.shape[0]**2*A.shape[1]*B.shape[0]*(size_of_datatype)/1e9

在你的问题中，你提到了A.shape=（10000,50），B=（40,50）。使用float64，您的结果大约是1600gb。如果您有足够的磁盘空间，这可以毫无问题地完成，但是您必须考虑下一步不能对数据做什么。也许这只是一个中间的结果，处理块中的数据是可能的。在

如果不是这样的话，这里有一个如何有效地处理1600GB数据的例子（RAM使用量大约为200MB）。根据实际数据，troughput应该在200MB/s左右。在

计算结果的代码来自@PaulPanzer。在

^{pr2}$

读取数据

File_Name_HDF5='Test.h5'
f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300)
dset=f['Data']
chunks_size=500
for i in range(0,dset.shape[0],chunks_size):
  #Iterate over the first column
  data=dset[i:i+chunks_size,:] #avoid excessive calls to the hdf5 library
  #Do something with the data

f.close()

f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300)
dset=f['Data']
for i in range(dset.shape[1]):
  # Iterate over the second dimension
  # fancy indexing e.g.[:,i] will be much slower
  # use np.expand_dims or in this case np.squeeze after the read operation from the dset
  # if you wan't to have the same result than [:,i] (1 dim array)
  data=dset[:,i:i+1] 
  #Do something with the data

f.close()

在这个测试示例中，我获得了大约550m/s的写入吞吐量，大约（500 M/s的第一个dim，1000M/s的第二个dim）的读吞吐量和50的压缩比。numpymemmap只会提供可接受的速度，如果你沿着变化最快的方向（在最后一个维度的C中）读写数据，这里使用HDF5使用的分块数据格式，这根本不是问题。使用numpymemmap也不可能进行压缩，这会导致更大的文件大小和较慢的速度。在

请注意，压缩过滤器和块形状必须根据您的需要进行设置。这取决于如何避免事后读取数据和实际数据。在

如果你做了完全错误的事情，与正确的方法相比，执行速度可能会慢10-100倍（例如，可以为第一个或第二个阅读示例优化chunkshape）。在