如何在numpy中反转排列数组

Posted 2023-02-23

【中文标题】如何在numpy中反转排列数组

【英文标题】：How to invert a permutation array in numpy

【发布时间】：2012-07-23 21:00:51

【问题描述】：

给定一个自索引（不确定这是否是正确的术语）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

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

【问题讨论】：

a = np.array([1, 1, 1, 1])应该返回什么？

@eumiro 你可以假设这种情况不会出现。

【参考方案1】：

在这里排序有点过头了。这只是一个单通道线性时间算法，需要恒定的内存：

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 快；有趣的是，上述for 循环的微不足道的矢量化（s[p] = xrange(p.size)，see index arrays）实际上比np.argsort 稍慢，只要p.size < 700 000（嗯，在我的机器上，你的里程会 em> 不同）：

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)

来自我的 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) 用于单通道算法）并且在足够大的 n = p.size 之后单通道算法将始终更快（阈值约为700k 在我的机器上）。

但是，使用np.put 对上述for 循环进行矢量化还有一种不太直接的方法：

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 对于较小的n 仍然优于np.put 方法（我的机器上的临界点在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_put 方法分配并填充了一个额外的数组（在np.arange() 调用中）。

---

虽然你没有要求 Cython 解决方案，只是出于好奇，我还用typed memoryviews 计时了以下 Cython 解决方案：

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 在下面的 cmets 中指出，花式索引的性能问题已得到解决。 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

确实是一个显着的改进！

---

---

结论

总而言之，我会选择

def invert_permutation(p):
    '''The argument p is assumed to be some permutation of 0, 1, ..., len(p)-1. 
    Returns an array s, where s[i] gives the index of i in p.
    '''
    s = np.empty_like(p)
    s[p] = np.arange(p.size)
    return s

代码清晰的方法。在我看来，它不像argsort 那样晦涩难懂，而且对于大输入大小也更快。如果速度成为问题，我会选择 Cython 解决方案。

【讨论】：

不错！我不知道np.put，它填补了我的慢速解决方案和 Cython 之间的重要利基。 +1。

+1 大智若愚！ ;-) 大约在你为一个两年前的问题写这个答案的同时，我发送了一个 PR 以在 numpy 的unique 函数中使用与此非常相似的技术，请参阅here。 FWIW，np.put 相对于花哨的索引的优势在 numpy 1.9 中大部分都消失了..

@Jaime 谢谢你的好消息！我发现花哨的索引是最干净的解决方案；很遗憾它曾经如此缓慢。很高兴知道它在 NumPy 1.9 中大部分都消失了。

s[p]=np.arange(p.size) 更不晦涩，在我的机器上运行速度是 np.put 的两倍（我知道，我知道）。

@Paul 感谢您的信息！事实上，显然从 NumPy 1.9 开始，使用 np.put() 就没有意义了。我会很快相应地更新我的答案！

【参考方案2】：

np.arange(n) 的排列 p 的倒数是排序 p 的索引数组 s，即

p[s] == np.arange(n)

必须全部为真。这样的s 正是np.argsort 返回的内容：

>>> p = np.array([3, 2, 0, 1])
>>> np.argsort(p)
array([2, 3, 1, 0])
>>> p[np.argsort(p)]
array([0, 1, 2, 3])

【讨论】：

@lazyr：解释一下。

@larsmans 还有一个更简单的单通算法：任务基本是s[p] = xrange(p.size)，请看我的回答。

【参考方案3】：

我想为 larsman 的正确答案提供更多背景知识。当您使用permutation by a matrix 的表示时，可以找到argsort 正确的原因。置换 matrix P 的数学优势在于矩阵“对向量进行运算”，即置换矩阵乘以向量可以置换向量。

你的排列看起来像：

import numpy as np
a   = np.array([3,2,0,1])
N   = a.size
rows = np.arange(N)
P   = np.zeros((N,N),dtype=int)
P[rows,a] = 1

[[0 0 0 1]
 [0 0 1 0]
 [1 0 0 0]
 [0 1 0 0]]

给定一个置换矩阵，我们可以通过乘以它的逆 P^-1 来“撤销”乘法。置换矩阵的美妙之处在于它们是正交的，因此P*P^(-1)=I，或者换句话说P(-1)=P^T，逆是转置。这意味着我们可以利用转置矩阵的索引来找到你的倒排向量：

inv_a = np.where(P.T)[1]
[2 3 1 0]

如果您考虑一下，这与查找对P 的列进行排序的索引完全相同！

【讨论】：

非常感谢您的解释！很有启发性。