python - 使用 scipys generic_filter 实现 "Kurtosis filter"

Question

我有一个 5000*5000 numpy 数组，我想在其上计算大小为 25 的窗口的峰度。我尝试将 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$ 这一事实。

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 ~ 470 亿 (!) 乘法。这甚至会太慢而无法在普通 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

鉴于任务的规模，这是一个非常合理的表现!