使用Cython优化NumPy

Question

我目前正在尝试优化我用纯Python编写的代码。这段代码非常依赖于NumPy，因为我正在使用NumPy数组。下面是我转换为Cython的最简单的类之一，它只是对两个Numpy数组进行乘法运算。如下：

bendingForces = self.matrixPrefactor * membraneHeight

我的问题是，我是否可以优化这个过程，以及如何优化。当我查看“cython -a”生成的C代码时，发现有很多NumPy调用，这看起来不太高效。

import numpy as np
cimport numpy as np
ctypedef np.float64_t dtype_t
ctypedef np.complex128_t cplxtype_t
ctypedef Py_ssize_t index_t

    cdef class bendingForcesClass( object ):
        cdef dtype_t bendingRigidity
        cdef np.ndarray matrixPrefactor
        cdef np.ndarray bendingForces

        def __init__( self, dtype_t bendingRigidity, np.ndarray[dtype_t, ndim=2] waveNumbersNorm ):
            self.bendingRigidity = bendingRigidity
            self.matrixPrefactor = -self.bendingRigidity * waveNumbersNorm**2

        cpdef np.ndarray calculate( self, np.ndarray membraneHeight ) :
            cdef np.ndarray bendingForces
            bendingForces = self.matrixPrefactor * membraneHeight
            return bendingForces

我想使用两个for循环来迭代数组的条目。也许我可以使用编译器通过SIMD操作来优化它？！我尝试过，可以编译，但结果奇怪并且需要很长时间。这是替换函数的代码：

cpdef np.ndarray calculate( self, np.ndarray membraneHeight ) :

    cdef index_t index1, index2 # corresponds to: cdef Py_ssize_t index1, index2
    for index1 in range( self.matrixSize ):
        for index2 in range( self.matrixSize ):
            self.bendingForces[ index1, index2 ] = self.matrixPrefactor.data[ index1, index2 ] * membraneHeight.data[ index1, index2 ]
    return self.bendingForces

然而，正如我所说的那样，这段代码非常慢，并且不能按预期运行。那么我做错了什么？最佳的优化方式是什么，如何消除NumPy调用操作？

- packoman

---

你应该尽可能避免在Cython中重写所有代码。我建议先检查哪些部分是瓶颈，然后只优化那些部分。不确定相对于NumPy的逐元素数组乘法，你是否会获得巨大的性能提升。 - JoshAdel

你能给我们提供一个这个函数的输入大小和输出的例子吗？我知道如何优化它，但需要知道您如何使用它以及它的速度有多快。 - fabrizioM

如果你总是渴望更好的性能，考虑离开Python。Cython并不是魔法。 - BatchyX

你认为自己能比numpy做得更好，是什么让你有这样的想法？ - David Heffernan

@David Heffernan：在某些情况下，使用Cython击败numpy并不难，例如，请参见我的答案https://dev59.com/Hm445IYBdhLWcg3wQX-G - jfs

Answer 1

对于简单的矩阵乘法，NumPy代码已经在本地执行循环和乘法，所以很难在Cython中超越它。 Cython非常适用于将Python中的循环替换为Cython中的循环的情况。您的代码比NumPy慢的原因之一是因为每次在数组中进行索引查找时，

self.bendingForces[ index1, index2 ] = self.matrixPrefactor.data[ index1, index2 ] * membraneHeight.data[ index1, index2 ]

它会执行更多的计算，例如边界检查（索引是否有效）。如果将索引强制转换为无符号整数，您可以在函数前使用修饰符@cython.boundscheck(False)。

有关加速Cython代码的更多详细信息，请参见此教程。

- highBandWidth

Answer 2

你可以通过使用（某些方法）来加速这个过程。

for index1 from 0 <= index1 < max1:

不要使用我不确定类型的范围。

你有检查过这个和那个吗？

- Joost Rekveld

---

这是旧的Cython语法。当在循环中使用时，range和xrange都会被优化掉。 - Fred Foo

是的，这现在是旧的Cython语法，但是在我发布这篇文章的时候它还不是。尽管如此，还是很值得指出的。 - Joost Rekveld