Speicherverbrauch der NumPy-Funktion für Standardabweichung

Question

Ich verwende derzeit die Python-Bindungen von GDAL, um mit ziemlich großen Raster-Datensätzen (> 4 GB) zu arbeiten. Da das Laden in den Speicher auf einmal keine praktikable Lösung für mich ist, lese ich sie in kleinere Blöcke und führe die Berechnungen Stück für Stück durch. Um eine neue Zuordnung für jeden gelesenen Block zu vermeiden, verwende ich das buf_obj Argument (here), um die Werte in ein vorbelegtes NumPy-Array einzulesen. An einem Punkt muss ich die mittlere und Standardabweichung des gesamten Rasters berechnen. Natürlich habe ich np.std für die Berechnung verwendet. Durch das Profiling des Speicherverbrauchs meines Programms habe ich jedoch festgestellt, dass bei jedem Aufruf von np.std zusätzlich Speicher zugewiesen und freigegeben wird.

Ein Minimum Beispiel arbeitet, die dieses Verhalten zeigt:

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 ergab im Quellbaum von NumPy auf GitHub, dass die np.std Funktion intern die _var Funktion von _methods.py ruft (here). An einem Punkt berechnet _var die Abweichungen vom Mittelwert und summiert sie. Daher wird eine temporäre Kopie des Eingabearrays erstellt. Die Funktion berechnet im Wesentlichen die Standardabweichung wie folgt:

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

Während dieser Ansatz OK für kleinere Eingabe-Arrays ist es definitiv keine Möglichkeit für größere zu gehen. Da ich, wie bereits erwähnt, kleinere Speicherblöcke benutze, ist diese zusätzliche Kopie in meinem Programm kein Speicherproblem. Was mich jedoch stört ist, dass für jeden Block eine neue Zuweisung gemacht und freigegeben wird, bevor der nächste Block gelesen wird.

Gibt es eine andere Funktion innerhalb von NumPy oder SciPy, die einen Ansatz mit konstantem Speicherverbrauch wie den Welford-Algorithmus (Wikipedia) für die Berechnung von Mittelwert und Standardabweichung verwendet?

Ein anderer Weg zu gehen wäre, eine benutzerdefinierte Version der _var Funktion mit einem optionalen out Argument für einen vorab zugeordneten Puffer (wie NumPy ufuncs) zu implementieren. Bei diesem Ansatz würde die zusätzliche Kopie nicht eliminiert werden, aber zumindest wäre der Speicherverbrauch konstant und die Laufzeit für die Zuweisungen in jedem Block wird gespeichert.

EDIT: Getestet die Cython-Implementierung des Welford-Algorithmus wie von Kezzos vorgeschlagen.

Cython Umsetzung (von kezzos modifiziert):

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 Umsetzung (Mittelwert und Standard kann 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) 

