Внедрение «фильтра Kurtosis» с использованием scipys generic_filter

У меня есть массив 5000*5000 numpy, на котором я хочу рассчитать Kurtosis для окон размером 25. Я попытался использовать функцию generic_filter в generic_filter найденную в ndimage.filters например:

 import numpy as np from scipy.stats import kurtosis from scipy.ndimage.filters import generic_filter mat = np.random.random_sample((5000, 5000)) kurtosis_filter = generic_filter(mat, kurtosis, size=25, mode='reflect') 

Это никогда не заканчивается, и я не уверен, что это дает правильный ответ. Поэтому мой первый вопрос заключается в том, что это правильный способ использования generic_filter с функцией scipy. Если это было правильно, то это слишком медленно для меня, чтобы оно было мне полезно. Итак, мой следующий вопрос будет, если есть более быстрый способ достичь этого? Например, думая о стандартном отклонении, вы можете просто сделать что-то вроде:

 usual_mean = uniform_filter(mat, size=25, mode='reflect') mean_of_squared = uniform_filter(np.multiply(mat,mat), size=25, mode='reflect') standard_deviation = (mean_of_squared - np.multiply(usual_mean,usual_mean))**.5 

Это быстро вспыхивает и просто исходит из того, что $ \ sigma ^ 2 = E [(X – \ mu) ^ 2] = E [X ^ 2] – (E [X]) ^ 2 $.

One Solution collect form web for “Внедрение «фильтра Kurtosis» с использованием scipys generic_filter”

Ваш подход правильный, но, как вы заметили, он слишком медленный для задачи. Подумайте, насколько велика ваша задача в наилучшей по размеру реализации (не беспокоясь о граничных значениях):

 def kurt(X, w): n, m = X.shape K = np.zeros_like(X) for i in xrange(w, nw): # 5000 iterations for j in xrange(w, mw): # 5000 iterations x = X[iw:i+w+1,jw:j+w+1].flatten() # copy 25*25=625 values x -= x.mean() # calculate and subtract mean x /= np.sqrt((x**2).mean()) # normalize by stddev (625 mult.) K[i,j] = (x**4).mean() - 3. # 2*625 = 1250 multiplications return K 

Таким образом, у нас есть 5000*5000*1875 ~ 47 billion (!) Умножений. Это будет слишком медленным, чтобы быть полезным в простой реализации C, не говоря уже о передаче generic_filter() функции Python kurtosis() во внутренний цикл generic_filter() . Последний на самом деле вызывает функцию расширения C, но есть незначительные преимущества, поскольку она должна переходить на Python на каждой итерации, что очень дорого.

Итак, актуальная проблема заключается в том, что вам нужен лучший алгоритм. Поскольку scipy не имеет этого, давайте развить его шаг за шагом здесь.

Главное наблюдение, которое позволяет ускорить эту проблему, состоит в том, что вычисления эксцессов для последовательных окон основаны на основном одних и тех же значениях, за исключением одной строки (25 значений), которая заменяется. Таким образом, вместо пересчета куртоза с нуля, используя все значения 625, мы пытаемся отслеживать ранее рассчитанные суммы и обновлять их таким образом, чтобы обрабатывать только 25 новых значений.

Это требует расширения фактора (x - mu)**4 , так как только текущие суммы по x , x**2 , x**3 и x**4 могут быть легко обновлены. Нет никакой хорошей отмены, как в формуле стандартного отклонения, которую вы упомянули, но это вполне осуществимо:

 def kurt2(X, w): n, m = X.shape K = np.zeros_like(X) W = 2*w + 1 for j in xrange(m-W+1): for i in xrange(n-W+1): x = X[i:i+W,j:j+W].flatten() x2 = x*x x3 = x2*x x4 = x2*x2 M1 = x.mean() M2 = x2.mean() M3 = x3.mean() M4 = x4.mean() M12 = M1*M1 V = M2 - M12; K[w+i,w+j] = (M4 - 4*M1*M3 + 3*M12*(M12 + 2*V)) / (V*V) - 3 return K 

Примечание . Алгоритм, написанный в этой форме, численно менее устойчив, так как мы даем числитель и знаменатель индивидуально очень большими, в то время как раньше мы раньше делились, чтобы предотвратить это (даже ценой sqrt). Однако я обнаружил, что для эксцессов это никогда не было проблемой для практических приложений.

В приведенном выше коде я попытался свести к минимуму количество умножений. Средство запуска M1 , M2 , M3 и M4 теперь может быть легко обновлено путем вычитания вкладов строки, которая больше не является частью окна, и добавления вкладов новой строки.

