numpy数组上的有效循环

Posted

技术标签:

【中文标题】numpy数组上的有效循环【英文标题】:efficient loop over numpy array 【发布时间】:2017-01-15 05:03:35 【问题描述】:

已经问过这个问题的不同版本,但我没有找到满意的答案。

问题:给定一个大的 numpy 向量,找到重复的向量元素的索引(其变化可以与容差进行比较)。

所以问题是〜O(N ^ 2)和内存限制(至少从当前算法的角度来看)。我想知道为什么我尝试过的 Python 比等效的 C 代码慢 100 倍或更多。

import numpy as np
N = 10000
vect = np.arange(float(N))
vect[N/2] = 1
vect[N/4] = 1
dupl = []
print("init done")
counter = 0
for i in range(N):
    for j in range(i+1, N):
        if vect[i] == vect[j]:
            dupl.append(j)
            counter += 1

print("counter =", counter)
print(dupl)
# For simplicity, this code ignores repeated indices 
# which can be trimmed later. Ref output is
# counter = 3
# [2500, 5000, 5000]

我尝试使用 numpy 迭代器,但它们更糟糕(~ x4-5) http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html

使用 N=10,000 我在 C 中得到 0.1 秒,在 Python 中得到 12 秒(上面的代码),在 Python 中使用 np.nditer 得到 40 秒,在 Python 中使用 np.ndindex 得到 50 秒。我把它推到 N=160,000 并且时间按预期缩放为 N^2。

【问题讨论】:

因为 Python 很慢? numpy 数组在使用内置 numpy 函数(在 C 中实现)时非常有效。无论您是否使用 numpy,Python 循环都很慢。尝试仅使用 numpy 函数来实现您的算法。使用内置 Python 函数和/或推导式还应该提高性能(低于 numpy 但高于普通循环)。 因此,Python 中的循环也不错。无论如何,循环有什么困难。我怀疑它的嵌套循环正在杀死这段代码(创建另一个上下文?) 我不想回答我自己的问题,但我最终通过求助于 Numba 解决了。一开始还记得,后来忘记了。时间几乎就是 C 编译代码给我的时间,而且还有一个 Python 循环。所以外循环不是问题。我感谢所有关于使用库的 cmets。实际上,必须使用它们才能获得最佳性能。但是,我发现很难记住所有这些调用,并且从代码而不是库的角度思考对我来说更容易。 nditer 页面以 cython 示例结尾。那是你获得了一些速度。否则nditer 只是处理多输入和输出广播的一种方式。 【参考方案1】:

方法#1

您可以使用triangular matrix 模拟矢量化解决方案的迭代器依赖条件。这是基于处理涉及iterator dependency 的乘法的this post。为了执行vect 中每个元素对其所有元素的元素相等性,我们可以使用NumPy broadcasting。最后,我们可以使用np.count_nonzero 来获取计数,因为它应该在对布尔数组求和时非常有效。

所以,我们会有这样的解决方案 -

mask = np.triu(vect[:,None] == vect,1)
counter = np.count_nonzero(mask)
dupl = np.where(mask)[1]

如果您只关心counter 的计数,我们可以提供以下两种方法。

方法 #2

我们可以避免使用三角矩阵,简单地获取整个计数,然后从对角线元素中减去贡献,并通过将剩余计数减半来仅考虑上三角区域的下一个区域,因为任何一个的贡献都会完全相同。

所以,我们会有一个修改后的解决方案,像这样 -

counter = (np.count_nonzero(vect[:,None] == vect) - vect.size)//2

方法#3

这是一种完全不同的方法,它利用每个唯一元素的计数对最终总数的累积贡献这一事实。

因此,考虑到这个想法,我们将采用第三种方法 -

count = np.bincount(vect) # OR np.unique(vect,return_counts=True)[1]
idx = count[count>1]
id_arr = np.ones(idx.sum(),dtype=int)
id_arr[0] = 0
id_arr[idx[:-1].cumsum()] = -idx[:-1]+1
counter = np.sum(id_arr.cumsum())

【讨论】:

【参考方案2】:

Python 本身是一种高度动态、缓慢的语言。 numpy 中的想法是使用vectorization,并避免显式循环。在这种情况下,您可以使用np.equal.outer。你可以从

a = np.equal.outer(vect, vect)

现在,例如,求和:

 >>> np.sum(a)
 10006

