如何在numpy中反转置换数组

Question

给定自我索引(不确定这是否是正确的术语)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

看起来很简单,但我想不出怎么做.我试过谷歌搜索,但没有找到任何相关的.

Answer 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; 有趣的是,上述循环的平凡向量化(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)用于argsort与O(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解决方案的类型化内存视图:

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 ms; argsort需要72.7 ms).

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

确实有显着改善!

结论

总而言之,我会选择

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 2

的置换的逆p的np.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 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]]

给定置换矩阵,我们可以通过乘以它的逆来"撤销"乘法P^-1.置换矩阵的美妙之处在于它们是正交的P*P^(-1)=I,或者换句话说P(-1)=P^T,逆是转置.这意味着我们可以使用转置矩阵的索引来找到您的反转置换向量:

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

如果你考虑它,那就像找到对列进行排序的索引完全相同P!