在没有 GIL 的情况下在 Cython 中并行化

Question

我正在尝试计算 numpy 数组的某些列，使用 cdef 函数在 for 循环中对 python 对象（numpy 数组）进行操作.

我想并行执行。但不确定如何执行。

这是一个玩具示例，一个 def 函数使用 prange 在 for 循环中调用 cdef 函数，这不是允许，因为 np.ndarray 是一个 python 对象。在我的实际问题中，一个矩阵和一个向量是 cdef 函数的参数，并且执行一些 numpy 矩阵运算，例如 np.linalg.pinv() （我猜这实际上是瓶颈）。

%%cython
import numpy as np
cimport numpy as np
from cython.parallel import prange
from c_functions import estimate_coef_linear_regression

DTYPE = np.float
ctypedef np.float_t DTYPE_t

def transpose_example(np.ndarray[DTYPE_t, ndim=2] data):
    """
    Transposes a matrix. It does each row independently and parallel
    """

    cdef Py_ssize_t n = data.shape[0]
    cdef Py_ssize_t t = data.shape[1]

    cdef np.ndarray[DTYPE_t, ndim = 2] results = np.zeros((t, n))

    cdef Py_ssize_t i

    for i in prange(n, nogil=True):
        results[i, :] = transpose_vector(data[:, i])

    return results

cdef transpose_vector(np.ndarray[DTYPE_t, ndim=1] vector):
    """
    transposes a np vector
    """
    return vector.transpose()

a = np.random.rand(100, 20)
transpose_example(a)

输出

Converting to Python object not allowed without gil

并行执行此操作的最佳方法是什么？

Answer 1

您可以在没有 GIL 的情况下传递类型化内存视图切片 (cdef transpose_vector(DTYPE_t[:] vector)) - 这是较新类型化内存视图语法优于 np.ndarray 的主要优势之一。

然而，

您不能在内存视图上调用 Numpy 成员函数（如转置），除非您转换回 Numpy 数组 (np.asarray(vector))。这需要 GIL。
调用任何类型的 Python 函数（例如 transpose）都需要 GIL。这可以在 with gil: 块内完成，但是当该块几乎是你的整个循环时就变得毫无意义了。
您没有为 transpose_vector 指定 return 类型，因此它将默认为 object，这需要 GIL。您可以指定一个 Cython return 类型，但我怀疑甚至 returning 内存视图切片可能需要在某处进行一些引用计数。
注意不要让多个线程覆盖您传递的 memoryview 切片中的相同数据。

总而言之：memoryview 切片，但请记住，如果没有 GIL，您可以做的事情非常有限。您当前的示例不可并行化（但这可能主要是因为它是一个玩具示例）。

Answer 2

Q : "What would be the best way to do this in parallel?"
+
""

让我从自由解释亨利福特的格言开始：汽车中缺陷最少的部分就是那个根本不存在的部分 - 它永远不会损坏。

那些知道 numpy 的内部数组对象表示如何工作的人肯定，
需要 _几乎零时间
并且
它需要 _几乎零memory-I/Os _{（至少在过去十年左右是这样）}

为什么？

Numpy 很聪明.
它根本不会为此移动任何数据。它只是以大约 15 [us]

的成本调整索引

答案：

有史以来最好的 np.transpose() 方法根本不需要任何并行化。

任何尝试这样做都会导致性能下降，因为人为地强制执行了许多无用的 memory-I/O-s，而原生 np.transpose() 从来没有这样做过 - 它只是交换索引方案，而不移动大量的单元格数据，全部位于其原始位置（包括保持所有缓存一致性有效 - 因此任何下一次访问都将从缓存中进行 0.5 ~ 5 [ns] ，无需再次支付大量费用 150-350 [ns] 来移动任何单元格数据 from/to 物理 RAM 位置并破坏缓存行）

证明：

______1x THE PROBLEM SIZE 100 x 20 x 8 [B] ( 16 kB RAM).transpose() ~ 17 [us]

>>> r =   100; c =   20; a = np.arange(r*c).reshape( (r,c) );a.itemsize*a.size/1E6
0.016
>>> aClk.start(); _ = a.transpose(); aClk.stop()  ####  16   [kB] RAM-footprint
17                                                ####  17   [us] !!! ZERO-COPY

____100x THE PROBLEM SIZE 1000 x 200 x 8 [B] (1.6 MB RAM).transpose() ~ 17 [us]

>>> r =  1000; c =  200; a = np.arange(r*c).reshape( (r,c) );a.itemsize*a.size/1E6
1.6
>>> aClk.start(); _ = a.transpose(); aClk.stop()  ####   1.6 [MB] RAM-footprint
17                                                ####  17   [us] !!! ZERO-COPY

__10000x THE PROBLEM SIZE 10000 x 2000 x 8 [B] (160 MB RAM).transpose() ~ 16 [us]

>>> r = 10000; c = 2000; a = np.arange(r*c).reshape( (r,c) );a.itemsize*a.size/1E6
160.0
>>> aClk.start(); _ = a.transpose(); aClk.stop()  #### 160.0 [MB] RAM-footprint
16                                                ####  16   [us] !!! ZERO-COPY

像任何 prange-d 或任何其他代码一样，移动或 ALAP 分配和复制那么多 RAM 存储的单元数据，需要很长时间，而不是像智能 numpy 设计那样的智能 16 [ns]。

>>> a.flags
  C_CONTIGUOUS    : True <-------------------------ORIGINAL [indirect] RAM-indexing
  F_CONTIGUOUS    : False
  OWNDATA         : False
  WRITEABLE       : True
  ALIGNED         : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY    : False
>>>
>>> _.flags
  C_CONTIGUOUS    : False
  F_CONTIGUOUS    : True <------------------------TRANSPOSE'd
  OWNDATA         : False
  WRITEABLE       : True
  ALIGNED         : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY    : False

Q.E.D.

结语：

现在的重点来了。希望如此。

在没有 GIL 的情况下在 Cython 中并行化

Parallelize in Cython without GIL

python

parallel-processing

hpc

numpy

cython

为什么？

答案：

证明：

结语：