使用scipys generic_filter实现"峰度过滤器"

Question

我有一个5000*5000numpy数组,我想要计算大小为25的窗口的Kurtosis.我尝试将scipys自己的kurtosis功能放入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_filterscipy函数的正确方法.如果它恰好是正确的,那么它对我来说太慢了.所以我的下一个问题是,如果有更快的方法来实现这一目标？例如,考虑标准偏差,您可以简单地执行以下操作:

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 $的事实.

Answer 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

考虑到任务的大小,这是一个非常合理的表现!