Prüfergebnisse eine Zufallsdatensatzes unter Verwendung von (Zeiten in Millisekunden, Best of 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 wie die iterative a pproach mit Cython ist schneller mit kleineren Datensätzen und der NumPy Vektor (möglicherweise SIMD beschleunigt) Ansatz für größere Datensätze mit 10000 Elemente. Alle Tests wurden mit Python 2.7.9 und NumPy Version 1.9.2 durchgeführt.

Beachten Sie, dass im realen Fall zu den oberen Funktionen würden verwendet werden, um die Statistiken für einen einzelnen Block des Rasters zu berechnen.Die Standardabweichungen und Mittelwerte für alle Blöcke sind mit der in Wikipedia vorgeschlagenen Methode (here) zu kombinieren. Es hat den Vorteil, dass nicht alle Elemente des Rasters summiert werden müssen und somit das Float-Overflow-Problem (zumindest bis zu einem gewissen Punkt) vermieden wird.

Answer 1

Cython zur Rettung! Dadurch wird eine schöne beschleunigen:

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

@cython.boundscheck(False) 
def std_welford(np.ndarray[np.float64_t, ndim=1] a): 
    cdef int n = 0 
    cdef float mean = 0 
    cdef float M2 = 0 
    cdef int a_len = len(a) 
    cdef int i 
    cdef float delta 
    cdef float result 
    for i in range(a_len): 
     n += 1 
     delta = a[i] - mean 
     mean += delta/n 
     M2 += delta * (a[i] - mean) 
    if n < 2: 
     result = np.nan 
     return result 
    else: 
     result = sqrt(M2/(n - 1)) 
     return result 

diesen Test zu verwenden:

a = np.random.rand(10000).astype(np.float) 
print std_welford(a) 
%timeit -n 10 -r 10 std_welford(a) 

Cython Code

0.288327455521 
10 loops, best of 10: 59.6 µs per loop 

Originalcode

0.289605617397 
10 loops, best of 10: 18.5 ms per loop 

Numpy std

0.289493223504 
10 loops, best of 10: 29.3 µs per loop 

So eine Geschwindigkeitssteigerung von 300x um. Immer noch nicht so gut wie die numpy Version ..

Answer 2

Ich bezweifle, dass Sie solche Funktionen in numpy finden werden. Die Existenzberechtigung von numpy ist, dass es vector processor Befehlssätze nutzt - die gleiche Anweisung von großen Datenmengen auszuführen. Grundsätzlich numpy Handel Speichereffizienz für Geschwindigkeitseffizienz. Aufgrund der speicherintensiven Natur von Python ist numpy jedoch auch in der Lage, bestimmte Speichereffizienzen zu erzielen, indem der Datentyp dem Array als Ganzes und nicht jedem einzelnen Element zugeordnet wird.

Ein Weg, um die Geschwindigkeit zu verbessern, aber immer noch etwas Speicheraufwand opfern, ist die Standardabweichung in Chunks berechnen z.

import numpy as np 

def std(arr, blocksize=1000000): 
    """Written for py3, change range to xrange for py2. 
    This implementation requires the entire array in memory, but it shows how you can 
    calculate the standard deviation in a piecemeal way. 
    """ 
    num_blocks, remainder = divmod(len(arr), blocksize) 
    mean = arr.mean() 
    tmp = np.empty(blocksize, dtype=float) 
    total_squares = 0 
    for start in range(0, blocksize*num_blocks, blocksize): 
     # get a view of the data we want -- views do not "own" the data they point to 
     # -- they have minimal memory overhead 
     view = arr[start:start+blocksize] 
     # # inplace operations prevent a new array from being created 
     np.subtract(view, mean, out=tmp) 
     tmp *= tmp 
     total_squares += tmp.sum() 
    if remainder: 
     # len(arr) % blocksize != 0 and need process last part of array 
     # create copy of view, with the smallest amount of new memory allocation possible 
     # -- one more array *view* 
     view = arr[-remainder:] 
     tmp = tmp[-remainder:] 
     np.subtract(view, mean, out=tmp) 
     tmp *= tmp 
     total_squares += tmp.sum() 

    var = total_squares/len(arr) 
    sd = var ** 0.5 
    return sd 

a = np.arange(20e6) 
assert np.isclose(np.std(a), std(a)) 

die Geschwindigkeit auftauchend --- desto größer ist die blocksize, desto größer ist die Geschwindigkeit auf. Und erheblich weniger Speicheraufwand. Nicht ganz der niedrigere Speicheraufwand ist 100% genau.

In [70]: %timeit np.std(a) 
10 loops, best of 3: 105 ms per loop 

In [71]: %timeit std(a, blocksize=4096) 
10 loops, best of 3: 160 ms per loop 

In [72]: %timeit std(a, blocksize=1000000) 
10 loops, best of 3: 105 ms per loop 

In [73]: %memit std(a, blocksize=4096) 
peak memory: 360.11 MiB, increment: 0.00 MiB 

In [74]: %memit std(a, blocksize=1000000) 
peak memory: 360.11 MiB, increment: 0.00 MiB 

In [75]: %memit np.std(a) 
peak memory: 512.70 MiB, increment: 152.59 MiB