在Numpy中创建笛卡尔积时的MemoryError

Question

我有3个numpy阵列,需要在它们之间形成笛卡尔积.阵列的尺寸不固定,因此它们可以采用不同的值,一个例子可以是A =(10000,50),B =(40,50),C =(10000,50).

然后,我执行一些处理(如a + bc)以下是我用于产品的功能.

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相比,它运行良好.

在此之后,我做到了

result[:, 0, :] + result[:, 1, :] - result[:, 2, :]

//array([[ 0.  ,  0.03],
       [-1.  ,  0.16],
       [ 1.  , -0.1 ],
       [ 0.  ,  0.03]])

所以这是最终的预期结果.

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

这些数组代表向量,必须存储float,所以我不能使用int.此外,它们是密集阵列,因此使用sparse不是一种选择.

我将使用这些数组进行进一步处理,理想情况下我不想在这个阶段将它们存储在文件中.所以memmap/ h5py格式可能没有帮助,虽然我不确定.

如果还有其他方法可以形成这个产品,那也没关系.

我确信有些应用程序的数据集大于此,我希望有人之前遇到过这样的问题,并希望知道如何处理这个问题.请帮忙.

Answer 1

如果至少你的结果适合记忆

以下代码会产生您的预期结果，而不依赖于结果大小三倍的中间值。它使用广播。

请注意，几乎所有 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 的第一行所获得的行。

一点解释

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

    [[2],                 [[7],   [[2],
7 +  [3],     equiv to     [7], +  [3],
     [4]]                  [7]]    [4]]

广播：

              [[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）。如果一个操作数的维度少于其他操作数，则其形状将在左侧填充1。请注意，图中显示的等效数组并未显式创建。

如果结果也不适合

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

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

这只是切片的问题：

a_no_id = a[:, 1:]

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