使用 Cython 比较 scipy.stats.t.sf 与 GSL答案

【问题标题】：Comparing scipy.stats.t.sf vs GSL using Cython使用 Cython 比较 scipy.stats.t.sf 与 GSL
【发布时间】：2016-10-05 22:26:32
【问题描述】：

我想计算一个大型 2D numpy t 值数组的 p 值。但是，这需要很长时间，我想提高它的速度。我尝试使用 GSL。尽管单个 gsl_cdf_tdist_P 比 scipy.stats.t.sf 快得多，但在遍历 ndarray 时，该过程非常缓慢。我想帮助改善这一点。请参阅下面的代码。

GSL_Test.pyx

import cython
cimport cython

import numpy
cimport numpy
from cython_gsl cimport *

DTYPE = numpy.float32
ctypedef numpy.float32_t DTYPE_t

cdef get_gsl_p(double t, double nu):
    return (1 - gsl_cdf_tdist_P(t, nu)) * 2

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef get_gsl_p_for_2D_matrix(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):

    cdef unsigned int rows = t_matrix.shape[0]
    cdef numpy.ndarray[DTYPE_t, ndim=2] out = numpy.zeros((rows, rows), dtype='float32')
    cdef unsigned int row, col

    for row in range(rows):
        for col in range(rows):
            out[row, col] = get_gsl_p(t_matrix[row, col], n-2)
    return out

def get_gsl_p_for_2D_matrix_def(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):
    return get_gsl_p_for_2D_matrix(t_matrix, n)

ipython

import GSL_Test
import numpy
import scipy.stats
a = numpy.random.rand(3544, 3544).astype('float32')
%timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix(a, 25)
1 loop, best of 3: 7.87 s per loop
%timeit -n 1 scipy.stats.t.sf(a, 25)*2
1 loop, best of 3: 4.66 s per loop

更新：添加 cdef 声明我能够减少计算时间，但仍不低于 scipy。我修改了代码以具有 cdef 声明。

%timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix_def(a, 25)
1 loop, best of 3: 6.73 s per loop

【问题讨论】：

你对stats.t.sf了解多少？为什么你认为cython 版本会更快？我寻找可以翻译成c的迭代之类的东西。
scipy.stats.t.sf 和其他方法有一些 python 开销，主要用于输入检查。如果你想要速度，那么你可以直接使用用 Fortran 或 C 编写的 scipy.special 函数。我怀疑有可能以任何可观的数量击败 scipy.special 除非有更有效的算法或交易算法关闭速度以确保准确性。
为 get_gsl_p_for_2D_matrix 使用单独的 def 和 cdef 函数没有任何好处。为 get_gsl_p 指定返回类型可能会有所帮助。

标签： python numpy scipy cython gsl

【解决方案1】：

通过使用原始特殊函数而不是stats.t.sf，您可以获得一些原始性能的小幅提升。看源码你发现(https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L3849)

def _sf(self, x, df):
      return sc.stdtr(df, -x)

这样你就可以直接使用stdtr：

np.random.seed(1234)
x = np.random.random((3740, 374))
t1 = stats.t.sf(x, 25)
t2 = stdtr(25, -x)

1 loop, best of 3: 653 ms per loop
1 loop, best of 3: 562 ms per loop

如果您确实使用 cython，则类型化的 memoryview 语法通常会比旧的 ndarray 语法为您提供更快的代码：

from scipy.special.cython_special cimport stdtr
from numpy cimport npy_intp
import numpy as np

def tsf(double [:, ::1] x, int df=25):
    cdef double[:, ::1] out = np.empty_like(x)
    cdef npy_intp i, j
    cdef double tmp, xx

    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            xx = x[i, j]
            out[i, j] = stdtr(df, -xx)

    return np.asarray(out)

这里我也使用cython_special接口，这个接口只在scipy的开发版（http://scipy.github.io/devdocs/special.cython_special.html#module-scipy.special.cython_special）中可用，但是你可以使用GSL如果你愿意。

最后，如果您怀疑迭代中存在瓶颈，请不要忘记检查 cython -a 的输出以查看热循环中是否存在一些 python 开销。

【讨论】：