Давайте реализуем это:

 def kurt3(X, w): n, m = X.shape K = np.zeros_like(X) W = 2*w + 1 N = W*W Xp = np.zeros((4, W, W), dtype=X.dtype) xp = np.zeros((4, W), dtype=X.dtype) for j in xrange(m-W+1): # reinitialize every time we reach row 0 Xp[0] = x1 = X[:W,j:j+W] Xp[1] = x2 = x1*x1 Xp[2] = x3 = x2*x1 Xp[3] = x4 = x2*x2 s = Xp.sum(axis=2) # make sure we sum along the fastest index S = s.sum(axis=1) # the running sums s = sTcopy() # circular buffer of row sums M = S / N M12 = M[0]*M[0] V = M[1] - M12; # kurtosis at row 0 K[w,w+j] = (M[3] - 4*M[0]*M[2] + 3*M12*(M12 + 2*V)) / (V*V) - 3 for i in xrange(nW): xp[0] = x1 = X[i+W,j:j+W] # the next row xp[1] = x2 = x1*x1 xp[2] = x3 = x2*x1 xp[3] = x4 = x2*x2 k = i % W # index in circular buffer S -= s[k] # remove cached contribution of old row s[k] = xp.sum(axis=1) # cache new row S += s[k] # add contributions of new row M = S / N M12 = M[0]*M[0] V = M[1] - M12; # kurtosis at row != 0 K[w+1+i,w+j] = (M[3] - 4*M[0]*M[2] + 3*M12*(M12 + 2*V)) / (V*V) - 3 return K 

Теперь, когда у нас есть хороший алгоритм, мы отмечаем, что результаты синхронизации все еще довольно разочаровывают. Наша проблема заключается в том, что Python + numpy – это неправильный язык для такого количества хрустящих заданий. Давайте напишем расширение C! Вот _kurtosismodule.c :

 #include <Python.h> #include <numpy/arrayobject.h> static inline void add_line(double *b, double *S, const double *x, size_t W) { size_t l; double x1, x2; b[0] = b[1] = b[2] = b[3] = 0.; for (l = 0; l < W; ++l) { b[0] += x1 = x[l]; b[1] += x2 = x1*x1; b[2] += x2*x1; b[3] += x2*x2; } S[0] += b[0]; S[1] += b[1]; S[2] += b[2]; S[3] += b[3]; } static PyObject* py_kurt(PyObject* self, PyObject* args) { PyObject *objK, *objX, *objB; int w; PyArg_ParseTuple(args, "OOOi", &objK, &objX, &objB, &w); double *K = PyArray_DATA(objK); double *X = PyArray_DATA(objX); double *B = PyArray_DATA(objB); size_t n = PyArray_DIM(objX, 0); size_t m = PyArray_DIM(objX, 1); size_t W = 2*w + 1, N = W*W, i, j, k, I, J; double *S = B + 4*W; double *x, *b, M, M2, V; for (j = 0, J = m*w + w; j < m-W+1; ++j, ++J) { S[0] = S[1] = S[2] = S[3] = 0.; for (k = 0, x = X + j, b = B; k < W; ++k, x += m, b += 4) { add_line(b, S, x, W); } M = S[0] / N; M2 = M*M; V = S[1] / N - M2; K[J] = ((S[3] - 4*M*S[2]) / N + 3*M2*(M2 + 2*V)) / (V*V) - 3; for (i = 0, I = J + m; i < nW; ++i, x += m, I += m) { b = B + 4*(i % W); // row in circular buffer S[0] -= b[0]; S[1] -= b[1]; S[2] -= b[2]; S[3] -= b[3]; add_line(b, S, x, W); M = S[0] / N; M2 = M*M; V = S[1] / N - M2; K[I] = ((S[3] - 4*M*S[2]) / N + 3*M2*(M2 + 2*V)) / (V*V) - 3; } } Py_RETURN_NONE; } static PyMethodDef methods[] = { {"kurt", py_kurt, METH_VARARGS, ""}, {0} }; PyMODINIT_FUNC init_kurtosis(void) { Py_InitModule("_kurtosis", methods); import_array(); } 

Построить с помощью:

 python setup.py build_ext --inplace 

где setup.py :

 from distutils.core import setup, Extension module = Extension('_kurtosis', sources=['_kurtosismodule.c']) setup(ext_modules=[module]) 

Обратите внимание: мы не выделяем память в расширении C. Таким образом, нам не нужно возиться со сбором ссылок / сбором мусора. Мы просто используем точку входа в Python:

 import _kurtosis def kurt4(X, w): # add type/size checking if you like K = np.zeros(X.shape, np.double) scratch = np.zeros(8*(w + 1), np.double) _kurtosis.kurt(K, X, scratch, w) return K 

Наконец, давайте сделаем выбор времени:

 In [1]: mat = np.random.random_sample((5000, 5000)) In [2]: %timeit K = kurt4(mat, 12) # 2*12 + 1 = 25 1 loops, best of 3: 5.25 s per loop 

Очень разумная производительность, учитывая размер задачи!

  • scipy.sparse.coo_matrix, как быстро найти нулевой столбец, заполнить 1 и нормализовать
  • Как эффективно вычислять огромное умножение матрицы (функции tfidf) в Python?
  • Многомерный нормальный CDF в Python с использованием scipy
  • scipy.misc.imshow RuntimeError ('Не удалось выполнить просмотр изображения')
  • правильный и эффективный способ сглаживания массива в numpy в python?
  • интерполяция на основе значений одного массива
  • Как ускорить создание матрицы перехода в Numpy?
  • преобразование списка строк python в их тип
  • Python - лучший язык программирования в мире.