在 Numpy 中创建笛卡尔积时出现 MemoryError
MemoryError while creating cartesian product in Numpy
我有 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 相比，它工作得很好。
之后，我做
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 格式可能没有帮助，虽然我不确定这一点。
如果有其他方法可以形成这个产品，那也可以。
因为我确信有比这更大的数据集的应用程序，我希望有人以前遇到过这样的问题并且想知道如何处理这个问题。请帮忙。
另一种解决方案是 create a cartesian product of indices（这更容易，因为存在一维数组的笛卡尔积的解决方案）：
idx = cartesian_product(
np.arange(len(a)),
np.arange(len(b)) + len(a),
np.arange(len(a))
)
然后使用花式索引创建输出数组：
x = np.concatenate((a, b))
result = x[idx.ravel(), :].reshape(*idx.shape, -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。如果一个操作数的维度少于其他维度，则其形状将填充为左。请注意，插图中显示的 equiv 数组并未显式创建。
如果结果也不合适
在那种情况下，我看不出如何避免使用存储，所以 h5py 或类似的东西。
从每个操作数中删除第一列
这只是切片的问题：
a_no_id = a[:, 1:]
等请注意，与 Python 列表不同，NumPy 数组在切片时不会 return 一个副本而是一个视图。因此效率（内存或运行时）在这里不是问题。
在磁盘上高效写入结果
首先对结果数据的大小有一些想法。
结果数据的大小
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 你的结果将是 aproximately 1600 GB。如果您有足够的磁盘 space，这可以毫无问题地完成，但是您必须考虑接下来您不想对数据做什么。也许这只是一个中间结果，可以按块处理数据。
如果不是这种情况，这里有一个如何有效处理 1600GB 数据的示例（RAM 使用量约为 200 MB）。真实数据的吞吐量应该在 200 左右 MB/s。
计算结果的代码来自@PaulPanzer。
import numpy as np
import tables #register blosc
import h5py as h5
import h5py_cache as h5c
a=np.arange(500*50).reshape(500, 50)
b=np.arange(40*50).reshape(40, 50)
# isn't well documented, have a look at https://github.com/Blosc/hdf5-blosc
compression_opts=(0, 0, 0, 0, 5, 1, 1)
compression_opts[4]=9 #compression level 0...9
compression_opts[5]=1 #shuffle
compression_opts[6]=1 #compressor (I guess that's lz4)
File_Name_HDF5='Test.h5'
f = h5.File(File_Name_HDF5, 'w',chunk_cache_mem_size=1024**2*300)
dset = f.create_dataset('Data', shape=(a.shape[0]**2*b.shape[0],a.shape[1]),dtype='d',chunks=(a.shape[0]*b.shape[0],1),compression=32001,compression_opts=(0, 0, 0, 0, 9, 1, 1), shuffle=False)
#Write the data
for i in range(a.shape[0]):
sh = a.shape[1:]
aba = (a[i] + b[:, None] - a).reshape(-1, *sh)
dset[i*a.shape[0]*b.shape[0]:(i+1)*a.shape[0]*b.shape[0]]=aba
f.close()
读取数据
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()
在这个测试示例中，我获得了大约 550 M/s 的写入吞吐量，大约（500 M/s 第一个暗淡，1000 M/s 第二个暗淡）的读取吞吐量和压缩比率为 50。如果您沿着最快变化的方向（在 C 中最后一个维度）读取或写入数据，Numpy memmap 将仅提供可接受的速度，这里使用 HDF5 使用的分块数据格式，这根本不是问题。 Numpy memmap 也无法进行压缩，导致文件更大，速度更慢。
请注意，必须根据您的需要设置压缩过滤器和块形状。这取决于您之后不想读取数据的方式和实际数据。
如果你做的事情完全错误，与正确的方法相比，性能可能会慢 10-100 倍（例如，可以针对第一个或第二个读取示例优化块形状）。
我有 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 相比，它工作得很好。
之后，我做
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 格式可能没有帮助，虽然我不确定这一点。
如果有其他方法可以形成这个产品，那也可以。
因为我确信有比这更大的数据集的应用程序，我希望有人以前遇到过这样的问题并且想知道如何处理这个问题。请帮忙。
另一种解决方案是 create a cartesian product of indices（这更容易，因为存在一维数组的笛卡尔积的解决方案）：
idx = cartesian_product(
np.arange(len(a)),
np.arange(len(b)) + len(a),
np.arange(len(a))
)
然后使用花式索引创建输出数组：
x = np.concatenate((a, b))
result = x[idx.ravel(), :].reshape(*idx.shape, -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。如果一个操作数的维度少于其他维度，则其形状将填充为左。请注意，插图中显示的 equiv 数组并未显式创建。
如果结果也不合适
在那种情况下，我看不出如何避免使用存储，所以 h5py 或类似的东西。
从每个操作数中删除第一列
这只是切片的问题：
a_no_id = a[:, 1:]
等请注意，与 Python 列表不同，NumPy 数组在切片时不会 return 一个副本而是一个视图。因此效率（内存或运行时）在这里不是问题。
在磁盘上高效写入结果
首先对结果数据的大小有一些想法。
结果数据的大小
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 你的结果将是 aproximately 1600 GB。如果您有足够的磁盘 space，这可以毫无问题地完成，但是您必须考虑接下来您不想对数据做什么。也许这只是一个中间结果，可以按块处理数据。
如果不是这种情况，这里有一个如何有效处理 1600GB 数据的示例（RAM 使用量约为 200 MB）。真实数据的吞吐量应该在 200 左右 MB/s。
计算结果的代码来自@PaulPanzer。
import numpy as np
import tables #register blosc
import h5py as h5
import h5py_cache as h5c
a=np.arange(500*50).reshape(500, 50)
b=np.arange(40*50).reshape(40, 50)
# isn't well documented, have a look at https://github.com/Blosc/hdf5-blosc
compression_opts=(0, 0, 0, 0, 5, 1, 1)
compression_opts[4]=9 #compression level 0...9
compression_opts[5]=1 #shuffle
compression_opts[6]=1 #compressor (I guess that's lz4)
File_Name_HDF5='Test.h5'
f = h5.File(File_Name_HDF5, 'w',chunk_cache_mem_size=1024**2*300)
dset = f.create_dataset('Data', shape=(a.shape[0]**2*b.shape[0],a.shape[1]),dtype='d',chunks=(a.shape[0]*b.shape[0],1),compression=32001,compression_opts=(0, 0, 0, 0, 9, 1, 1), shuffle=False)
#Write the data
for i in range(a.shape[0]):
sh = a.shape[1:]
aba = (a[i] + b[:, None] - a).reshape(-1, *sh)
dset[i*a.shape[0]*b.shape[0]:(i+1)*a.shape[0]*b.shape[0]]=aba
f.close()
读取数据
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()
在这个测试示例中，我获得了大约 550 M/s 的写入吞吐量，大约（500 M/s 第一个暗淡，1000 M/s 第二个暗淡）的读取吞吐量和压缩比率为 50。如果您沿着最快变化的方向（在 C 中最后一个维度）读取或写入数据，Numpy memmap 将仅提供可接受的速度，这里使用 HDF5 使用的分块数据格式，这根本不是问题。 Numpy memmap 也无法进行压缩，导致文件更大，速度更慢。
请注意，必须根据您的需要设置压缩过滤器和块形状。这取决于您之后不想读取数据的方式和实际数据。
如果你做的事情完全错误，与正确的方法相比，性能可能会慢 10-100 倍（例如，可以针对第一个或第二个读取示例优化块形状）。