logo标签列表

使用NumPy进行向量化分组

pythonnumpy
17
Pandas拥有广泛使用的groupby功能，可以根据相应的映射将DataFrame分割成子组，您可以在每个子组上应用计算并重新合并结果。
在NumPy中是否可以灵活地完成此操作而不需要本机Python for循环？ 通过Python循环实现该操作将如下所示：
>>> import numpy as np

>>> X = np.arange(10).reshape(5, 2)
>>> groups = np.array([0, 0, 0, 1, 1])

# Split up elements (rows) of `X` based on their element wise group
>>> np.array([X[groups==i].sum() for i in np.unique(groups)])
array([15, 30])

15以上是X的前三行的总和，30是剩余两行的总和。

通过“灵活地”，我只是指我们不专注于一种特定的计算，例如求和、计数、最大值等，而是将任何计算传递给分组数组。

如果没有，请问是否有比上述方法更快的方法？

@juanpa.arrivillaga 对的...而且itertools的key参数需要是可调用的（不是一组组）。你不能使用像groups.__getitem__这样的东西，因为它需要是一个应用于第一个参数的每个元素的函数。 - Brad Solomon
你有没有想过pandas是如何实现groupby的？它使用Python结构（字典？itertools groupby？），C或Cython代码吗？itertools版本将连续的组分组在一起。因此，输入（键）必须排序。 - hpaulj
这有点像垃圾箱问题：[如何使用numpy数组有效地获取特定值选定的索引列表？]（https://stackoverflow.com/questions/48686381）。我认为很难得到真正“向量化”的解决方案，因为箱（或组）的大小不同。如果您只想要一个索引，@ Divikar提供了一个searchsorted解决方案，以向量化的方式检索多个值的索引。 - hpaulj
@hpaulj 这个 Python包装器 基本上构建了一个（键，数组）对的迭代器。它确实会为一些实际的计算做 Cython 调用。 - Brad Solomon
是的，始终为1d无符号uint8，具体而言@DanielF（但不保证已排序） - Brad Solomon
显示剩余2条评论
5个回答
9

使用scipy稀疏矩阵如何？

import numpy as np
from scipy import sparse
import time

x_len = 500000
g_len = 100

X = np.arange(x_len * 2).reshape(x_len, 2)
groups = np.random.randint(0, g_len, x_len)

# original
s = time.time()

a = np.array([X[groups==i].sum() for i in np.unique(groups)])

print(time.time() - s)

# using scipy sparse matrix
s = time.time()

x_sum = X.sum(axis=1)
b = np.array(sparse.coo_matrix(
    (
        x_sum,
        (groups, np.arange(len(x_sum)))
    ),
    shape=(g_len, x_len)
).sum(axis=1)).ravel()

print(time.time() - s)

#compare
print(np.abs((a-b)).sum())

我的电脑上的结果

0.15915322303771973
0.012875080108642578
0

比原来快了10倍以上。

更新！

让我们对 @Paul Panzer 和 @Daniel F. 的答案进行基准测试。这只是一个求和基准测试。

import numpy as np
from scipy import sparse
import time

# by @Daniel F
def groupby_np(X, groups, axis = 0, uf = np.add, out = None, minlength = 0, identity = None):
    if minlength < groups.max() + 1:
        minlength = groups.max() + 1
    if identity is None:
        identity = uf.identity
    i = list(range(X.ndim))
    del i[axis]
    i = tuple(i)
    n = out is None
    if n:
        if identity is None:  # fallback to loops over 0-index for identity
            assert np.all(np.in1d(np.arange(minlength), groups)), "No valid identity for unassinged groups"
            s = [slice(None)] * X.ndim
            for i_ in i:
                s[i_] = 0
            out = np.array([uf.reduce(X[tuple(s)][groups == i]) for i in range(minlength)])
        else:
            out = np.full((minlength,), identity, dtype = X.dtype)
    uf.at(out, groups, uf.reduce(X, i))
    if n:
        return out

x_len = 500000
g_len = 200

X = np.arange(x_len * 2).reshape(x_len, 2)
groups = np.random.randint(0, g_len, x_len)

print("original")
s = time.time()

a = np.array([X[groups==i].sum() for i in np.unique(groups)])

print(time.time() - s)

print("use scipy coo matrix")
s = time.time()

x_sum = X.sum(axis=1)
b = np.array(sparse.coo_matrix(
    (
        x_sum,
        (groups, np.arange(len(x_sum)))
    ),
    shape=(g_len, x_len)
).sum(axis=1)).ravel()

print(time.time() - s)

#compare
print(np.abs((a-b)).sum())


print("use scipy csr matrix @Daniel F")
s = time.time()
x_sum = X.sum(axis=1)
c = np.array(sparse.csr_matrix(
    (
        x_sum,
        groups,
        np.arange(len(groups)+1)
    ),
    shape=(len(groups), g_len)
).sum(axis=0)).ravel()

print(time.time() - s)

#compare
print(np.abs((a-c)).sum())


print("use bincount @Paul Panzer @Daniel F")
s = time.time()
d = np.bincount(groups, X.sum(axis=1), g_len)
print(time.time() - s)

#compare
print(np.abs((a-d)).sum())

print("use ufunc @Daniel F")
s = time.time()
e = groupby_np(X, groups)
print(time.time() - s)

#compare
print(np.abs((a-e)).sum())

标准输出

original
0.2882847785949707
use scipy coo matrix
0.012301445007324219
0
use scipy csr matrix @Daniel F
0.01046299934387207
0
use bincount @Paul Panzer @Daniel F
0.007468223571777344
0.0
use ufunc @Daniel F
0.04431319236755371
0

