两个大向量之间的元素比较，高度稀疏答案

【问题标题】：Elementwise comparison between two large vectors, high degree of sparsity两个大向量之间的元素比较，高度稀疏
【发布时间】：2021-07-23 21:45:02
【问题描述】：

需要一个执行类似于numpy.where 函数的函数，但不会遇到由布尔数组的密集表示引起的内存问题。因此，该函数应该能够返回一个极其稀疏的布尔数组。

虽然下面给出的示例适用于小型数据集/向量，但不可能使用 numpy.where 函数，例如，my_sample 的形状为 (10.000.000, 1) 和 my_population 的形状为 @987654327 @。阅读其他线程后，numpy.where 在评估表达式numpy.where((my_sample == my_population.T)) 时显然创建了一个形状为(10.000.000, 100.000) 的密集布尔数组。这个密集的(10.000.000, 100.000) 数组无法放入我的机器/大多数机器的内存中。

生成的数组非常稀疏。就我而言，要知道每行最多有两个 1！使用上面的规范，稀疏度等于 0.002%。这绝对应该适合记忆。

尝试为数值模拟创建类似于模型/设计矩阵的东西。生成的矩阵将用于一些线性代数运算。

最小的工作示例：请注意向量中的位置/坐标很重要。

# import packages
import numpy as np

# my_sample is the vector of observations
my_sample = ['a', 'b', 'c', 'a']

# my_population is the lookup vector
my_population = ['a', 'b', 'c']

# initalise the matrix (dense matrix for this exampe)
my_zero = np.zeros((len(my_sample), len(my_population)))

# reshape to arrays
my_sample = np.array(my_sample).reshape(-1, 1)
my_population = np.array((my_population)).reshape(-1, 1)

# THIS STEP CAUSES THE MEMORY ISSUES
my_indices = np.where((my_sample == my_population.T))

# set the matches to equal one
my_zero[my_indices] = 1

# show matrix
my_zero
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.],
       [1., 0., 0.]])

【问题讨论】：

您在寻找“更好的性能”吗？如果是这样，您可以尝试使用numba。这样做有很多好处
@JohnBrookfields 我正在寻找一种适合普通机器内存的解决方案（比如 16GB 的 RAM）。但是感谢您指出这一点，我将编辑我的问题。
numba 也有内存管理。见here
需要明确的是，== 操作占用了大量内存。 where 只是找到该数组的 True 元素。
@hpaulj 非常感谢您的澄清！

标签： python numpy scipy sparse-matrix elementwise-operations

【解决方案1】：

首先，让我们将其编码为整数，而不是字符串。字符串很烂。

pop_levels, pop_idx = np.unique(my_population, return_inverse=True)
sample_levels, sample_idx = np.unique(my_sample, return_inverse=True)

pop_levels 和 sample_levels 必须相同，这一点很重要，但如果它们是相同的，那么您就完成了 - 将它们打包成稀疏掩码：

sample_mask = sps.csr_matrix((np.ones_like(sample_idx), sample_idx, range(len(sample_idx) + 1)))

我们完成了：

>>> sample_mask.A
array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1],
       [1, 0, 0]])

您可能需要重新排序您的因子水平，以便它们在您的样本和总体之间是相同的，但只要您可以统一这些标签，这非常简单，只需矩阵分配即可。

【讨论】：

非常感谢@CJR 的回答。我已将其标记为解决方案，因为它几乎完全符合我的需要。我想知道如果pop_levels 和sample_levels 不同，您是否知道该怎么做。假设pop_levels 中的另一个元素不会出现在my_sample/sample_levels 中。

【解决方案2】：

更直接的路线：

In [125]: my_sample = ['a', 'b', 'c', 'a']
     ...: my_population = ['a', 'b', 'c']
     ...: 
     ...: 
In [126]: np.array(my_sample)[:,None]==np.array(my_population)
Out[126]: 
array([[ True, False, False],
       [False,  True, False],
       [False, False,  True],
       [ True, False, False]])

这是一个布尔数据类型。如果你想要 0/1 整数矩阵：

In [128]: (np.array(my_sample)[:,None]==np.array(my_population)).astype(int)
Out[128]: 
array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1],
       [1, 0, 0]])

如果你有空间来创建my_zero，你应该有空间来创建这个数组。如果大的临时缓冲区仍然存在问题，您可以尝试转换为“uint8”，这样会占用更少的空间。

在您的版本中，您创建了两个大数组，my_zero 和 my_sample == my_population.T。但请注意，即使您通过了这一步，您也可能没有空间使用my_zero 做任何其他事情。

创建稀疏矩阵可能会节省空间，但稀疏度必须非常高才能保持任何速度。虽然矩阵乘法是 scipy.sparse 矩阵的一个相对强大的领域。

时间测试

In [134]: %%timeit
     ...: pop_levels, pop_idx = np.unique(my_population, return_inverse=True)
     ...: sample_levels, sample_idx = np.unique(my_sample, return_inverse=True)
     ...: sample_mask = sparse.csr_matrix((np.ones_like(sample_idx), sample_idx, range(len(s
     ...: ample_idx) + 1)))

247 µs ± 195 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [135]: timeit (np.array(my_sample)[:,None]==np.array(my_population)).astype(int)
9.61 µs ± 9.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

并从稠密中制作一个稀疏矩阵：

In [136]: timeit sparse.csr_matrix((np.array(my_sample)[:,None]==np.array(my_population)).as
     ...: type(int))
332 µs ± 1.44 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

大型阵列的缩放可能完全不同。

【讨论】：