0

我一直在尝试加速创建和操作非常大的数据矩阵的一段代码（大约 15,000 x 15,000；双精度类型）。目前，我认为矩阵的大小并不那么重要，因为即使对于小的 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 上的指示）。简单的“hello world”cython 示例在 64 位模式下编译和运行良好，上面的代码也编译和运行没有错误。对于整个 15,000 x 15,000 矩阵的操作，需要 64 位架构，或者至少我相信是这样，因为在编译为 32 位后运行代码会导致内存错误。对于这个问题，请假设将矩阵分解成更小的块是不可能的。请让我知道是否需要任何其他信息来回答此问题。

干杯，科学家R

更新

我认为避免 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？

4

2 回答 2

0

我并不是说听起来天真，但我们都知道 C、C++ 和 Python 是“行主要”语言，对吧？Matlab（和 Fortran）是“列专业的”。我确定您已经尝试过反转iand j，但只是想提一下，以防万一没有人想过尝试。

于 2013-11-08T08:06:50.300 回答
0

这可能是一个老问题，但任何问题都不应悬而未决:)。通过使用简单的 for 循环（在 Cython 中实际上速度很快），我能够将您的 Cython 代码加速约 8 倍，数组大小为 7000。请注意，您的实现使用np.roll不会产生您想要的数组（！），但我使用该功能来比较时间。

编辑代码以使用 Typed 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 的工作只比在常规 Python 会话中运行以下代码快 1.2 倍：

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 中根本不需要这样做！

于 2013-11-01T11:41:33.667 回答