C阵列与NumPy阵列

Question

在性能（代数操作、查找、缓存等）方面，C数组（可以公开为C数组或cython.view.array[Cython数组]或上述两个数组的memoryview）与NumPy数组（Cython中不应该有Python开销）之间是否存在差异

编辑：

我应该提到，在NumPy数组中使用Cython进行静态类型化，并且dtype是NumPy编译时数据类型（例如cdef np.int_t或cdef np.float32_t），在C情况下的类型是C等价的（cdef int_t和cdef float）

编辑2：

下面是Cython Memoryview documentation中的示例，以进一步说明我的问题：

from cython.view cimport array as cvarray
import numpy as np

# Memoryview on a NumPy array
narr = np.arange(27, dtype=np.dtype("i")).reshape((3, 3, 3))
cdef int [:, :, :] narr_view = narr

# Memoryview on a C array
cdef int carr[3][3][3]
cdef int [:, :, :] carr_view = carr

# Memoryview on a Cython array
cyarr = cvarray(shape=(3, 3, 3), itemsize=sizeof(int), format="i")
cdef int [:, :, :] cyarr_view = cyarr

坚持aC array与aCython array与aNumPy array有什么区别吗？

Answer 1

不要使用cython.view.array，使用cpython.array.array。

有关详细信息，请参见this answer of mine，尽管这只涉及速度。建议将cython.view.array视为“演示”材料，并将cpython.array.array视为实际的可靠实现。这些阵列非常轻量级，而且在用作暂存空间时效果更好。

此外，如果您曾经受到malloc的诱惑，那么对它们的原始访问不会变慢，实例化只需要两倍的时间。

---

关于伊恩的

If you need both the array slicing and the NumPy functionality for a given array, you can make a memory view that points to the same memory as the NumPy array.

值得注意的是，memoryView有一个“base”属性，许多Numpy函数也可以使用memoryView，因此这些不必是分离变量。

Answer 2

我对这方面的知识还不完善，但这可能会有所帮助。 我运行了一些非正式的基准测试来展示每种数组类型的优点，并对我的发现很感兴趣。

尽管这些数组类型在许多方面不同，但如果您使用大型数组进行大量计算，那么您应该能够从其中任何一个数组中获得类似的性能，因为逐项访问在总体上应该大致相同。

NumPy数组是使用Python的C API实现的Python对象。 NumPy数组确实在C级别提供了API，但是它们不能独立于Python解释器创建。 它们特别有用，因为NumPy和SciPy中提供了所有不同的数组操作例程。

Cython内存视图也是Python对象，但它是作为Cython扩展类型生成的。 它似乎不是为纯Python设计的，因为它不是Cython的一部分，不能直接从Python导入，但是可以从Cython函数返回Python的视图。 您可以在https://github.com/cython/cython/blob/master/Cython/Utility/MemoryView.pyx查看实现

C数组是C语言中的本机类型。 它像指针一样被索引，但是数组和指针是不同的。 在http://c-faq.com/aryptr/index.html有一些很好的讨论 它们可以在堆栈上分配，并且更容易让C编译器进行优化，但是它们在Cython之外更难访问。 我知道你可以从其他程序动态分配的内存中创建一个NumPy数组，但这样做似乎困难得多。 特拉维斯·奥列芬特在http://blog.enthought.com/python/numpy-arrays-with-pre-allocated-memory/上发布了一个这样的例子 如果在程序中使用C数组或指针作为临时存储，它们应该非常适合您。 它们对于切片或任何其他类型的矢量计算都不太方便，因为您必须自己使用显式循环来完成所有工作，但是它们应该更快地分配和释放，并且应该为速度提供一个良好的基线。

Cython还提供了一个数组类。 看起来它是为内部使用而设计的。 复制memoryview时创建实例。 见http://docs.cython.org/src/userguide/memoryviews.html#view-cython-arrays

在Cython中，您还可以分配内存并为指针编制索引，使分配的内存有点像数组。 见http://docs.cython.org/src/tutorial/memory_allocation.html

下面是一些基准测试，它们在索引大型数组时显示了一些类似的性能。 这是Cython文件。

from numpy cimport ndarray as ar, uint64_t
cimport cython
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t i, j
    for i in range(n):
        for j in range(size):
            A[j] = n

