这个 NumPy 操作的执行速度能否与其 Cython 等效操作一样快，或者更快？

Question

TLDR; 我正在执行数组运算（没有数学运算），我发现 Cython 的速度要快得多。有什么办法可以在 NumPy 中加快速度吗？还是 Cython？

上下文

我正在编写一个函数，用于从 index 两个方向（其顶角沿对角线）获取 NxN 数组的子集并将其移动一个位置向上沿对角线。其次，我需要将第一行从 index 向左移动一位。最后，我需要在操作后将数组中的最后一列设置为零。

该数组是一个严格的上三角矩阵，这意味着从对角线向下的所有内容都设置为0。这是我尝试以一种优雅的方式来存储对之间的历史碰撞数据对象（其索引由矩阵中的索引表示）。这类似于制作一个大小为 n!/(2(n-2)!) 的嵌套列表，它表示长度为 n 的索引列表的有序对。在这个算法中，我希望从碰撞配对矩阵中“移除”一个物体。

我在这个实现中发现的优点是，从矩阵中“删除碰撞对”比从嵌套列表中删除对并将索引成对移动到“要删除的索引”点的计算强度要小得多。

整个项目的中心是将 3D 模型自动“打包”到用于粉末床熔合增材制造的构建体积中。该算法使用模拟退火，因此 p运行e 碰撞集、存储历史信息、add/remove 几何的能力是最重要的，需要很好地优化。

例子

假设我们的数组采用这种形式（不代表实际数据）。

