将Numpy数组以FORTRAN顺序转换为ctypes

Question

将Numpy数组以FORTRAN顺序转换为ctypes

9

有没有高效的方法将numpy数组转换为FORTRAN序列化的ctypes数组，最好不需要进行拷贝并避免与步幅相关的问题？

import numpy as np

# Sample data
n = 10000
A = np.zeros((n,n), dtype=np.int8)
A[0,1] = 1

def slow_conversion(A):
    return np.ctypeslib.as_ctypes(np.ascontiguousarray(A.T))

assert slow_conversion(A)[1][0] == 1

仅对调用 as_ctypes 进行性能分析：

%%timeit
np.ctypeslib.as_ctypes(A)

每次循环3.35微秒±10.5纳秒(7次运行的平均值和标准差，每次循环100000次)

所提供（缓慢）转换的性能分析

%%timeit
slow_conversion(A)

每个循环大约需要 206 毫秒 ± 10.4 毫秒（7 次运行的平均值 ± 标准差）

理想结果是获得与仅调用 as_ctypes 相似的性能。

- Andrew Walker

2

如何以及何时进行转置以加速并不是一个简单的问题，而且很大程度上取决于底层硬件。请参考这个类似的问题 https://dev59.com/tGQn5IYBdhLWcg3wj3zz - Carlos Horn

2

你想将一个“C-ordered”的numpy数组转换为FORTRAN ordered的ctypes数组吗？因为如果你已经有一个“F-ordered”的numpy数组，那么这个转换可以很快地完成（请参见我对@Stephan答案的评论）。 - Vladimir Fokow

输入是否为正方形可以假定吗？ - norok2

1

@norok2，无论是正方形还是非正方形情况都是相关的。不过，我仍然很想知道是否有一种特殊情况可以处理输入为正方形的情况。 - Andrew Walker

3个回答

4

如果您的起始数组必须是C序列：

正如您所指出的，np.ctypeslib.as_ctypes(...)很快。

您计算的瓶颈是np.ascontiguousarray(A.T)，它等同于np.asfortranarray(A)，同样很慢。

这让我相信仅使用numpy无法使其更快。我的意思是，既然已经存在一个专门用于此的函数（np.asfortranarray），我们可以假设其已经被优化，以具有最佳可能的性能。

如果您可以从F序列数组开始：

在这种情况下，问题解决了！因为数组已经按照所需的内存顺序排列，所以您可以像as_ctypes调用那样获得性能：

n = 10000
A = np.zeros((n, n), dtype=np.int8, order='F')
A[0, 1] = 1

def conversion(A):
    return np.ctypeslib.as_ctypes(np.ascontiguousarray(A.T))

assert conversion(A)[1][0] == 1

%%timeit
conversion(A)

# 4.76 µs ± 68.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

感谢@JérômeRichard和@StephanSchlecht的启发，这个想法。阅读@JérômeRichard的回答以获取更多见解！

- Vladimir Fokow

2

我只是好奇，Fortran的order数组不是连续数组的转置吗？ - user16836078

@Kevin，@jkr，print(np.asfortranarray(A)) 打印的结果与 print(A) 相同。它们之间的区别在于它们元素在内存中排列的顺序不同。 - Vladimir Fokow

1

这个方面还有待改进。

操作不仅是在复制，而是因为它加载和存储到主内存而不是使用缓存。

通常情况下处理器会访问多个块的多个连续字节，然后从缓存中使用它。如果缓存空间不足，则会删除一些旧的块。

为了论证，假设您的CPU每个块包含8个字节，并且行是相邻的。在一个矩阵中，您将访问列，在另一个矩阵中您将访问行。当您写入一列时，您正在加载多个列但只更新一个。通过复制几列可以看出这个单列的开销。

n = 2**14
A = np.random.randint(0, 100, (n,n), dtype=np.int8)
B = np.empty_like(A)
%%timeit
B[:1,:] = A[:1,:]
%%timeit
B[:4,:] = A[:4,:]

如果你在行上做同样的操作，你应该会发现大致呈线性。如果你复制列，复制一列的成本非常接近于复制两列甚至8列或16列，这取决于硬件。

我将使用 n=2**14 来简化事情，但这个原则适用于任何维度。

如果你有足够小的矩阵，比如 8 x 8，整个矩阵都可以放入缓存中，因此你可以在不访问任何缓存的情况下转置它。
如果你正在复制大块连续的数据，即使你不能在缓存中完成整个操作，也可以减少给定数据从/到内存加载的次数。

