ein „Kurtosis-Filter“ mit scipys Implementierung generic_filter

Question

ich einen 5000*5000 numpy Array, auf das ich die Kurtosis für Fenster der Größe berechnen möge 25. Ich in den generic_filter eigene Kurtosis Funktion scipys versucht setzen wie so in ndimage.filters gefunden:

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') 

Dies endet nie und ich bin mir nicht sicher, ob es die richtige Antwort gibt. Also meine erste Frage ist, ob dies eine korrekte Möglichkeit ist, die generic_filter mit einer scipy-Funktion zu verwenden. Wenn es richtig war, dann ist es zu langsam, als daß es mir nützen könnte. Meine nächste Frage wäre also, ob es einen schnelleren Weg dafür gibt? Zum Beispiel, um eine Standardabweichung denken können Sie einfach so etwas wie:

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 

Dies wird schnell pralle und kommen einfach aus der Tatsache, dass $ \ sigma^2 = E [(X - \ mu)^2] = E [X^2] - (E [X])^2 $.

Answer 1

Ihr Ansatz ist korrekt, aber wie Sie bemerken, ist es für die anstehende Aufgabe viel zu langsam. Überlegen Sie, wie groß Ihre Aufgabe in der numerisch beste Umsetzung ist (nicht über Grenzwerte stört):

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 

So haben wir 5000*5000*1875 ~ 47 billion Multiplikationen (!). Dies ist sogar zu langsam, um in einer einfachen C-Implementierung nützlich zu sein, geschweige denn, indem eine Python-Funktion kurtosis() an die innere Schleife von generic_filter() übergeben wird. Letzteres ruft tatsächlich eine C-Erweiterungsfunktion auf, aber es gibt vernachlässigbare Vorteile, da es bei jeder Iteration zurück in Python aufgerufen werden muss, was sehr teuer ist.

Also, das eigentliche Problem ist, dass Sie einen besseren Algorithmus benötigen. Da scipy es nicht hat, lasst es uns hier Schritt für Schritt entwickeln.

Die Hauptbeobachtung, die eine Beschleunigung dieses Problems ermöglicht, besteht darin, dass die Kurtosis-Berechnungen für aufeinanderfolgende Fenster auf größtenteils den gleichen Werten basieren, mit Ausnahme einer Zeile (25 Werte), die ersetzt wird. Statt also die Kurtosis mit allen 625 Werten von Grund auf neu zu berechnen, versuchen wir, zuvor berechnete Summen zu verfolgen und sie so zu aktualisieren, dass nur die 25 neuen Werte verarbeitet werden müssen.

Dies erfordert den (x - mu)**4 Faktor erweitert, da nur die laufenden Summen über x, x**2, x**3 und x**4 können leicht aktualisiert werden. Es gibt keine schöne Stornierung wie in der Formel für die Standardabweichung, die Sie erwähnen, aber es ist durchaus möglich:

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 

Hinweis:Der Algorithmus in dieser Form geschrieben ist numerisch weniger stabil, da wir Zähler lassen und Nenner werden einzeln sehr groß, während wir uns vorher schon früh geteilt haben, um dies zu verhindern (selbst auf Kosten eines sqrt). Ich fand jedoch, dass dies für die Kurtosis nie ein Problem für praktische Anwendungen war.

Im obigen Code habe ich versucht, die Anzahl der Multiplikationen zu minimieren. Die läuft bedeutetM1, M2, M3 und M4 kann jetzt ziemlich einfach aktualisiert werden, durch Subtraktion der Beiträge der Zeile, die nicht länger Teil des Fensters und Hinzufügen der Beiträge der neuen Zeile.

Lassen Sie uns dies umzusetzen:

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 

Nun, da wir einen guten Algorithmus haben, stellen wir fest, dass die Timing-Ergebnisse noch eher enttäuschend sind. Unser Problem ist jetzt, dass Python + numpy die falsche Sprache für solch eine Nummerierung ist. Lass uns eine C-Erweiterung schreiben! Hier ist _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(); 
} 

Build with:

python setup.py build_ext --inplace 

wo setup.py ist:

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

Bitte beachte, dass wir keine Speicher in der C-Erweiterung zuordnen. Auf diese Weise müssen wir uns nicht mit der Anzahl der Referenzen/Garbage Collection herumschlagen. Wir haben gerade einen Einstieg in Python verwenden:

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 

Schließlich lassen Sie uns das Timing tun:

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 

Eine sehr vernünftige Leistung die Größe der Aufgabe gegeben!