我正在寻找一种简洁有效的方法来获取多维数组切片,对这些切片执行标量和矩阵运算,然后最终将结果数组保存为另一个数组中的切片.

您可以使用以下语法在fortran中做得很好：

real*8, dimension(4,4,4,4)               ::  matrix_a
real*8, dimension(4,4)               ::  matrix_b
...
matrix_a(:, 2, :, 4) =  matrix_a(:, 2, :, 4) + (2 * matrix_b(:, :))

我试图在Cython中找到这样做的方法.这是我能想到的最好的：

cdef double matrix_a[4][4][4][4]
cdef double matrix_b[4][4]
...
cdef int i0, i1
for i0 in range(4):
  for i1 in range(4):
    matrix_a[i0][1][i1][3] += (2 * matrix_b[i0][i1])

显然,您可以使用Numpy数组以非常简洁的方式执行此操作,但这似乎比上面的代码要慢几个数量级.有什么方法可以让我在这两个世界中都做到最好吗？可能某种方式我可以以更快的方式调用Numpy？或者也许某种C/C++ API可以从Cython中使用.谢谢

解决方法:

好问题.最简洁的答案是不.

答案很长……这取决于你想要什么以及你愿意付出多少努力.没有一个答案,因为问题非常广泛,不是因为操作不可行.

在撰写本文时,没有默认方法可以对cython的内存视图执行fortran-speed算法.
在Cython邮件列表上有关于此的多个讨论.
这个recent thread很好地总结了这种情况.

现在,有很多方法可以使用NumPy数组进行快速算术运算.没有一个具有其他功能,但有些可能适合您的需求.

第一：

> NumPy的矢量化写得非常好.在类型推断中预先存在一些开销,但是,对于紧密循环内的小数组上的操作以外的任何事情,大部分计算时间都在C中非常优化的循环中运行.如果您正在做一个简单的 – 像这样放置操作,即使在Cython中,最好的方法就是做一些像np.add(a,b,out = c)这样的东西.如果您有更多特定需求,np.einsum,scipy的blas / lapack包装器和numpy的数组缩减方法提供了更大的灵活性.

在大多数情况下,如果a和b是numpy数组,就像

a[:,2,:,4] += 2 * b

应该足够了.

即使你在Cython中运行内存视图而不是numpy数组,你也可以这样做：

import numpy as np
np.add(a, b, out=c)

调用像这样的Python函数有一个固定的成本,但是如果数组很大,那就不重要了.

在你的情况下,numpy仍然会慢一两个原因.
首先,启动循环操作的固定成本很高,因为numpy会在runtine处理很多类型和形状元数据.
对于较大的数组,这种差异消失了,所以,除非你在一个紧密的循环中,否则无关紧要.
其次,numpy和fortran的默认内存布局是不同的.
这种差异可能只会出现在具有较大阵列的操作中.
如果以fortran内存顺序初始化numpy数组并且数组很大,那么numpy应该与fortran相提并论.

现在,对于那些真正处于numpy不够快或者对阵列没有足够控制的情况下的人来说,还有很多其他的选择.
以下是一些更有希望的：

>在Cython中用cdef函数内的循环编写算术运算.让它接受内存视图作为参数并传入预先分配的输出数组.
签名看起来像cdef inline void add(double [：,] a,double [：,：] b,double [：,]] out)nogil：
Memoryview可以在没有Python调用开销的情况下进行切片,因此您可以将切片传递给您定义的函数,而不会产生Python调用的开销.
这种方法不是最快的,但老实说,它通常都足够好.
它还避免了必须手动管理阵列的内存布局所带来的大部分痛苦.
这比使用原始C数组做同样的事情更容易,因为内存视图会动态管理它们的内存.
>存在几个numpy-aware auto-vecotrizers.
考虑尝试numba,parakeet或pythran.
这些包主要集中在编译特定于类型的函数版本,这些函数在numpy函数上运行,然后在适用的情况下调用这些例程.
我和numba和长尾小鹦鹉的结果好坏参半,并且没有和pythran玩过多.
据我所知,numba和长尾小鹦鹉目前没有进行任何形式的静态表达分析,尽管它不一定超出其范围.
Pythran看起来很有希望.
它的主要卖点是对数组表达式进行静态分析(如Fortran).
这些库易于使用,因此可能值得尝试它们以确定它们是否有用.
他们特别擅长优化循环.
>您还可以使用一些Python级别的表达式优化器.
见numexpr和Theano.
粗略地说,这些包专注于将算术表达式编译为机器代码,而不是整个函数.
据我所知,他们没有针对切片进行优化,但是你可以将它们传递给numpy数组.
这些最适合优化大型数组上的操作,但它们通常可以为表达式提供比使用numpy数组的直接算法更快的代码.
>您可以在fortran中实现该操作并使用Cython将其包装.在fortran90 website有关于如何使用fortran执行此操作的说明.我在this answer写了一个如何执行此操作的示例.
你必须自己管理内存布局,但这不是太难,而且会非常快.
>您还可以使用各种C库中的任何一个进行数组操作.
不幸的是,他们中的许多人(eigen,armadillo,blaze-lib)目前都没有对高维数组的内置支持,但你所拥有的特定表达式实际上只是使用跨步的2D数组进行算术运算,因此它可以正常工作.
如果你对这个约束没有问题,那么最初的努力就是为犰狳制作Cython绑定(cy-arma,尽管它们仍然不发达.
一个支持n维数组操作的c库是libdynd.
它还有一套更好的python bindings,用Cython编写.
在C中实现操作并通过Cython将它们包装为Python可能是最简单的,但是您可以从包装器中获取C pxd声明并尝试在Cython级别编写函数.
>有一个名为ceygen的库,它使用eigen作为后端,将内存视图上的算术运算写为函数调用.
我在编译时遇到了麻烦,但从概念上讲,这可能就是你要找的东西.
我不确定它是否能很好地处理跨步信息.
>你可以在C中自己实现这些操作.唯一一次,当你想将sse内在函数应用于数组时(或者某种优于纯循环的优化).写这样的东西非常耗时,但你可以做到.如果你只是循环,你最好在Cython中编写循环并在内存视图上操作.
>您可以使用NumPy C API for ufuncs或gufuncs.这允许您进行标量(或数组)操作并生成将其应用于更通用数组的函数.这些文档现在还不是很完整.
>这对于Cython中的基本算术运算不起作用,但是很快将添加到scipy的一个新功能是BLAS和LAPACK的Cython级包装器.这将使您可以进行各种线性代数运算,而且每个函数调用的成本都很低.目前,这是相关的pull request.

其中一些将要求您自己管理内存布局,但这并不是非常困难.

如果我在你的鞋子,我会尝试选项1,然后选项4.选项1将相当快,但如果这还不够,选项4可能会更快,但现在你将不得不担心内存你自己的阵列布局.