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.