ajz*_*z34 5 python parallel-processing multithreading transpose numpy
由于 NumPy 的转置函数，我目前遇到了巨大的开销。我发现这个函数实际上总是在单线程中运行，无论转置矩阵/数组有多大。我可能需要避免这种巨大的时间成本。
\n据我了解，np.dot如果 numpy 数组足够大，其他函数（例如向量增量）将并行运行。一些按元素操作似乎在 numexpr 包中可以更好地并行化，但 numexpr 可能无法处理转置。
我想了解解决问题的更好方法。为了详细说明这个问题，
\nB = A.T），因为转置后的张量不用于计算或转储，并且在此阶段不需要真正转置数据。调用时B[:] = A.T，确实会转置数据。在此先感谢您的任何建议！
\n在我的 Linux 个人计算机上处理模型转置问题（大小约为A763MB），可用 4 核（总共 400% CPU）。
A = np.random.random((9999, 10001))\nB = np.random.random((10001, 9999))\nD = np.random.random((9999, 10001))\n当前的解决方法似乎不够有效。一般来说，如果在 4 核 CPU 上完全并行化，它应该会看到大约 3 倍~4 倍的加速，但我编写的代码只获得了大约 1.5 倍。
\nB[:] = A.T # ~255 ms, time for transpose\nD[:] = A[:] # ~ 65 ms, time for coping data\n有趣的是，如果A是 10000 * 10000 方阵，那么转置时间将增加到约 310ms。我不知道这里会发生什么。如果矩阵是方阵，甚至 C/ctypes 绑定函数的时间成本也会受到影响（更慢）。
使用以下 OpenMP/BLAS（基本，未优化）编写：
\nA = np.random.random((9999, 10001))\nB = np.random.random((10001, 9999))\nD = np.random.random((9999, 10001))\n然后执行python代码（4核线程）
\nB[:] = A.T # ~255 ms, time for transpose\nD[:] = A[:] # ~ 65 ms, time for coping data\n使用 C/ctype 绑定可能会有点麻烦，因为它不是纯 python，并且也使用 CBLAS 作为外部包。
\n// transpose.c\n#include <stdio.h>\n#include <cblas.h>\nvoid matrix_transpose(double* A, double* B, size_t d1, size_t d2) {\n size_t i;\n #pragma omp parallel for shared(A, B, d1, d2) private(i)\n for (i = 0; i < d1; ++i) {\n cblas_dcopy(d2, A + i*d2, 1, B + i, d1);\n }\n}\n我猜想对于大型集群（HPC）来说，更多的进程/线程并不一定会获得加速。那么如何设置进程/线程的数量也可能是一个问题。
\n这个问题不仅与并行化有关，还与缓存感知和基于分片的转置算法有关。可能的解决方案可能是
\npyscf.lib.transpose( python link , c link ) 进行并行 3 索引张量转置M.transpose(0,2,1)。我是pyscf的粉丝。它的主要目的是量子或半经验化学计算，但它的一些数值计算实用程序可能为此目的进行了优化。在我在服务器上的测试中，转置 (1502, 601, 601) 张量可能比 MKL mkl_domatcopy（9.19 秒）快两倍（4.09 秒）。相关算法文章：
\n\n相关 stackoverflow 和 github 问题页面：
\n\n有效地计算转置是很困难的。该原语不受计算限制，而是受内存限制。对于存储在 RAM（而不是 CPU 缓存）中的大矩阵尤其如此。
Numpy 使用基于视图的方法，当只需要数组的一部分并且计算速度非常慢时（例如，当执行复制时），这种方法非常有用。在这种情况下执行复制时，Numpy 的实现方式会导致大量缓存未命中，从而严重降低性能。
我发现这个函数实际上总是在单线程中运行，无论转置矩阵/数组有多大。
这显然取决于所使用的 Numpy 实现。AFAIK，一些优化的软件包（例如英特尔的软件包）效率更高，并且更经常地利用多线程。
我认为并行转置函数应该是一个解决方案。问题是如何实施。
是和不是。使用更多线程可能并不需要更快。至少不多，而且不是在所有平台上。使用的算法远比使用多线程重要。
在现代桌面 x86-64 处理器上，每个内核都可以受其缓存层次结构的限制。但这个限制是相当高的。因此，两个核心通常足以使 RAM 吞吐量接近饱和。例如，在我的（4 核）机器上，顺序复制可以达到 20.4 GiB/s（Numpy 成功达到此限制），而我的（实际）内存吞吐量接近 35 GiB/s。复制A需要 72 毫秒，而简单的 Numpy 转置A需要 700 毫秒。即使使用我的所有内核，并行实现也不会快于 175 毫秒，而最佳理论时间为 42 毫秒。实际上，由于缓存未命中和 L3 缓存饱和，简单的并行实现会比 175 毫秒慢得多。
简单的转置实现不会连续写入/读取数据。内存访问模式是跨步的，并且大多数高速缓存行都被浪费了。因此，数据会在巨大的矩阵上从内存中读取/写入多次（在当前使用双精度的 x86-64 平台上通常为 8 次）。基于图块的转置算法是防止此问题的有效方法。它还大大减少了缓存未命中的情况。理想情况下，应使用分层切片或忽略缓存的 Z 切片模式，尽管这很难实现。
以下是基于之前信息的基于 Numba 的实现：
@nb.njit('void(float64[:,::1], float64[:,::1])', parallel=True)
def transpose(mat, out):
blockSize, tileSize = 256, 32 # To be tuned
n, m = mat.shape
assert blockSize % tileSize == 0
for tmp in nb.prange((m+blockSize-1)//blockSize):
i = tmp * blockSize
for j in range(0, n, blockSize):
tiMin, tiMax = i, min(i+blockSize, m)
tjMin, tjMax = j, min(j+blockSize, n)
for ti in range(tiMin, tiMax, tileSize):
for tj in range(tjMin, tjMax, tileSize):
out[ti:ti+tileSize, tj:tj+tileSize] = mat[tj:tj+tileSize, ti:ti+tileSize].T
如果您想要更快的代码，您可以使用非常优化的本机库来实现转置，例如英特尔 MKL。此类库通常利用低级处理器特定指令（SIMD 指令和非临时存储）来更有效地使用缓存/RAM。
以下是计时结果（假设输出矩阵已填充到内存中）：
Naive Numpy: 700 ms
Above code without Numba (sequential): 287 ms
Numba version (sequential): 157 ms
Numba version (parallel): 104 ms
Very-optimized C code (parallel): 68 ms
Theoretical optimum: 42 ms