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

У меня есть массив 5000*5000 numpy, на котором я хочу рассчитать Kurtosis для окон размером 25. Я попытался поместить функцию curtsus scipys в 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 $.

Ответ 1

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

def kurt(X, w):
    n, m = X.shape
    K = np.zeros_like(X)

    for i in xrange(w, n-w):                       # 5000 iterations
        for j in xrange(w, m-w):                   # 5000 iterations
            x = X[i-w:i+w+1,j-w: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, не говоря уже о передаче функции 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 = s.T.copy()           # 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(n-W):
            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 < n-W; ++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

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