Consumo di memoria della funzione NumPy per la deviazione standard

Question

Attualmente sto usando i binding Python di GDAL per lavorare su set di dati raster abbastanza grandi (> 4 GB). Dal momento che caricarli in memoria in una sola volta non è una soluzione fattibile per me li leggo in blocchi più piccoli e faccio i calcoli pezzo dopo pezzo. Per evitare una nuova allocazione per ogni blocco letto sto usando l'argomento buf_obj (here) per leggere i valori in un array NumPy preallocato. A un certo punto devo calcolare la media e la deviazione standard dell'intero raster. Naturalmente ho usato np.std per il calcolo. Tuttavia, analizzando il consumo di memoria del mio programma, mi sono reso conto che con ogni chiamata di np.std la memoria è stata allocata e rilasciata.

Un esempio di lavoro minima che illustra questo comportamento:

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 

Una ricerca all'interno dell'albero sorgente di NumPy su GitHub rivelato che la funzione np.std richiama internamente la funzione _var da _methods.py (here). A un certo punto, _var calcola le deviazioni dalla media e le somma. Pertanto viene creata una copia temporanea dell'array di input. La funzione calcola essenzialmente la deviazione standard come segue:

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

Mentre questo approccio è OK per array in input piccole è sicuramente modo di andare per quelli più grandi. Dato che sto usando piccoli blocchi di memoria come accennato prima, questa copia aggiuntiva non è un problema di gioco dal punto di vista della memoria nel mio programma. Tuttavia, ciò che mi infastidisce è che per ogni blocco viene fatta e rilasciata una nuova allocazione prima di leggere il blocco successivo.

C'è qualche altra funzione in NumPy o SciPy che utilizza un aggiornamento con consumo di memoria costante come l'algoritmo di Welford (Wikipedia) per un calcolo di un passaggio della media e deviazione standard?

Un altro modo per andare sarebbe implementare una versione personalizzata della funzione _var con un argomento opzionale out per un buffer preallocato (come gli ufunc NumPy). Con questo approccio la copia aggiuntiva non verrebbe eliminata, ma almeno il consumo di memoria sarebbe costante e il tempo di esecuzione per le allocazioni in ogni blocco verrà salvato.

MODIFICA: Testata l'implementazione Cython dell'algoritmo di Welford come suggerito da kezzos.

attuazione Cython (modificato da 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)) 

attuazione NumPy (media e std può eventualmente essere calcolato in una funzione):

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

risultati mediante un set di dati casuali (volte in millisecondi, meglio di 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 | 
---------------------------------- 

Sembra il iterativo un pproach con Cython è più veloce con set di dati più piccoli e l'approccio del vettore NumPy (possibilmente SIMD accelerato) per set di dati più grandi con oltre 10000 elementi. Tutti i test sono stati condotti con Python 2.7.9 e NumPy versione 1.9.2.

Nota che nel caso reale alle funzioni superiori verrebbe utilizzato per calcolare le statistiche per un singolo blocco del raster.Le deviazioni standard e i mezzi per tutti i blocchi devono essere combinati con la metodologia suggerita in Wikipedia (here). Ha il vantaggio che non tutti gli elementi del raster devono essere riassunti e quindi evita il problema dell'overflow del float (almeno fino ad un certo punto).

Answer 1

Cython al salvataggio! Questo raggiunge una bella velocità di fino:

%%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 

Usando questo a prova:

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

codice Cython

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

Codice originale

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 

Quindi un aumento di velocità di circa 300x. Ancora non buono come la versione numpy ..

Answer 2

Dubito che troverete tali funzioni in numpy. La ragion d'essere di numpy è che sfrutta i set di istruzioni vector processor - eseguendo la stessa istruzione di grandi quantità di dati. Fondamentalmente lo numpy scambia l'efficienza della memoria per l'efficienza della velocità. Tuttavia, a causa della natura intensiva della memoria di Python, numpy è anche in grado di ottenere determinate efficienze di memoria associando il tipo di dati con l'array nel suo complesso e non ogni singolo elemento.

Un modo per migliorare la velocità, ma comunque sacrificare parte della memoria, è calcolare la deviazione standard in blocchi, ad es.

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

Mostrando la velocità fino --- più grande è il blocksize, maggiore è la velocità del programma. E un sovraccarico di memoria considerevolmente inferiore. Non interamente il sovraccarico di memoria inferiore è accurato al 100%.

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