Numpy：通过组合先前数组的子集来创建新数组

Question

我有一个二进制稀疏 CSR 数组。我想通过组合来自这个原始数组的列来创建一个新数组。也就是说，我有一个 "column groups" 的列表：[[1,10,3], [5,54,202], [12,199], [5], ...]

对于这些 "column groups" 中的每一个，我想将原始数组中的列与 OR 运算相结合（np.max 适用于此）并将组合的列添加到新矩阵中。

我目前的解决方案是使用 hstack，但它很慢：

for cg in column_groups:
    tmp = np.max(data_orig[:,cg].toarray(), axis=1, keepdims=True)
    data = np.hstack((data, tmp))

Answer 1

好吧，您基本上是在每次迭代时 selecting max 列。所以，我们可以 select 所有列，然后使用 np.maximum.reduceat 得到 "intervaled-maximum" 列，从而给我们一个向量化的解决方案，就像这样 -

def grouped_max(data_orig, column_groups):
    cols = np.hstack((column_groups))
    clens = np.hstack((0,np.cumsum(map(len,column_groups))[:-1]))
    all_data = data_orig[:,cols].toarray()
    return np.maximum.reduceat(all_data, clens,axis=1)

对于python 3.x版本，我们需要计算clens，像这样-

clens = np.hstack((0,np.cumsum(list(map(len,column_groups)))[:-1]))

由于循环版本沿组迭代，因此这种矢量化解决方案将在处理大量组时显示其优势。

样本运行-

In [303]: # Setup sample csr matrix
     ...: a = np.random.randint(0,3,(12,28))
     ...: data_orig = sparse.csr_matrix(a)
     ...: 
     ...: # Random column IDs
     ...: column_groups = [[1,10,3], [5,14],[2]]
     ...: 
     ...: data = np.empty((12,0),dtype=int)
     ...: for cg in column_groups:
     ...:     tmp = np.max(data_orig[:,cg].toarray(), axis=1, keepdims=True)
     ...:     data = np.hstack((data, tmp))
     ...:     

In [304]: out = grouped_max(data_orig, column_groups)

In [305]: # Verify results between original and propsed ones
     ...: print np.allclose(out, data)
True

Answer 2

我猜，主要问题是 hstack 在每次迭代中都会构建一个新矩阵 - 复制大量数据。

我没有使用过稀疏矩阵，所以我可能不合时宜，但据我从文档中了解到，可以将它们切片为普通的 numpy 数组。在这种情况下，最好预先分配数组并逐列添加结果：

rows = data_orig.shape[0]
cols = len(column_groups)
data = scipy.sparse.csr_matrix((rows, cols))
for cg in enumerate(column_groups):
    tmp = np.max(data_orig[:,cg[1]].toarray(), axis=1, keepdims=True)
    data[:, cg[0]] = tmp

Answer 3

In [412]: data_orig=sparse.random(10,300,.2,'csr')
In [413]: data_orig = (data_orig>.5).astype(int)
In [414]: cg = [[1,10,3], [5,54,202], [12,199], [5]]
In [420]: data_orig
Out[420]: 
<10x300 sparse matrix of type '<class 'numpy.int32'>'
   with 299 stored elements in Compressed Sparse Row format>

对于重复加入，最好附加到列表中，并堆叠一次：

def test1(data_orig, cg):
    data = []
    for g in cg:
        temp=np.max(data_orig[:,g].A, axis=1,keepdims=True)
        data.append(temp)
    return np.hstack(data)

In [429]: test1(data_orig, cg)
Out[429]: 
array([[0, 1, 0, 1],
       [1, 1, 0, 0],
       [0, 0, 0, 0],
       [0, 1, 0, 0],
       [0, 1, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 1, 0, 0],
       [1, 0, 0, 0],
       [0, 0, 0, 0]], dtype=int32)

我可以将稀疏矩阵转换为密集矩阵，并得到相同的结果

In [431]: dataM=data_orig.todense()
In [432]: test1(dataM, cg)

我本可以使用 data_orig.A，但我不得不在函数中省略 .A。

In [433]: timeit test1(data_orig, cg)
100 loops, best of 3: 2.52 ms per loop
In [434]: timeit test1(dataM, cg)
10000 loops, best of 3: 118 µs per loop

这些时间证实了我的猜测，即稀疏列索引相对较慢。

使用 Divakar 的 reduceat 版本：

In [451]: timeit grouped_max(data_orig, cg)
1000 loops, best of 3: 706 µs per loop
In [452]: timeit grouped_max(dataM, cg)
10000 loops, best of 3: 90.3 µs per loop

使用我自己的 Py3 改编（差异只是风格）：

def grouped_max(data_orig, column_groups):
    cols = np.hstack((column_groups))
    clens = np.hstack((0,np.cumsum([len(i) for i in cg])[:-1]))
    all_data = data_orig[:,cols].A
    return np.maximum.reduceat(all_data, clens,axis=1)

稀疏矩阵的加速非常显着。显然，一个较大的列 selection 比许多较小的列更快。密集矩阵的加速不是那么显着。

像这样的稀疏矩阵列select其实是用矩阵乘法来进行的。它为所需的列构造另一个 1s 的稀疏矩阵，并执行 dot 乘积。 Row/column 求和也是用矩阵乘积进行的。

这是一个纯稀疏版本：

def test2(data_orig, cg):
    data = []
    for g in cg:
        temp=data_orig[:,g].max(axis=1)
        data.append(temp)          # alt append(temp.A)
    return sparse.hstack(data)     #      np.hstack

In [465]: timeit test2(data_orig, cg).A
100 loops, best of 3: 3.21 ms per loop  (2.99 with np.hstack)

Numpy：通过组合先前数组的子集来创建新数组

Numpy: create new array by combining subsets of previous array

python

numpy

scipy

sparse-matrix