在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

另一种解决方案是create a cartesian product of indices（这更简单，因为存在一维数组的笛卡尔积的解决方案）：

idx = cartesian_product(
    np.arange(len(a)),
    np.arange(len(b)) + len(a),
    np.arange(len(a))
)

然后使用花式索引创建输出数组：

^{pr2}$

Answer 2

在磁盘上高效地写入结果

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

结果数据的大小

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）。在

Answer 3

如果你的结果符合记忆的话

在不依赖三倍于结果大小的中间值的情况下，下面的结果将生成您的预期结果。它使用广播。在

请注意，几乎任何NumPy操作都可以这样广播，因此在实践中可能不需要显式的笛卡尔积：

#shared dimensions:
sh = a.shape[1:]
aba = (a[:, None, None] + b[None, :, None] - a[None, None, :]).reshape(-1, *sh)
aba
#array([[ 0.  ,  0.03],
#       [-1.  ,  0.16],
#       [ 1.  , -0.1 ],
#       [ 0.  ,  0.03]])

按“ID”寻址结果行

您可以考虑省略reshape。这将允许您通过组合索引对结果中的行进行寻址。如果你的组件ID只有0,1,2，。。。与您的示例一样，这将与组合ID相同。例如aba[1,0,0]将对应于a的第二行+b的第一行-a的第一行

一点解释

广播：例如，当添加两个数组时，它们的形状不必完全相同，只是由于广播的关系而兼容。广播在某种意义上是向数组添加标量的一种概括：

^{pr2}$

广播：

              [[4],            [[1, 2, 3],   [[4, 4, 4],
[[1, 2, 3]] +  [5],  equiv to   [1, 2, 3], +  [5, 5, 5],
               [6]]             [1, 2, 3]]    [6, 6, 6]]

为此，每个操作数的每个维度必须为1或等于其他操作数中的相应维度，除非它是1。如果一个操作数的维数比其他操作数少，它的形状将用左边的维数填充。请注意，图中显示的equiv数组不是显式创建的。在

如果结果也不符合

在这种情况下，我不知道如何才能避免使用存储，所以h5py或类似的东西。在

从每个操作数中删除第一列

这只是一个切片问题：

a_no_id = a[:, 1:]

注意，与Python列表不同，NumPy数组在切片时不返回副本，而是返回一个视图。因此，效率（内存或运行时）在这里不是问题。在