在物理上转置一个大的非方数值矩阵

Question

有没有比array.transpose.copy()更快的方法来物理转置一个大的2D numpy矩阵？有没有什么例程可以有效地使用内存呢？

Answer 1

我假设，如果行在内存中是连续的，并且没有足够的可用内存来制作副本，那么您需要执行一个更有效地使用CPU缓存的逐行操作。

维基百科上有一篇关于in-place matrix transposition的文章。事实证明，这样的换位是不平凡的。下面是一个跟随循环的算法，如下所述：

import numpy as np
from numba import njit

@njit # comment this line for debugging
def transpose_in_place(a):
    """In-place matrix transposition for a rectangular matrix.
    
    https://stackoverflow.com/a/62507342/6228891 
    
    Parameter:
    
    - a: 2D array. Unless it's a square matrix, it will be scrambled
      in the process.
    
    Return:
        
    - transposed array, using the same in-memory data storage as the 
      input array.
      
    This algorithm is typically 10x slower than a.T.copy().
    Only use it if you are short on memory.
    """
    if a.shape == (1, 1):
        return a # special case

    n, m = a.shape
    
    # find max length L of permutation cycle by starting at a[0,1].
    # k is the index in the flat buffer; i, j are the indices in
    # a.
    L = 0
    k = 1
    while True:
        j = k % m
        i = k // m
        k = n*j + i
        L += 1
        if k == 1:
            break
    permut = np.zeros(L, dtype=np.int32)

    # Now do the permutations, one cycle at a time    
    seen = np.full(n*m, False)
    aflat = a.reshape(-1) # flat view
    
    for k0 in range(1, n*m-1):
        if seen[k0]:
            continue
        # construct cycle
        k = k0
        permut[0] = k0
        q = 1 # size of permutation array
        while True:
            seen[k] = True
            # note that this is slightly faster than the formula
            # on Wikipedia, k = n*k % (n*m-1)      
            i = k // m
            j = k - i*m            
            k = n*j + i
            if k == k0:
                break
            permut[q] = k
            q += 1
            
        # apply cyclic permutation
        tmp = aflat[permut[q-1]]
        aflat[permut[1:q]] = aflat[permut[:q-1]]
        aflat[permut[0]] = tmp            
    
    aT = aflat.reshape(m, n)
    return aT

def test_transpose(n, m):
    a = np.arange(n*m).reshape(n, m)
    aT = a.T.copy()
    assert np.all(transpose_in_place(a) == aT)

def roundtrip_inplace(a):
    a = transpose_in_place(a)
    a = transpose_in_place(a)
    
def roundtrip_copy(a):
    a = a.T.copy()
    a = a.T.copy()

if __name__ == '__main__':
    test_transpose(1, 1)
    test_transpose(3, 4)
    test_transpose(5, 5)
    test_transpose(1, 5)
    test_transpose(5, 1)
    test_transpose(19, 29)

即使我在这里使用numba.njit来编译转置函数中的循环，它仍然比复制转置慢得多。

n, m = 1000, 10000
a_big = np.arange(n*m, dtype=np.float64).reshape(n, m)

%timeit -r2 -n10 roundtrip_copy(a_big)
54.5 ms ± 153 µs per loop (mean ± std. dev. of 2 runs, 10 loops each)

%timeit -r2 -n1 roundtrip_inplace(a_big)
614 ms ± 141 ms per loop (mean ± std. dev. of 2 runs, 1 loop each) 

Answer 2

无论您做什么，都需要O(n^2)时间和内存。我假设对于您的应用程序，.transpose和.copy (用C编写)将是最有效的选择。

编辑:假设您确实需要复制矩阵

Answer 3

也许有必要看看转置是做什么的，这样我们才能清楚你所说的“物理转置”是什么意思。

从一个小的(4,3)数组开始：

In [51]: arr = np.array([[1,2,3],[10,11,12],[22,23,24],[30,32,34]])             
In [52]: arr                                                                    
Out[52]: 
array([[ 1,  2,  3],
       [10, 11, 12],
       [22, 23, 24],
       [30, 32, 34]])

这与一维数据缓冲区一起存储，我们可以使用ravel显示该缓冲区

In [53]: arr.ravel()                                                            
Out[53]: array([ 1,  2,  3, 10, 11, 12, 22, 23, 24, 30, 32, 34])

和strides，它告诉它按8字节步进列，按24 (3*8)步进行数：

In [54]: arr.strides                                                            
Out[54]: (24, 8)

我们可以用"F“的顺序拉开--这是顺行进行的：

In [55]: arr.ravel(order='F')                                                   
Out[55]: array([ 1, 10, 22, 30,  2, 11, 23, 32,  3, 12, 24, 34])

53是view，55是副本。

现在转置：

In [57]: arrt=arr.T                                                             
In [58]: arrt                                                                   
Out[58]: 
array([[ 1, 10, 22, 30],
       [ 2, 11, 23, 32],
       [ 3, 12, 24, 34]])

这是一个view；我们可以遍历53数据缓冲区，以8字节的步长向下移动。使用arrt进行计算基本上和使用arr一样快。使用strided迭代，订单'F‘和订单'C’的速度一样快。

In [59]: arrt.strides                                                           
Out[59]: (8, 24)

原始顺序：

In [60]: arrt.ravel(order='F')                                                  
Out[60]: array([ 1,  2,  3, 10, 11, 12, 22, 23, 24, 30, 32, 34])

但是做一个'C‘ravel会创建一个副本，和55一样

In [61]: arrt.ravel(order='C')                                                  
Out[61]: array([ 1, 10, 22, 30,  2, 11, 23, 32,  3, 12, 24, 34])

复制转置将生成一个转置为'C‘顺序的数组。这是你的“物理转置”：

In [62]: arrc = arrt.copy()                                                     
In [63]: arrc.strides                                                           
Out[63]: (32, 8)

像61一样重塑转置确实会生成副本，但通常我们不需要显式地生成副本。我认为这样做的唯一原因是为了避免在以后的计算中出现几个多余的副本。