【问题标题】:Memory consumption of NumPy function for standard deviation标准差的 NumPy 函数的内存消耗
【发布时间】:2015-08-21 08:19:38
【问题描述】:

我目前正在使用 GDAL 的 Python 绑定来处理相当大的栅格数据集 (> 4 GB)。由于一次将它们加载到内存中对我来说不是可行的解决方案,我将它们读入较小的块并逐个进行计算。为了避免每次读取的块都分配新的分配,我使用 buf_obj 参数 (here) 将值读入预分配的 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 函数在内部调用来自 _methods.py (here) 的 _var 函数。在某一时刻_var 计算与平均值的偏差并将它们相加。因此创建了输入数组的临时副本。该函数本质上计算标准差如下:

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

虽然这种方法适用于较小的输入数组,但绝对不适用于较大的输入数组。因为我使用的是前面提到的较小的内存块,所以从我的程序的内存角度来看,这个额外的副本并不是一个破坏游戏的问题。然而让我烦恼的是,对于每个块,在读取下一个块之前都会进行新的分配并释放。

在 NumPy 或 SciPy 中是否有其他函数利用诸如 Welford 算法 (Wikipedia) 等具有恒定内存消耗的方法来一次性计算均值和标准差?

另一种方法是实现 _var 函数的自定义版本,并为预分配的缓冲区(如 NumPy ufunc)提供可选的 out 参数。使用这种方法不会消除额外的副本,但至少内存消耗是恒定的,并且每个块中分配的运行时间都被节省了。

编辑:按照 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 实现(可以在一个函数中计算均值和标准差):

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 版进行。

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

【问题讨论】:

  • 您是否尝试过在缓冲区对象上使用 Welford 算法(在纯 Python 中)?
  • @kezzos 我已经用纯 Python 实现的测试结果更新了这个问题。 Python 实现比 NumPy 版本慢很多。

标签: python numpy memory scipy


【解决方案1】:

我怀疑你会在numpy 中找到任何此类功能。 numpy 存在的理由是它利用了vector processor 指令集——对大量数据执行相同的指令。基本上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 [75]: %memit np.std(a)
peak memory: 512.70 MiB, increment: 152.59 MiB

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

【讨论】:

  • 实际上,长期以来,SIMD(向量)指令仅由 Numpy 在调用 BLAS / LAPACK 函数时使用(取决于后端)。只有从 v1.6 开始,Numpy 中才有实际的 SSE 代码,用于einsum。并且仅从 v1.8 开始,基本数学运算得到了加速。
  • @Dunes 实际上我已经在使用这种方法,因为我正在计算整个栅格的标准差。但是,我使用的是数据集的 GDAL 驱动程序建议的自然块大小。这通常是一次一个图像行。您自己实现标准差函数的方式似乎是可行的方法,因为我可以通过使用常量预分配缓冲区(在您的示例中为tmp)来避免每个块的额外分配。
  • @moarningsun 感谢您的提示。我假设 NumPy 至少会使用 SSE2 的 SIMD 指令更长的时间。尤其是标准差(减法、乘法和求和)的计算应该得到很好的加速。
【解决方案2】:

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

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

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

所以速度提高了大约 300 倍。仍然不如numpy版本..

【讨论】:

  • 现在您的测试数据和临时数据仍然适合 CPU 缓存。我认为随着输入数组变大,你会看到 Numpy 相对于 Cython 函数的性能更差。
  • 当我尝试上面实现的 welford 算法时,我得到的结果与 np.std 的结果大不相同。给定arr = np.arange(10, dtype=float),然后是std_welford(arr) == 3.02...,但是np.std(arr) == 2.87...
  • @Dunes,只是定义问题。试试np.std(arr, ddof=1)
  • @moarningsun 谢谢我知道我错过了一些东西。我的统计数据有点生疏了。
  • @kezzos 我已经用测试结果更新了这个问题。它表明存在对数组大小的依赖性。使用较小的数组,您的实现会更快,而使用较大的数组 NumPy 会更快,尽管所有 NumPy 函数一起迭代数组多次。但是 SIMD 指令集带来的加速似乎超过了这个因素。
猜你喜欢
  • 2021-11-16
  • 1970-01-01
  • 2021-05-05
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2022-11-01
  • 1970-01-01
相关资源
最近更新 更多