scipy.stats.t.sfとCythonを使用したGSLの比較

Question

大きな2Dのnumpy値配列のp値を計算したいと思います。しかし、これには時間がかかり、スピードを上げたいと思います。 GSLを使ってみました。 単一のgsl_cdf_tdist_Pはscipy.stats.t.sfよりはるかに高速ですが、ndarrayを反復処理すると処理が非常に遅くなります。私はこれを改善するために援助をしたいと思います。 下記のコードを参照してください。

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 

ipython

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) 

UPDATE：が、私はまだscipyのダウンロードよりも低い計算時間を短縮することができましたではなく、CDEF宣言を追加します。私は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 

Answer 1

代わりstats.t.sfの生の特殊機能を使用することにより、生のパフォーマンスのいくつかの小さなゲインを得ることができます。ソースを見ると、あなたが直接stdtrを使用することができます（https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L3849）

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

ので見つける：あなたはcythonのために手を差し伸べるない場合

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 

は、入力されたmemoryview構文は、多くの場合、あなたはより高速なコードを提供します古いndarray構文：ここでは私もののdevのバージョンでのみavaialbleですcython_specialインタフェースを、使用してい

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) 

scipy（http://scipy.github.io/devdocs/special.cython_special.html#module-scipy.special.cython_special）ですが、必要に応じてGSLを使用できます。

最後に、反復でボトルネックが疑われる場合は、cython -aの出力を調べて、ホットループにPythonオーバーヘッドがあるかどうかを確認してください。