如何反转在​​numpy的置换阵列(How to invert a permutation arra

2019-06-26 13:16发布

鉴于自索引(不知道这是正确的术语)numpy的阵列,例如:

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

这代表本置换 ( =>是箭头):

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

我试图使代表逆变换的阵列,而不在Python“手动”做这件事,那就是我想要一个纯粹的 numpy的解决方案。 我希望在上述情况的结果是:

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

这相当于

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

这似乎很简单,但我想不出如何做到这一点。 我试过google搜索,但没有发现任何有关。

Answer 1:

的置换的逆pnp.arange(n)是索引的阵列s那种p ,即

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])


Answer 2:

排序是这里矫枉过正。 这仅仅是一个单通,具有恒定的存储器需求线性时间算法:

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循环实际上比稍慢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)

从我的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)用于argsortO(n)为单遍算法)和单通算法将是始终如一更快后的足够大的n = p.size (阈值是我的机器上70万左右)。

然而,向量化上述不太直接的方式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解决方案类型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

杰米,安德里斯和保罗在低于花式索引的性能问题就解决了评论指出。 杰米说,它在与NumPy 1.9已经解决了。 我使用Python 3.5和1.11 NumPy的测试它,我用回在2014年的机器上。

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(p.size, p.dtype)
    s[p] = np.arange(p.size)
    return s

方法代码清晰。 在我看来,它是小于晦涩argsort ,也快于大尺寸的输入。 如果速度成为一个问题,我会去用用Cython解决方案。



Answer 3:

我想提供一点点更多的背景larsmans正确答案。 之所以 argsort是正确的可以当您使用的表述中找到由矩阵排列 。 数学优势的置换矩阵 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]]

给定一个置换矩阵,我们可以利用“撤销” multipication由它的逆P^-1 。 置换矩阵的美妙之处在于它们是正交的,因此P*P^(-1)=I ,或换言之P(-1)=P^T ,逆是转置。 这意味着我们可以把转置矩阵的索引来找到你倒置换矢量:

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

其中,如果你想想看,是完全一样的发现的列进行排序的指数P



文章来源: How to invert a permutation array in numpy