使用 scipys generic_filter 实现“峰度过滤器”

Posted 2023-02-14

【中文标题】使用 scipys generic_filter 实现“峰度过滤器”

【英文标题】：Implementing a "Kurtosis filter" using scipys generic_filter

【发布时间】：2014-08-09 09:05:09

【问题描述】：

我有一个5000*5000 numpy 数组，我想在该数组上计算大小为 25 的窗口的峰度。我尝试将 scipys 自己的峰度函数放在 ndimage.filters 中的 generic_filter 中，如下所示：

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$。

【问题讨论】：

您需要注意与您提出的其他方法类似的方法的数值稳定性，尤其是对于涉及 4 次方的峰态。 pandas 有一个滚动峰态函数，pd.stats.moments.rolling_kurt`，但实现也没有做好稳定工作，它只能在一个维度上工作......

你需要第四个时刻来计算峰度。您可以这样计算，峰度 = mu_4/sigma^4 - 3。sigma 是标准差，mu_4 是均值附近的第四个矩。

关键字是“围绕均值” - 从非居中的四阶矩（在滚动窗口样式中很容易获得）到居中的四阶矩不太容易，如问题中所述，而不是从非中心二阶矩到中心二阶矩（您必须编写中心版本的完整多项式展开式）。

【参考方案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

考虑到任务的大小，性能非常合理！

【讨论】：

感谢您的精彩回答。

我目前正在尝试在 Windows 64 位工作站上执行此操作。我安装了蟒蛇。当我执行“python setup.py build_ext --inplace”时，出现以下错误：致命错误：numpy/arrayobject.h：没有这样的文件或目录。编译终止。错误：命令“gcc”失败，退出状态为 1。关于我可以做些什么来解决这个问题的任何线索？再次感谢。

我把 setup.py 换成了： from distutils.core import setup from distutils.core import Extension import numpy setup( ext_modules=[ Extension("_kurtosis", ["_kurtosismodule.c"], include_dirs= [numpy.get_include()]), ], )

我不能再编译了。我正在为 python 2.7 使用 VC：microsoft.com/en-us/download/details.aspx?id=44266 但不记得我以前用过什么。会是这样吗？