arr = 
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 0. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 0. 0. 3. 4. 5. 6. 7. 8. 9.]
 [0. 0. 0. 0. 4. 5. 6. 7. 8. 9.]
 [0. 0. 0. 0. 0. 5. 6. 7. 8. 9.]
 [0. 0. 0. 0. 0. 0. 6. 7. 8. 9.]
 [0. 0. 0. 0. 0. 0. 0. 7. 8. 9.]
 [0. 0. 0. 0. 0. 0. 0. 0. 8. 9.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 9.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

然后使用 index = 3 我们应该获取子集 index+1:n, index+1:n 中的所有内容并将其设置为等于 index:n-1, index:n-1。然后将顶行向左移动；再次从 index 开始。然后将最后一列设置为 0.

fun(3, arr)
[[0. 1. 2. 4. 5. 6. 7. 8. 9. 0.]
 [0. 0. 2. 3. 4. 5. 6. 7. 8. 0.]
 [0. 0. 0. 3. 4. 5. 6. 7. 8. 0.]
 [0. 0. 0. 0. 5. 6. 7. 8. 9. 0.]
 [0. 0. 0. 0. 0. 6. 7. 8. 9. 0.]
 [0. 0. 0. 0. 0. 0. 7. 8. 9. 0.]
 [0. 0. 0. 0. 0. 0. 0. 8. 9. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 9. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

实现一：纯NumPy

再次假设 arr 是一个 NxN 矩阵。

def fun(index, n, arr):
     arr[index:-1, index:-1] = arr[index + 1:, index + 1:]
     arr[0, index:-1] = arr[0, index + 1:]
     arr[:, n-1:] = 0
     return arr

实现二：Cython

请耐心等待，因为这是我第一次实现 Cython。

@cython.boundscheck(False)
def remove_from_collision_array(int index, int n, double[:,:] arr):

    cdef int i, j, x_shape, y_shape
    x_shape = arr.shape[0]

    for i in range(index, x_shape):
        for j in range(index, x_shape):
            if j <= i:
                # We are below the diagonal, do nothing
                continue
            elif i >= n-1 or j >= n-1:
                arr[i, j] = 0
            else:
                arr[i, j] = arr[i+1, j+1]

    arr[0, index:-1] = arr[0, index+1:]
    arr[:, n-1:] = 0

    return np.asarray(arr)

讨论

在有人不高兴之前，是的，我不知道我在 Cython 中做什么。我禁用了 bounds_checking 因为它确实加快了速度。我正在使用我的 elif 语句之一在循环中执行边界检查。

我最初认为在循环中执行此操作不会比 NumPy 更快。我预先分配了一个大小为 5000x5000 的 NumPy 数组，以避免需要动态追加等。我什至使用与 Numpy 相同的 3 行代码测试了 Cython 实现，但它的性能也很差。

您可以看到使用 index=0 将需要最多的计算。所以我用它作为基准。在循环测试时，我发现 Cython 实现比 Numpy 版本快 50% 以上。也许这是因为我没有充分使用 NumPy 提供的工具？

我绝不是计算机科学家，我也不知道这是不是最好的路线。我是一名设计系统原型的设计师。如果有人知道如何使尖叫声更快，请告诉我！

更新答案

感谢杰罗姆今天教了我一些东西！这将有助于以闪电般的速度制作此包运行。我已将他的见解添加到我的代码中，导致性能大幅提升，原因有二：

我通过在对角线上方开始 j 循环，将循环迭代次数减少了 n*(n-1)/2。
我删除了所有条件语句。

这是更新后的 Cython：

@cython.boundscheck(False)
@cython.wraparound(False)
def remove_from_collision_arrayV2(int index, int n, double[:,:] arr):

    cdef int i, j

    # Shift the diagonal matrix
    for i in range(index, n-1):
        for j in range(i, n-1):
                arr[i, j] = arr[i+1, j+1]

    # Shift the rop row
    for j in range(index, n-1):
        arr[0, j] = arr[0, j+1]

    # Set Column column n-1 to zero
    for i in range(n):
        arr[i, n-1] = 0

    return np.asarray(arr)

用于基准测试目的。在 500x500 矩阵上使用 index=0 执行此迭代 500 次：

原始 NumPy 代码：52.8s

原始 Cython 代码：16.47s - 3.2 倍加速

更新的 Cython 代码：0.014s - 3550x 加速

Answer 1

表达式 arr[index:-1, index:-1] = arr[index + 1:, index + 1:] 在 Numpy 和 Cython 中都很慢而 Cython 代码快得多的原因有点违反直觉：此表达式在 Numpy 和 Cython.

中均未有效实现

事实上，Numpy 将右侧 (arr[index + 1:, index + 1:]) 复制到动态分配的 临时数组 中。然后将临时数组复制到左侧 (arr[index:-1, index:-1])。这意味着执行了两次内存复制，而只能使用一次。更糟糕的是：复制的内存非常大，无法放入缓存，从而导致更大的开销（在某些处理器上，如主流 x86/x86-64 处理器，write-back policy 会导致额外的缓慢读取）。此外，新的临时数组会导致许多页面错误，从而减慢复制速度。

Numpy 这样做是因为左侧和右侧可能重叠（这里就是这种情况），因此复制内存字节的顺序很重要。 Numpy 使用缓慢保守的方法而不是优化的实现。这是一个 遗漏的优化 。 Cython 做的完全一样。

您的 Cython 代码不会受到所有这些开销的影响：它直接相对高效地就地复制数组。读取的值保存在缓存中，然后立即写入，这样回写策略就不是问题了。此外，没有临时数组或页面错误。最后，您的 Cython 代码不会复制三角矩阵的下半部分，因此与前面提到的表达式相比，要复制的字节更少。

减少 Numpy 表达式开销的一种方法是逐块复制矩阵并为此分配一个小的临时缓冲区（通常是矩阵的几行）。然而，这远非易事，因为 CPython 循环通常非常慢，并且块大小应适合缓存，因此该方法很有用...

进一步优化： 条件语句很慢。您可以通过在 i+1 处开始基于 j 的循环并在 n-1 处结束来删除它们。然后另一个基于 j 的循环可以填充大于 n-1 的值。出于同样的原因，基于 i 的循环应该在 n-1 处结束，然后另一个循环可以填充数组的剩余部分。一个好的编译器应该使用更快的 SIMD 指令.

这个 NumPy 操作的执行速度能否与其 Cython 等效操作一样快，或者更快？

Can this NumPy operation perform as fast, or faster, than its Cython Equivalent?

python

arrays

numpy

cython