基于这个，我尝试重新排列矩阵，将其变成较小的连续块的矩阵，首先我转置一个块中的元素，然后转置矩阵中的块。

对于基准测试：

B = np.ascontiguousarray(A.T)

3.12 s ± 446 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

使用8x8块

T0 = A.reshape(2048,8,2048,8)
T1 = np.ascontiguousarray(T0.transpose(0,2,3,1))
T2 = np.ascontiguousarray(T1.transpose(1,0,2,3))
T3 = np.ascontiguousarray(T2.transpose(0,2,1,3))
B = T3.reshape(A.shape)

786 ms ± 54.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

assert np.all(B == A.T) # 2.8s

它仍然比简单复制慢200倍，但已经比原始方法快4倍。

只分配两个临时数组而不是3个也有帮助。

T0 = np.empty_like(A)
T1 = np.empty_like(A)
T0.reshape(2048,2048,8,8)[:] = A.reshape(2048,8,2048,8).transpose(0,2,3,1)
T1.reshape(2048,2048,8,8)[:] = T0.reshape(2048,2048,8,8).transpose(1,0,2,3)
T0.reshape(2048,8,2048,8)[:] = T1.reshape(2048,2048,8,8).transpose(0,2,1,3)
B = T0

686 ms ± 60.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

- Bob

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jérôme Richard · Accepted Answer

默认情况下，NumPy创建C-ordered数组（因为它是用C编写的），即行主序数组。A.T执行转置操作时会创建一个视图，其中步幅被反转（即没有复制）。但该数组不再连续，所以np.ascontiguousarray会执行复制，而复制操作很耗费时间。这就是为什么slow_conversion很慢的原因。请注意，可以使用yourarray.flags['F_CONTIGUOUS']和检查yourarray.strides来测试连续性。同时，yourarray.flags和yourarray.__array_interface__提供有关数组是否已复制以及步幅的信息。

根据文档，np.asfortranarray返回按Fortran顺序布局的数组。必要时，它可以执行复制。实际上，np.asfortranarray(A)会执行复制，而np.asfortranarray(A.T)不会。您可以检查该函数的C代码以了解更多信息。由于两者都被视为FORTRAN-contiguous，因此在这种情况下最好使用np.asfortranarray(A.T)，这样不需要进行任何复制。

关于ctypes，它处理按行主序存储的C数组，与按列主序存储的FORTRAN数组不同。此外，C数组不支持步幅，而FORTRAN数组则原生支持步幅。这意味着一行基本上是存储在内存中的连续数据的内存视图。由于需要满足slow_conversion(A)[1][0] == 1，这意味着返回的对象应该具有第二行的第一个项为1，并且因此列必须在内存中强制连续存储。问题在于初始数组不是FORTRAN-contiguous，而是C-contiguous，因此需要进行转置。转置操作非常耗费时间（尽管NumPy实现不太优化）。

假设您不想支付复制/转置的开销，则需要放宽要求。有几种可能的选择：

使用Numpy直接创建FORTRAN有序数组，例如 np.zeros((n,n), dtype=np.int8, order='F')。这会创建一个C数组，其步幅已经转置，以便像FORTRAN数组一样快速处理操作列（请记住，Numpy是用C编写的，因此按行主顺序排序的数组是参考）。因此，在ctypes中，第一行实际上是一列。请注意，在性能方面混合使用C和FORTRAN有序数组时，应该非常小心，因为非连续访问会慢得多。
使用跨步的FORTRAN数组。这种解决方案基本上意味着基本的基于列的计算将会慢得多，因此需要编写几乎不寻常的基于行的计算方法。您可以使用A.ctypes.data_as(POINTER(c_double))提取C连续数组的指针，使用A.strides提取步幅，使用A.shape提取形状。也就是说，这种解决方案似乎并不是真正可移植/标准化的。标准方法似乎是在FORTRAN中使用C绑定。我对此不太熟悉，但您可以在此答案中找到完整的示例。

最后一种解决方案是手动使用快速转置算法转置数据。原地转置比就地转置更快（通常是2倍），但这需要一个平方数组，并且它会改变不应该再使用的输入数组（除非对于操作转置的数组没问题）。不能直接使用Numpy进行快速转置。一种解决方案是在Numba中执行它（这是一种次优的就地解决方案，但已经相当快了），或者在C中执行它或者直接在FORTRAN中执行它（在所有情况下都使用包装函数）。这应该比Numpy做的要快得多，但仍远远慢于基本的ctypes封装。