NumPy函数的内存消耗为标准偏差

Question

我目前正在使用GDAL的Python绑定来处理相当大的栅格数据集(> 4 GB).因为一次将它们加载到内存中对我来说是不可行的解决方案我将它们读成更小的块并逐个进行计算.为了避免每个块读取的新分配,我使用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

在GitHub上的NumPy源代码树中进行的搜索显示该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中是否有一些其他函数利用具有恒定内存消耗的方法,如Welford算法(维基百科),用于平均值和标准差的一次通过计算？

另一种方法是使用预分配缓冲区_var的可选out参数(如NumPy ufuncs)实现该函数的自定义版本.使用这种方法,不会消除附加副本,但至少内存消耗将是恒定的,并且保存每个块中的分配的运行时间.

编辑:测试kezzos建议的Welford算法的Cython实现.

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实现(可以在一个函数中计算mean和std):

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的迭代方法对于较小的数据集更快,而对于具有10000+元素的更大数据集,NumPy向量(可能是SIMD加速)方法更快.所有测试均使用Python 2.7.9和NumPy 1.9.2版进行.

请注意,在实际情况下,上层函数将用于计算栅格的单个块的统计信息.所有块的标准偏差和均值将与维基百科(此处)中建议的方法结合使用.它的优点是不需要对光栅的所有元素求和,从而避免浮点溢出问题(至少在某一点上).

Answer 1

我怀疑你会发现任何这样的功能numpy.存在的理由numpy是,它利用了矢量处理器指令集 - 执行大量数据的相同指令.基本上numpy换取内存效率以提高速度效率.但是,由于Python的内存密集型特性,numpy通过将数据类型与整个数组相关联而不是每个单独的元素,也能够实现某些内存效率.

提高速度但仍牺牲一些内存开销的一种方法是计算块中的标准偏差,例如.

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

显示加速 - 越大blocksize,加速越大.并且大大降低了内存开销.并非完全较低的内存开销是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

Answer 2

Cython来救援!这实现了很好的加速:

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

用它来测试:

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

Cython代码

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

原始代码

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

所以速度提高了大约300倍.仍然没有numpy版本那么好..