Speicherverbrauch der NumPy-Funktion für Standardabweichung

Ich verwende derzeit die Python-Bindungen von GDAL, um mit recht großen Rasterdatensätzen (> 4 GB) zu arbeiten. Da es für mich keine praktikable Lösung ist, sie sofort in den Speicher zu laden, lese ich sie in kleinere Blöcke und berechne sie Stück für Stück. Um eine neue Zuordnung für jeden gelesenen Block zu vermeiden, verwende ich dasbuf_obj Streit Hie), um die Werte in ein vorab zugewiesenes NumPy-Array einzulesen. An einem Punkt muss ich den Mittelwert und die Standardabweichung des gesamten Rasters berechnen. Natürlich habe ich @ verwendnp.std für die Berechnung. Durch die Profilerstellung des Speicherverbrauchs meines Programms wurde mir jedoch klar, dass bei jedem Aufruf vonnp.std zusätzlich wird Speicher reserviert und freigegeben.

Ein minimales Arbeitsbeispiel, das dieses Verhalten demonstriert:

In [1]  import numpy as np
In [2]  a = np.random.rand(20e6)  # Approx. 150 MiB of memory
In [3]  %memit np.mean(a)
peak memory: 187.30 MiB, increment: 0.48 MiB
In [4]  %memit np.std(a)
peak memory: 340.24 MiB, increment: 152.91 MiB

Eine Suche im Quelltext von NumPy auf GitHub ergab, dass dasnp.std -Funktion ruft intern das @ a_var Funktion von_methods.py (Hie). An einer Stelle_var berechnet die Abweichungen vom Mittelwert und summiert sie auf. Daher wird eine temporäre Kopie des Eingabearrays erstellt. Die Funktion berechnet die Standardabweichung im Wesentlichen wie folgt:

mu = sum(arr) / len(arr)
tmp = arr - mu
tmp = tmp * tmp
sd = np.sum(tmp) / len(arr)

Während dieser Ansatz für kleinere Eingabearrays in Ordnung ist, ist er für größere definitiv nicht geeignet. Da ich, wie bereits erwähnt, kleinere Speicherblöcke verwende, ist diese zusätzliche Kopie aus Speichersicht in meinem Programm kein Problem, das das Spiel bricht. Was mich jedoch stört, ist, dass für jeden Block eine neue Zuordnung vorgenommen und freigegeben wird, bevor der nächste Block gelesen wird.

ibt es eine andere Funktion in NumPy oder SciPy, die einen Ansatz mit konstantem Speicherverbrauch verwendet, wie der Welford-Algorithmus Wikipedia) für eine Durchgangsberechnung von Mittelwert und Standardabweichung?

Ein weiterer Weg wäre die Implementierung einer benutzerdefinierten Version des_var Funktion mit einem optionalenout Argument für einen vorab zugewiesenen Puffer (wie NumPy ufuncs). Mit diesem Ansatz würde die zusätzliche Kopie nicht eliminiert, aber zumindest der Speicherverbrauch wäre konstant und die Laufzeit für die Zuweisungen in jedem Block wird gespeicher

BEARBEITEN Die von kezzos vorgeschlagene Cython-Implementierung des Welford-Algorithmus wurde getestet.

Cython-Implementierung (geändert von kezzos):

cimport cython
cimport numpy as np
from libc.math cimport sqrt

@cython.boundscheck(False)
def iterative_approach(np.ndarray[np.float32_t, ndim=1] a):
    cdef long n = 0
    cdef float mean = 0
    cdef float M2 = 0
    cdef long i
    cdef float delta
    cdef float a_min = 10000000  # Must be set to Inf and -Inf for real cases
    cdef float a_max = -10000000
    for i in range(len(a)):
        n += 1
        delta = a[i] - mean
        mean += delta / n
        M2 += delta * (a[i] - mean)
        if a[i] < a_min:
            a_min = a[i]
        if a[i] > a_max:
            a_max = a[i]
    return a_min, a_max, mean, sqrt(M2 / (n - 1))

NumPy-Implementierung (mean und std können möglicherweise in einer Funktion berechnet werden):

def vector_approach(a):
    return np.min(a), np.max(a), np.mean(a), np.std(a, ddof=1)

Testergebnisse unter Verwendung eines zufälligen Datensatzes (Zeiten in Millisekunden, beste von 25):

----------------------------------
| Size |  Iterative |     Vector |
----------------------------------
|  1e2 |    0.00529 |    0.17149 |
|  1e3 |    0.02027 |    0.16856 |
|  1e4 |    0.17850 |    0.23069 |
|  1e5 |    1.93980 |    0.77727 |
|  1e6 |   18.78207 |    8.83245 |
|  1e7 |  180.04069 |  101.14722 |
|  1e8 | 1789.60228 | 1086.66737 |
----------------------------------

Es scheint, als ob der iterative Ansatz mit Cython bei kleineren Datensätzen und der NumPy-Vektor-Ansatz (möglicherweise SIMD-beschleunigt) bei größeren Datensätzen mit mehr als 10000 Elementen schneller ist. Alle Tests wurden mit Python 2.7.9 und NumPy Version 1.9.2 durchgeführt.

Beachten Sie, dass im realen Fall zur Berechnung der Statistiken für einen einzelnen Block des Rasters die oberen Funktionen verwendet werden. Die Standardabweichungen und Mittelwerte für alle Blöcke sind mit der in Wikipedia vorgeschlagenen Methodik zu kombinieren Hie). Dies hat den Vorteil, dass nicht alle Elemente des Rasters aufsummiert werden müssen, wodurch das Problem des Überlaufs des Floats (zumindest bis zu einem gewissen Punkt) vermieden wird.