В настоящее время я использую привязки Python GDAL для работы с довольно большими наборами растровых данных (> 4 ГБ). Поскольку загрузка их в память для меня нереально, я читаю их в меньшие блоки и выполняю вычисления по частям. Чтобы избежать нового выделения для каждого прочитанного блока, я используюbuf_obj аргумент (Вот) для чтения значений в предварительно выделенный массив NumPy. В какой-то момент мне нужно вычислить среднее и стандартное отклонение всего растра. Естественно я использовалnp.std для расчета. Однако, профилируя потребление памяти моей программой, я понял, что при каждом вызовеnp.std дополнительно выделяется и освобождается память.

Минимальный рабочий пример, который демонстрирует это поведение:

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

Поиск в дереве исходников NumPy на GitHub показал, чтоnp.std Функция внутренне вызывает_var функция от_methods.py (Вот). В одной точке_var вычисляет отклонения от среднего и суммирует их. Поэтому создается временная копия входного массива. Функция по существу вычисляет стандартное отклонение следующим образом:

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

Хотя этот подход подходит для меньших входных массивов, он явно не подходит для больших массивов. Поскольку я использую меньшие блоки памяти, как упоминалось ранее, эта дополнительная копия не является проблемой с точки зрения памяти в моей программе. Однако, что меня беспокоит, так это то, что для каждого блока производится новое выделение и освобождается перед чтением следующего блока.

Есть ли какая-то другая функция в NumPy или SciPy, которая использует подход с постоянным потреблением памяти, такой как алгоритм Уэлфорда (Википедия) за один проход вычисление среднего значения и стандартного отклонения?

Другим способом было бы реализовать пользовательскую версию_var функция с дополнительнымout аргумент для предварительно выделенного буфера (например, NumPy ufuncs). При таком подходе дополнительная копия не будет исключена, но по крайней мере потребление памяти будет постоянным, а время выполнения для выделений в каждом блоке будет сохранено.

РЕДАКТИРОВАТЬ: Протестировал реализацию Cython алгоритма Уэлфорда, как это было предложено kezzos.

Реализация на Cython (модифицированная из 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 (среднее и стандартное вычисление могут быть вычислены в одной функции):

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

Результаты теста с использованием случайного набора данных (время в миллисекундах, лучше всего из 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 |
----------------------------------

Кажется, что итеративный подход с использованием Cython быстрее с меньшими наборами данных и векторным подходом NumPy (возможно, с ускорением SIMD) для больших наборов данных с более чем 10000 элементами. Все тесты проводились на Python 2.7.9 и NumPy версии 1.9.2.

Обратите внимание, что в реальном случае функции верхнего уровня будут использоваться для вычисления статистики для одного блока растра. Стандартные отклонения и средства для всех блоков должны быть объединены с методологией, предложенной в Википедии (Вот). Преимущество состоит в том, что не все элементы растра необходимо суммировать, что позволяет избежать проблемы переполнения поплавка (по крайней мере, до некоторой точки).