如何在numpy中反转排列数组

7cjasjjr  于 2023-02-08  发布在  其他
关注(0)|答案(3)|浏览(133)

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

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

这表示该排列(=>是箭头):

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

这看起来很简单,但我就是想不出该怎么做。我试着用谷歌搜索,但没有找到任何相关的东西。

deyfvvtc

deyfvvtc1#

简短回答

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
    • 排序在这里是一种过度使用。**这只是一个单遍、线性时间算法,需要恒定的内存:
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),参见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)

从我的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

最终,渐近复杂性开始显现(argsortO(n log n)与单遍算法的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倍的速度几乎没有!**

公平地说,对于较小的nnp.argsort仍然优于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.8ms;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

确实是一个重大的进步!

结论

总而言之,为了代码的清晰性,我会采用上面提到的Short answer方法。在我看来,它比argsort更容易理解,而且对于大输入大小也更快。如果速度成为一个问题,我会采用Cython解决方案。

ct2axkht

ct2axkht2#

np.arange(n)的置换p的逆是对p排序的索引s的阵列,即

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

必须全为true。np.argsort返回的就是这样一个s

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

jhiyze9q3#

我想为larsmans的正确答案提供一点背景知识。当你使用permutation by a matrix的表示时,可以找到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的列进行排序的索引完全相同!

相关问题