要找到相等的 i 索引,您可以这样做

np.fill_diagonal(a, 0)

>>> np.nonzero(np.any(a, axis=0))[0]
array([   1, 2500, 5000])

时机

def find_vec():
    a = np.equal.outer(vect, vect)
    s = np.sum(a)
    np.fill_diagonal(a, 0)
    return np.sum(a), np.nonzero(np.any(a, axis=0))[0]

>>> %timeit find_vec()
1 loops, best of 3: 214 ms per loop

def find_loop():
    dupl = []
    counter = 0
    for i in range(N):
        for j in range(i+1, N):
             if vect[i] == vect[j]:
                 dupl.append(j)
                 counter += 1
    return dupl

>>> % timeit find_loop()
1 loops, best of 3: 8.51 s per loop

【讨论】:

使用这种方法的加速增益是多少? @vz0 相对于循环? OP 在他的问题中发布了运行纯 python 版本需要 12 秒的问题。运行您的版本需要多长时间? @vz0 查看编辑 - 我的循环大约需要 8.5 秒,而矢量化版本需要 214 毫秒。 矢量化版本还是 O(N**2) 吗?【参考方案3】:

我想知道为什么我尝试过的 Python 比等效的 C 代码慢 100 倍或更多。

因为 Python 程序通常比 C 程序慢 100 倍。

您可以在 C 中实现关键代码路径并提供 Python-C 绑定,或者更改算法。您可以使用 dict 将数组从值反转为索引来编写 O(N) 版本。

import numpy as np
N = 10000
vect = np.arange(float(N))
vect[N/2] = 1
vect[N/4] = 1
dupl = 
print("init done")
counter = 0
for i in range(N):
    e = dupl.get(vect[i], None)
    if e is None:
        dupl[vect[i]] = [i]
    else:
        e.append(i)
        counter += 1

print("counter =", counter)
print([(k, v) for k, v in dupl.items() if len(v) > 1])

编辑:

如果您需要使用 abs(vect[i] - vect[j])

abs(vect[i] - vect[j]) < eps ->
abs(vect[i] - vect[j]) / eps < (eps / eps) ->
abs(vect[i]/eps - vect[j]/eps) < 1
int(abs(vect[i]/eps - vect[j]/eps)) = 0

像这样:

import numpy as np
N = 10000
vect = np.arange(float(N))
vect[N/2] = 1
vect[N/4] = 1
dupl = 
print("init done")
counter = 0
eps = 0.01
for i in range(N):
    k = int(vect[i] / eps)
    e = dupl.get(k, None)
    if e is None:
        dupl[k] = [i]
    else:
        e.append(i)
        counter += 1

print("counter =", counter)
print([(k, v) for k, v in dupl.items() if len(v) > 1])

【讨论】:

不幸的是,字典方法对我不起作用,因为实际上我需要在一定的公差范围内比较数组元素 abs(vect[i] - vect[j]) 然后您可以将 dict 中的值标准化为 eps。 abs(vect[i] - vect[j]) abs(vect[i] - vect[j]) / eps abs(vect[i]/eps - vect[j]/每股收益) 【参考方案4】:

这个使用numpy_indexed包的解决方案复杂度为n Log n,并且是完全向量化的;所以很可能与 C 的性能没有太大区别。

import numpy_indexed as npi
dpl = np.flatnonzero(npi.multiplicity(vect) > 1)

【讨论】:

【参考方案5】:

显而易见的问题是您为什么要以这种方式执行此操作。 NumPy 数组旨在成为不透明的数据结构——我的意思是 NumPy 数组旨在在 NumPy 系统内部创建,然后将操作发送到 NumPy 子系统以传递结果。即 NumPy 应该是一个黑盒子,您可以向其中抛出请求并输出结果。

所以鉴于上面的代码,我一点也不惊讶 NumPy 的性能比可怕的还要糟糕。

我相信,以下内容应该是您想要的,但使用 NumPy 方式完成:

import numpy as np

N = 10000
vect = np.arange(float(N))
vect[N/2] = 1
vect[N/4] = 1

print([np.where(a == vect)[0] for a in vect][1])

# Delivers [1, 2500, 5000]

【讨论】:

不错的一个。我想到了 np.where 并尝试了它。问题是,我不需要答案中的 1,只需要 2500 和 5000,即重复。其次,近似相似的泛化不是那么整齐而且有点慢。但是,是的,你的一个班轮和 C 一样快。【参考方案6】:

作为 Ami Tavory 答案的替代方案,您可以使用 collections 包中的 Counter 来检测重复项。在我的电脑上,它似乎更快。请参阅下面的函数,该函数也可以找到不同的重复项。

import collections
import numpy as np

def find_duplicates_original(x):
    d = []
    for i in range(len(x)):
        for j in range(i + 1, len(x)):
            if x[i] == x[j]:
                d.append(j)
    return d

def find_duplicates_outer(x):
    a = np.equal.outer(x, x)
    np.fill_diagonal(a, 0)
    return np.flatnonzero(np.any(a, axis=0))

def find_duplicates_counter(x):
    counter = collections.Counter(x)
    values = (v for v, c in counter.items() if c > 1)
    return v: np.flatnonzero(x == v) for v in values


n = 10000
x = np.arange(float(n))
x[n // 2] = 1
x[n // 4] = 1

>>>> find_duplicates_counter(x)
1.0: array([   1, 2500, 5000], dtype=int64)

>>>> %timeit find_duplicates_original(x)
1 loop, best of 3: 12 s per loop

>>>> %timeit find_duplicates_outer(x)
10 loops, best of 3: 84.3 ms per loop

>>>> %timeit find_duplicates_counter(x)
1000 loops, best of 3: 1.63 ms per loop

【讨论】:

不幸的是,这与 Ami Tavory 的建议不适合。试试 N=160000【参考方案7】:

与您的代码 18 秒相比,这在 8 毫秒内运行,并且不使用任何奇怪的库。它类似于@vs0 的方法,但我更喜欢defaultdict。它应该是大约 O(N)。

from collections import defaultdict
dupl = []
counter = 0
indexes = defaultdict(list)
for i, e in enumerate(vect):
    indexes[e].append(i)
    if len(indexes[e]) > 1:
        dupl.append(i)
        counter += 1

【讨论】:

【参考方案8】:

由于答案已经停止出现并且没有一个是完全令人满意的,为了记录,我发布了我自己的解决方案。

我的理解是,在这种情况下,让 Python 变慢的是赋值,而不是我最初认为的嵌套循环。使用库或编译代码消除了分配的需要,性能显着提高。

from __future__ import print_function
import numpy as np
from numba import jit

N = 10000
vect = np.arange(N, dtype=np.float32)

vect[N/2] = 1
vect[N/4] = 1
dupl = np.zeros(N, dtype=np.int32)

print("init done")
# uncomment to enable compiled function
#@jit
def duplicates(i, counter, dupl, vect):
    eps = 0.01
    ns = len(vect)
    for j in range(i+1, ns):
        # replace if to use approx comparison
        #if abs(vect[i] - vect[j]) < eps:
        if vect[i] == vect[j]:
            dupl[counter] = j
            counter += 1
    return counter

counter = 0
for i in xrange(N):
    counter = duplicates(i, counter, dupl, vect)

print("counter =", counter)
print(dupl[0:counter])

测试

# no jit
$ time python array-test-numba.py
init done
counter = 3
[2500 5000 5000]

elapsed 10.135 s

# with jit
$ time python array-test-numba.py
init done
counter = 3
[2500 5000 5000]

elapsed 0.480 s

编译版本(@jit 未注释)的性能接近 C 代码性能 ~0.1 - 0.2 秒。也许消除最后一个循环可以进一步提高性能。使用 eps 进行近似比较时,性能差异甚至更大,而编译版本的差异很小。

# no jit
$ time python array-test-numba.py
init done
counter = 3
[2500 5000 5000]

elapsed 109.218 s

# with jit
$ time python array-test-numba.py
init done
counter = 3
[2500 5000 5000]

elapsed 0.506 s

这是 ~ 200 倍的差异。在实际代码中,我必须将两个循环都放在函数中,并使用具有变量类型的函数模板,所以它有点复杂但不是很复杂。

【讨论】:

以上是关于numpy数组上的有效循环的主要内容,如果未能解决你的问题,请参考以下文章

有效地减去不同形状的numpy数组

附加numpy数组的有效方法

如何有效地操作一个大的numpy数组

使用 numpy 数组和共享内存并行化 python 循环

一维 NumPy 数组上的滑动标准差

在 NumPy 数组的每个单元格处对函数进行有效评估