胜利者是bincount方案。但csr矩阵方案也非常有趣。
2
仍然比np.bincount(groups, X.sum(axis=1), g_len)慢；-D 不过，稀疏矩阵是一个好主意（+1） - Paul Panzer
1
快速，除了 .sum 不太灵活。 - Daniel F
6

@klim提出的稀疏矩阵解决方案乍一看似乎与求和有关。然而，我们可以通过在csr和csc格式之间转换来在一般情况下使用它：

让我们看一个小例子：

>>> m, n = 3, 8                                                                                                     
>>> idx = np.random.randint(0, m, (n,))
>>> data = np.arange(n)
>>>                                                                                                                 
>>> M = sparse.csr_matrix((data, idx, np.arange(n+1)), (n, m))                                                      
>>>                                                                                                                 
>>> idx                                                                                                             
array([0, 2, 2, 1, 1, 2, 2, 0])                                                                                     
>>> 
>>> M = M.tocsc()
>>> 
>>> M.indptr, M.indices
(array([0, 2, 4, 8], dtype=int32), array([0, 7, 3, 4, 1, 2, 5, 6], dtype=int32))

在转换后，我们可以看到稀疏矩阵的内部表示将索引分组并排序：

>>> groups = np.split(M.indices, M.indptr[1:-1])
>>> groups
[array([0, 7], dtype=int32), array([3, 4], dtype=int32), array([1, 2, 5, 6], dtype=int32)]
>>>

我们可以使用稳定的 argsort 来得到相同的结果：
>>> np.argsort(idx, kind='mergesort')
array([0, 7, 3, 4, 1, 2, 5, 6])
>>>

但是稀疏矩阵实际上更快，即使我们允许 argsort 使用更快的非稳定算法：

>>> m, n = 1000, 100000
>>> idx = np.random.randint(0, m, (n,))
>>> data = np.arange(n)
>>> 
>>> timeit('sparse.csr_matrix((data, idx, np.arange(n+1)), (n, m)).tocsc()', **kwds)
2.250748165184632
>>> timeit('np.argsort(idx)', **kwds)
5.783584725111723

如果我们要求argsort保持分组排序，差异会更大：
>>> timeit('np.argsort(idx, kind="mergesort")', **kwds)
10.507467685034499
这确实限制了data必须是1维的，对吗？也就是说，在这种情况下，求和的技巧是将输入的维度从2降至1。如果我理解有误，请告诉我。在进行类似于X轴零轴均值的操作时，似乎你只能使用2d。 - Brad Solomon
如果你想在稀疏矩阵 M 中有 data，那么是的。但这并非绝对必要。可以将 data 保持分开（并为 M 的值放置一些虚拟值），然后使用 M.indices 进行高级索引。在这种情况下，data 可以具有更多的维度。 - Paul Panzer
5

如果你想要一个更灵活的groupby实现，可以使用任何numpy的ufunc进行分组：

def groupby_np(X, groups, axis = 0, uf = np.add, out = None, minlength = 0, identity = None):
    if minlength < groups.max() + 1:
        minlength = groups.max() + 1
    if identity is None:
        identity = uf.identity
    i = list(range(X.ndim))
    del i[axis]
    i = tuple(i)
    n = out is None
    if n:
        if identity is None:  # fallback to loops over 0-index for identity
            assert np.all(np.in1d(np.arange(minlength), groups)), "No valid identity for unassinged groups"
            s = [slice(None)] * X.ndim
            for i_ in i:
                s[i_] = 0
            out = np.array([uf.reduce(X[tuple(s)][groups == i]) for i in range(minlength)])
        else:
            out = np.full((minlength,), identity, dtype = X.dtype)
    uf.at(out, groups, uf.reduce(X, i))
    if n:
        return out

groupby_np(X, groups)
array([15, 30])

groupby_np(X, groups, uf = np.multiply)
array([   0, 3024])

groupby_np(X, groups, uf = np.maximum)
array([5, 9])

groupby_np(X, groups, uf = np.minimum)
array([0, 6])
不错！只有一个小问题：初始值应该可以被设置，你不想从零开始 - 比如说 - 乘法。 - Paul Panzer
@PaulPanzer，没错，应该使用np.full和uf.identity。 - Daniel F
1
.identity，嗯？这些ufunc比我想象的要聪明。 - Paul Panzer
@PaulPanzer 没有我期望的那么聪明。maximum和minimum没有有效的“identity”。不得不稍微虚构一下。 - Daniel F
这太棒了。我不想加上这个问题，但是我想这会排除 np.mean，因为它不是一个 ufunc？ - Brad Solomon
1
True。但你总是可以执行 groupby_np(X, groups) / groupby_np(np.ones_like(X), groups)。 - Daniel F
2

可能有比这更快的方法（现在两个操作数都在制作副本），但是：

np.bincount(np.broadcast_to(groups, X.T.shape).ravel(), X.T.ravel())

array([ 15.,  30.])
0

如果您想将答案扩展到 ndarray，并且仍然希望进行快速计算，则可以扩展 Daniel 的解决方案：

x_len = 500000
g_len = 200
y_len = 2

X = np.arange(x_len * y_len).reshape(x_len, y_len)
groups = np.random.randint(0, g_len, x_len)

# original
a = np.array([X[groups==i].sum(axis=0) for i in np.unique(groups)])

# alternative
bins = [0] + list(np.bincount(groups, minlength=g_len).cumsum())
Z = np.argsort(groups)
d = np.array([X.take(Z[bins[i]:bins[i+1]],0).sum(axis=0) for i in range(g_len)])

在这个例子中，相比原来的方式，创建bins花费了15ms，求和也花费了15ms，总共大约花费了30ms，而不是280ms。

d.shape
>>> (1000, 2)
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接