Atualmente, estou usando as ligações Python do GDAL para trabalhar em conjuntos de dados raster muito grandes (> 4 GB). Como carregá-los na memória de uma vez não é uma solução viável para mim, eu os leio em blocos menores e faço os cálculos peça por peça. Para evitar uma nova alocação para cada bloco lido, estou usando obuf_obj argumento (aqui) para ler os valores em uma matriz NumPy pré-alocada. Em um ponto, eu tenho que calcular a média e o desvio padrão de toda a varredura. Naturalmente eu useinp.std para o cálculo. No entanto, ao traçar o perfil do consumo de memória do meu programa, percebi que a cada chamada denp.std adicionalmente, a memória é alocada e liberada.

Um exemplo de trabalho mínimo que demonstra esse 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

Uma pesquisa na árvore de origem do NumPy no GitHub revelou que onp.std função invoca internamente o_var função de_methods.py (aqui) Em um ponto_var calcula os desvios da média e os soma. Portanto, uma cópia temporária da matriz de entrada é criada. A função calcula essencialmente o desvio padrão da seguinte maneira:

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

Embora essa abordagem seja adequada para matrizes de entrada menores, definitivamente não há como ir para matrizes maiores. Como estou usando blocos de memória menores, como mencionado anteriormente, esta cópia adicional não é um problema de quebra de jogo do ponto de vista da memória no meu programa. No entanto, o que me incomoda é que, para cada bloco, uma nova alocação seja feita e liberada antes de ler o próximo bloco.

Existe alguma outra função no NumPy ou SciPy que utiliza uma abordagem com consumo constante de memória, como o algoritmo Welford (Wikipedia) para uma passagem de cálculo da média e desvio padrão?

Outro caminho a seguir seria implementar uma versão personalizada do_var função com um opcionalout argumento para um buffer pré-alocado (como NumPy ufuncs). Com essa abordagem, a cópia adicional não seria eliminada, mas pelo menos o consumo de memória seria constante e o tempo de execução das alocações em cada bloco é salvo.

EDITAR: Testou a implementação Cython do algoritmo Welford, conforme sugerido por kezzos.

Implementação do Cython (modificada a partir de 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))

Implementação do NumPy (a média e o std podem ser calculados em uma função):

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

Resultados do teste usando um conjunto de dados aleatórios (vezes em milissegundos, melhor de 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 |
----------------------------------

Parece que a abordagem iterativa usando o Cython é mais rápida com conjuntos de dados menores e a abordagem do vetor NumPy (possivelmente acelerada pelo SIMD) para conjuntos de dados maiores com mais de 10000 elementos. Todos os testes foram conduzidos com o Python 2.7.9 e o NumPy versão 1.9.2.

Observe que, no caso real, as funções superiores seriam usadas para calcular as estatísticas de um único bloco da varredura. Os desvios padrão e os meios para todos os blocos devem ser combinados com a metodologia sugerida na Wikipedia (aqui) Tem a vantagem de que nem todos os elementos da varredura precisam ser resumidos e, assim, evitam o problema de estouro de flutuação (pelo menos até certo ponto).