用于大型矩阵创建/操作的高效Cython

Question

我一直在试图加速创建和操作非常大的数据矩阵（约15,000 x 15,000;类型为double）的代码段。就目前而言，我认为矩阵的大小并不重要，因为即使对于一个10 x 10的小矩阵我也没有看到加速（事实上，对于小矩阵，编译后的cython代码比纯python慢​​，而时间在大型矩阵的cython和python之间几乎相同）。请耐心等待，因为我只编写了python一个星期（最近由Matlab转换而来），而我只是一个不起眼的化学工程师。

代码的目标是将一维数组（长度L）作为输入，对于一个示例：

[ 16.66 16.85 16.93 16.98 17.08 17.03 17.09 16.76 16.67 16.72] 

而产生的矩阵（高度L，宽度L-1）作为输出：

[[ 16.66 16.85 16.93 16.98 17.08 17.03 17.09 16.76 16.67] 
[ 16.85 16.93 16.98 17.08 17.03 17.09 16.76 16.67 16.72] 
[ 16.93 16.98 17.08 17.03 17.09 16.76 16.67 16.72 0. ] 
[ 16.98 17.08 17.03 17.09 16.76 16.67 16.72 0.  0. ] 
[ 17.08 17.03 17.09 16.76 16.67 16.72 0.  0.  0. ] 
[ 17.03 17.09 16.76 16.67 16.72 0.  0.  0.  0. ] 
[ 17.09 16.76 16.67 16.72 0.  0.  0.  0.  0. ] 
[ 16.76 16.67 16.72 0.  0.  0.  0.  0.  0. ] 
[ 16.67 16.72 0.  0.  0.  0.  0.  0.  0. ] 
[ 16.72 0.  0.  0.  0.  0.  0.  0.  0. ]] 

我希望从上面的例子和下面的代码清楚我试图实现什么。该算法需要扩展到非常大的矩阵，它目前没有错误，它只是很慢！

这是我用Cython代码：

from scipy.sparse import spdiags 
import numpy as np 
cimport numpy as np 
cimport cython 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def sfmat(np.ndarray[double, ndim=1] data): 
    cdef int h = data.shape[0] 
    cdef np.ndarray[double, ndim=2] m = np.zeros([h, h-1]) 
    m = np.flipud(spdiags(np.tril(np.tile(data,[h-1,1]).T,0),range(1-h,1), h, h-1).todense()) 
    return m 

我也试图更冗长的代码，这可能是更清楚地写着：

from scipy.sparse import spdiags 
import numpy as np 
cimport numpy as np 
cimport cython 

DTYPE = np.float 
ctypedef np.float_t DTYPE_t 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def sfmat(np.ndarray[DTYPE_t, ndim=1] data): 
    assert data.dtype == DTYPE 
    cdef int h = data.shape[0] 
    cdef np.ndarray[DTYPE_t, ndim=2] m = np.zeros([h, h-1], dtype=DTYPE) 
    cdef np.ndarray[DTYPE_t, ndim=2] s1 = np.zeros([h, h-1], dtype=DTYPE) 
    cdef np.ndarray[DTYPE_t, ndim=2] s2 = np.zeros([h, h-1], dtype=DTYPE) 
    cdef np.ndarray[DTYPE_t, ndim=2] s3 = np.zeros([h, h-1], dtype=DTYPE) 

    s1 = np.tile(data,[h-1,1]).T 
    s2 = np.tril(s1,0) 
    s3 = spdiags(s2,range(1-h,1), h, h-1).todense() 
    m = np.flipud(s3) 
    return m 

与用Cython实行任何帮助将是非常赞赏。如果还有其他方法可以加速这个算法，那也会有所帮助。 谢谢你的帮助！

因为我是新手，这里有更多的细节，这可能会或可能不会阻止我加快速度。 我正在运行64位Windows 7 Pro，并且使用Windows SDK C/C++编译器成功编译了cython代码。 （我按照github here上的指示成功）。简单的“hello world”cython示例编译得很好，在64位模式下运行良好，上面的代码也编译并运行时没有错误。对于整个15,000 x 15,000矩阵的处理，需要64位体系结构，或者至少我认为是这样，因为在编译32位后运行代码导致内存错误。对于这个问题，请假设将矩阵分成更小的块是不可能的。 请让我知道是否有任何其他信息需要回答这个问题。

干杯，scientistR

UPDATE

我认为避免for循环将是最好的方法，但是，spdiags是主要瓶颈。于是，一种新的算法效果更好（4倍提高我的电脑上）：

import numpy as np 
cimport numpy as np 
cimport cython 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def sfmat(np.ndarray[double, ndim=1] data): 
    cdef int i 
    cdef np.ndarray[double, ndim=2] m = np.zeros([data.shape[0], data.shape[0]-1]) 
    for i in range(data.shape[0]-1): 
     m[:,i] = np.roll(data,-i); 
    return m 

但用Cython不会提供超过纯Python任何改善。请帮忙。正如评论家指出的那样，除了更优化算法之外，可能还没有办法改进这一点，但我充满希望。谢谢！另外，还有更快的算法，cython或python？

Answer 1

这可能是一个老问题，但不应该留下任何问题没有回答:)。我能够通过使用简单的for循环（Cython中实际上是快速的）来加速您的Cython代码大约8倍，对于数组大小为7000 ..注意你的实现使用np.roll不会产生你想要的数组（！），但我使用该函数来比较时序。

编辑的代码中使用类型Memoryviews和np.empty代替np.zeros

def sfmat(double[:] data): 
    cdef int n = data.shape[0] 
    cdef np.ndarray[double, ndim=2] out = np.empty((n, n-1)) 
    cdef double [:, :] out_v = out # "typed memoryview" 

    cdef int i, j 
    for i in range(n-1): 
     out_v[0, i] = data[i] 

    for i in range(1, n): 
     for j in range(n-i): 
      out_v[i, j] = data[i+j] 
     for j in range(n-i, n-1): 
      out_v[i, j] = 0. 
    return out 

不幸的是，用Cython努力是只有约1.2倍比在常规的Python会话中运行下面的代码更快：

def sfmat(data): 
    n = len(data) 
    out = np.empty((n, n-1)) 
    out[0, :] = data[:n-1] 
    for i in xrange(1, n): 
     out[i, :n-i] = data[i:] 
     out[i, n-i:] = 0 
    return out 

但是，正如评论中已经讨论的那样，以这种方式炸毁原来相当小的矩阵可能不是解决实际问题的最有效方法，总体上问题无论如何。如果你最初想做的事情是避免使用for循环，那么在Cython中根本就没有必要这么做！

Answer 2

我并不是说听起来很天真，但我们都知道C，C++和Python是“行 - 主”语言，对吗？ Matlab（和Fortran）是“列专业”。我相信你已经尝试了扭转i和j，但只是想提及它，以防万一没有人想到这一点。