python - 如何在numpy中反转排列数组

Question

给定一个自索引(不确定这是否是正确的术语)numpy 数组，例如:

a = np.array([3, 2, 0, 1])

这表示 permutation (=> 是一个箭头):

0 => 3
1 => 2
2 => 0
3 => 1

我正在尝试制作一个表示逆变换的数组，而不是在 python 中“手动”执行它，也就是说，我想要一个 pure numpy 解决方案。在上述情况下我想要的结果是:

array([2, 3, 1, 0])

相当于

0 <= 3                0 => 2
1 <= 2       or       1 => 3
2 <= 0                2 => 1
3 <= 1                3 => 0

看起来很简单，但我就是想不出怎么做。我试过谷歌搜索，但没有找到任何相关信息。

Answer 1

简答

def invert_permutation(p):
    """Return an array s with which np.array_equal(arr[p][s], arr) is True.
    The array_like argument p must be some permutation of 0, 1, ..., len(p)-1.
    """
    p = np.asanyarray(p) # in case p is a tuple, etc.
    s = np.empty_like(p)
    s[p] = np.arange(p.size)
    return s

---

---

在这里排序是多余的。这只是一个单 channel 线性时间算法，具有恒定的内存需求:

from __future__ import print_function
import numpy as np

p = np.array([3, 2, 0, 1])
s = np.empty(p.size, dtype=np.int32)
for i in np.arange(p.size):
    s[p[i]] = i

print('s =', s)

上面的代码打印出来

 s = [2 3 1 0]

根据需要。

答案的其余部分与上述 for 的有效矢量化有关。环形。如果你只是想知道结论，请跳到这个答案的末尾。

---

(2014 年 8 月 27 日的原始答案；时间对 NumPy 1.8 有效。稍后会更新 NumPy 1.11。)

单次线性时间算法预计比 np.argsort 更快;有趣的是，上述 s[p] = xrange(p.size) 的微不足道的矢量化( for ， see index arrays )循环实际上比 np.argsort 稍慢只要p.size < 700 000 (嗯，在我的机器上，您的里程会有所不同):

import numpy as np

def np_argsort(p):
    return np.argsort(p)
    
def np_fancy(p):
    s = np.zeros(p.size, p.dtype) # np.zeros is better than np.empty here, at least on Linux
    s[p] = xrange(p.size) 
    return s

def create_input(n):
    np.random.seed(31)
    indices = np.arange(n, dtype = np.int32)
    return np.random.permutation(indices)

来 self 的 IPython 笔记本:

p = create_input(700000)
%timeit np_argsort(p)
10 loops, best of 3: 72.7 ms per loop
%timeit np_fancy(p)
10 loops, best of 3: 70.2 ms per loop

最终，渐近复杂度开始出现(O(n log n) 用于 argsort 与 O(n) 用于单 channel 算法)，并且在足够大的 n = p.size 之后单 channel 算法将始终更快(我的机器上的阈值约为 700k)。

但是，有一种不太直接的方法可以对上述 for 进行矢量化。循环 np.put :

def np_put(p):
    n = p.size
    s = np.zeros(n, dtype = np.int32)
    i = np.arange(n, dtype = np.int32)
    np.put(s, p, i) # s[p[i]] = i 
    return s

这给出了 n = 700 000 (和上面一样大小):

p = create_input(700000)
%timeit np_put(p)
100 loops, best of 3: 12.8 ms per loop

这是一个不错的 5.6 倍加速，几乎没有!

公平地说，np.argsort仍然胜过np.put较小的方法n (在我的机器上，临界点在 n = 1210 左右):

p = create_input(1210)
%timeit np_argsort(p)
10000 loops, best of 3: 25.1 µs per loop
%timeit np_fancy(p)
10000 loops, best of 3: 118 µs per loop
%timeit np_put(p)
10000 loops, best of 3: 25 µs per loop

这很可能是因为我们使用 np.arange() 分配并填充了一个额外的数组(在 np_put 调用中)接近。

---

虽然您没有要求 Cython 解决方案，但出于好奇，我还将以下 Cython 解决方案计时为 typed memoryviews :

import numpy as np
cimport numpy as np

def in_cython(np.ndarray[np.int32_t] p):    
    cdef int i
    cdef int[:] pmv
    cdef int[:] smv 
    pmv = p
    s = np.empty(p.size, dtype=np.int32)
    smv = s
    for i in xrange(p.size):
        smv[pmv[i]] = i
    return s

时间安排:

p = create_input(700000)
%timeit in_cython(p)
100 loops, best of 3: 2.59 ms per loop

所以，np.put解决方案仍然没有尽可能快(这个输入大小运行了 12.8 毫秒；argsort 花了 72.7 毫秒)。

---

2017 年 2 月 3 日更新 NumPy 1.11

Jamie、Andris 和 Paul 在下面的评论中指出，花式索引的性能问题已得到解决。 Jamie 说它已经在 NumPy 1.9 中解决了。我在 2014 年使用的机器上使用 Python 3.5 和 NumPy 1.11 对其进行了测试。

def invert_permutation(p):
    s = np.empty(p.size, p.dtype)
    s[p] = np.arange(p.size)
    return s

时间安排:

p = create_input(880)
%timeit np_argsort(p)
100000 loops, best of 3: 11.6 µs per loop
%timeit invert_permutation(p)
100000 loops, best of 3: 11.5 µs per loop

确实是一个显着的改进!

---

结论

总而言之，为了代码清晰，我会采用顶部提到的简答方法。在我看来，它没有 argsort 那样晦涩难懂。 ，并且对于大输入大小也更快。如果速度成为问题，我会选择 Cython 解决方案。