首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >用于大型矩阵创建/操作的高效Cython

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

Stack Overflow用户
提问于 2013-10-18 09:12:37
回答 2查看 1K关注 0票数 0

我一直试图加快一段代码,以创建和操作一个非常大的数据矩阵(大约。15,000 x 15,000;双型)。现在,我不认为矩阵的大小是那么重要,因为即使对于一个小的10×10矩阵,我也看不到加速比(实际上,编译的cython代码比小矩阵的纯python要慢,而对于大型矩阵而言,cython和python的时间几乎相同)。请耐心地对待我,因为我只编写了一个星期的python代码(刚从Matlab转换而来),而且我只是一个卑微的化学工程师。

代码的目标是以一个一维数组(长度L)作为输入,例如:

代码语言:javascript
运行
复制
[ 16.66  16.85  16.93  16.98  17.08  17.03  17.09  16.76  16.67  16.72]

并产生一个矩阵(高度L,宽度L1)作为输出:

代码语言:javascript
运行
复制
[[ 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代码:

代码语言:javascript
运行
复制
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

我还尝试了更详细的代码,这些代码可能更清晰:

代码语言:javascript
运行
复制
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 /C++编译器成功地编译了cython代码。(我成功地遵循了github 这里上的指示)。简单的"hello“cython示例在64位模式下编译并运行良好,上面的代码也编译和运行,没有错误。对于整个15,000 x 15,000矩阵的操作,需要64位架构,至少我相信是这样的,因为在编译32位后运行代码会导致内存错误。对于这个问题,请假设将矩阵分解成较小的块是不可能的。请让我知道是否有任何其他资料需要回答这个问题。

干杯,scientistR

更新

我认为避免for循环将是最好的方法,但是spdiags是主要的瓶颈。因此,一种新的算法工作得更好(我的计算机上改进了4倍):

代码语言:javascript
运行
复制
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?

EN

回答 2

Stack Overflow用户

回答已采纳

发布于 2013-11-01 11:41:33

这可能是一个老问题,但是:)不应该留下任何问题。通过使用简单的for -循环(实际上在Cython中速度很快),数组大小为7000,我能够将您的Cython代码速度提高大约8倍。注意,使用np.roll的实现不会产生您想要的数组(!),但是我使用这个函数来比较时间。

编辑后的代码使用类型化内存视图和np.empty代替np.zeros

代码语言:javascript
运行
复制
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的工作速度仅比在常规Python会话中运行以下代码快1.2倍:

代码语言:javascript
运行
复制
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中根本不需要这样做!

票数 0
EN

Stack Overflow用户

发布于 2013-11-08 08:06:50

我不想听起来很天真,但我们都知道C、C++和Python是“主要的”语言,对吗?Matlab (和Fortran)是“列专业”。我相信您已经尝试过逆转ij,但只想提到它,以防有人想到尝试。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/19445820

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档