def carr_time(uint64_t n=1000000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t AC[10000]
        uint64_t a
        int i, j
    for i in range(n):
        for j in range(10000):
            AC[j] = n

@cython.boundscheck(False)
@cython.wraparound(False)
def ptr_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t* AP = &A[0]
        uint64_t a
        int i, j
    for i in range(n):
        for j in range(size):
            AP[j] = n

@cython.boundscheck(False)
@cython.wraparound(False)
def view_time(uint64_t n=1000000, uint64_t size=10000):
    cdef:
        ar[uint64_t] A = np.empty(n, dtype=np.uint64)
        uint64_t[:] AV = A
        uint64_t i, j
    for i in range(n):
        for j in range(size):
            AV[j] = n

使用IPython计时

%timeit -n 10 ndarr_time()
%timeit -n 10 carr_time()
%timeit -n 10 ptr_time()
%timeit -n 10 view_time()

10 loops, best of 3: 6.33 s per loop
10 loops, best of 3: 3.12 s per loop
10 loops, best of 3: 6.26 s per loop
10 loops, best of 3: 3.74 s per loop

这些结果让我觉得有点奇怪，因为根据Efficiency: arrays vs pointers，数组不太可能比指针快得多。 似乎某种编译器优化正在使纯C数组和类型化内存视图更快。 我试着关闭C编译器上的所有优化标志，并得到了时间安排

1 loops, best of 3: 25.1 s per loop
1 loops, best of 3: 25.5 s per loop
1 loops, best of 3: 32 s per loop
1 loops, best of 3: 28.4 s per loop

在我看来，逐项访问在总体上几乎是一样的，只是C数组和Cython内存视图似乎更容易让编译器进行优化。

在我前段时间发现的这两篇博文中，可以看到更多关于这一点的评论： http://jakevdp.github.io/blog/2012/08/08/memoryview-benchmarks/http://jakevdp.github.io/blog/2012/08/16/memoryview-benchmarks-2/

在第二篇博文中，他评论了如果内存视图切片是内联的，那么它们可以提供与指针算法类似的速度。 我在自己的一些测试中注意到，显式地内联使用内存视图切片的函数并不总是必要的。 作为这个例子，我将计算数组中两行的每一个组合的内积。

from numpy cimport ndarray as ar
cimport cython
from numpy import empty

# An inlined dot product
@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline double dot_product(double[:] a, double[:] b, int size):
    cdef int i
    cdef double tot = 0.
    for i in range(size):
        tot += a[i] * b[i]
    return tot

# non-inlined dot-product
@cython.boundscheck(False)
@cython.wraparound(False)
cdef double dot_product_no_inline(double[:] a, double[:] b, int size):
    cdef int i
    cdef double tot = 0.
    for i in range(size):
        tot += a[i] * b[i]
    return tot

# function calling inlined dot product
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_slicing(ar[double,ndim=2] A):
    cdef:
        double[:,:] Aview = A
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = dot_product(Aview[i], Aview[j], A.shape[1])
    return res

# function calling non-inlined version
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_slicing_no_inline(ar[double,ndim=2] A):
    cdef:
        double[:,:] Aview = A
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = dot_product_no_inline(Aview[i], Aview[j], A.shape[1])
    return res

# inlined dot product using numpy arrays
@cython.boundscheck(False)
@cython.boundscheck(False)
cdef inline double ndarr_dot_product(ar[double] a, ar[double] b):
    cdef int i
    cdef double tot = 0.
    for i in range(a.size):
        tot += a[i] * b[i]
    return tot

# non-inlined dot product using numpy arrays
@cython.boundscheck(False)
@cython.boundscheck(False)
cdef double ndarr_dot_product_no_inline(ar[double] a, ar[double] b):
    cdef int i
    cdef double tot = 0.
    for i in range(a.size):
        tot += a[i] * b[i]
    return tot

# function calling inlined numpy array dot product
@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_dot_rows_slicing(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = ndarr_dot_product(A[i], A[j])
    return res

# function calling nun-inlined version for numpy arrays
@cython.boundscheck(False)
@cython.wraparound(False)
def ndarr_dot_rows_slicing_no_inline(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            res[i,j] = ndarr_dot_product(A[i], A[j])
    return res

# Version with explicit looping and item-by-item access.
@cython.boundscheck(False)
@cython.wraparound(False)
def dot_rows_loops(ar[double,ndim=2] A):
    cdef:
        ar[double,ndim=2] res = empty((A.shape[0], A.shape[0]))
        int i, j, k
        double tot
    for i in range(A.shape[0]):
        for j in range(A.shape[0]):
            tot = 0.
            for k in range(A.shape[1]):
                tot += A[i,k] * A[j,k]
            res[i,j] = tot
    return res

我们看到的时间

A = rand(1000, 1000)
%timeit dot_rows_slicing(A)
%timeit dot_rows_slicing_no_inline(A)
%timeit ndarr_dot_rows_slicing(A)
%timeit ndarr_dot_rows_slicing_no_inline(A)
%timeit dot_rows_loops(A)

1 loops, best of 3: 1.02 s per loop
1 loops, best of 3: 1.02 s per loop
1 loops, best of 3: 3.65 s per loop
1 loops, best of 3: 3.66 s per loop
1 loops, best of 3: 1.04 s per loop

使用显式内联的结果和没有内联的结果一样快。 在这两种情况下，类型化内存视图都可以与未经切片编写的函数版本相比。

在博客文章中，他必须编写一个特定的示例来强制编译器不内联函数。 一个好的C编译器（我使用的是MinGW）似乎能够处理这些优化，而不必被告知要内联某些函数。 内存视图可以是faster用于在Cython模块中的函数之间传递数组切片，即使没有显式内联。

然而，在这种特殊的情况下，即使将循环推到C，也不会真正达到通过适当使用矩阵乘法所能达到的速度。 布拉斯仍然是做这种事情的最好方法。

%timeit A.dot(A.T)
10 loops, best of 3: 25.7 ms per loop

也有从NumPy数组到memoryView的自动转换，如

cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def cysum(double[:] A):
    cdef tot = 0.
    cdef int i
    for i in range(A.size):
        tot += A[i]
    return tot

一个问题是，如果希望函数返回NumPy数组，则必须使用np.asarray将内存视图对象再次转换为NumPy数组。 这是一个相对便宜的操作，因为内存视图符合http://www.python.org/dev/peps/pep-3118/

结论

对于Cython模块内部使用的NumPy数组，类型化内存视图似乎是一个可行的替代方案。 对于内存视图，数组切片速度会更快，但对于内存视图，没有为NumPy数组编写的那么多函数和方法。 如果您不需要调用一堆NumPy数组方法，并且想要简单的数组切片，那么可以使用内存视图代替NumPy数组。 如果需要对给定数组同时使用数组切片和的NumPy功能，则可以创建指向与NumPy数组相同内存的内存视图。 然后可以使用视图在函数之间传递切片，并使用数组调用NumPy函数。 这种方法仍然有一定的局限性，但是如果您使用一个数组来完成大部分处理，那么它会很好地工作。

C数组和/或动态分配的内存块对于中间计算可能很有用，但是它们不容易传递回Python以便在那里使用。 在我看来，动态分配多维C数组也比较麻烦。 我知道的最好的方法是分配一大块内存，然后使用整数算法将其作为多维数组进行索引。 如果您希望轻松地动态分配数组，这可能是一个问题。 另一方面，对于C数组，分配时间可能要快一些。 其他数组类型的设计速度和使用方便得多，因此我建议您使用它们，除非有令人信服的理由要使用它们。

更新：正如@Veedrac在回答中提到的，您仍然可以将Cython内存视图传递给大多数NumPy函数。 当您这样做时，NumPy通常需要创建一个新的NumPy数组对象来处理内存视图，因此这会稍微慢一些。 对于大型阵列，这种影响可以忽略不计。 无论数组大小如何，对内存视图的np.asarray调用都相对较快。 然而，为了证明这一效果，这里有另一个基准：

Cython文件：

def npy_call_on_view(npy_func, double[:] A, int n):
    cdef int i
    for i in range(n):
        npy_func(A)

def npy_call_on_arr(npy_func, ar[double] A, int n):
    cdef int i
    for i in range(n):
        npy_func(A)

在伊普顿：

from numpy.random import rand
A = rand(1)
%timeit npy_call_on_view(np.amin, A, 10000)
%timeit npy_call_on_arr(np.amin, A, 10000)

输出：

10 loops, best of 3: 282 ms per loop
10 loops, best of 3: 35.9 ms per loop

我试着选择一个能很好地展示这种效果的例子。 除非涉及到对相对较小数组的许多NumPy函数调用，否则这不应该改变整个时间。 请记住，无论我们以何种方式调用NumPy，Python函数调用仍然会发生。

这只适用于NumPy中的函数。 大多数数组方法都不能用于memoryView（一些属性仍然是，比如size和shape和T）。 例如，使用NumPy数组的A.dot(A.T)将变成np.dot(A, A.T)。