title: Memory consumption of NumPy function for standard deviation
body: <p>I'm currently using the Python bindings of GDAL to work on quite large raster data sets (> 4 GB). Since loading them into memory at once is no feasible solution for me I read them into smaller blocks and do the computations piece by piece. To avoid a new allocation for every block read I'm using the <code>buf_obj</code> argument (<a href="http://gdal.org/python/osgeo.gdal_array-module.html#BandReadAsArray" rel="noreferrer">here</a>) to read the values into an preallocated NumPy array. At one point I have to compute the mean and standard deviation of the entire raster. Naturally I've used <code>np.std</code> for the computation. However by profiling the memory consumption of my program I realized that with each invocation of <code>np.std</code> additionally memory is allocated and released.</p>
  
  <p>A minimum working example which demonstrates this behavior:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>A search within the source tree of NumPy on GitHub revealed that the <code>np.std</code> function internally invokes the <code>_var</code> function from <code>_methods.py</code> (<a href="https://github.com/numpy/numpy/blob/master/numpy/core/_methods.py" rel="noreferrer" title="here">here</a>). At one point <code>_var</code> computes the deviations from the mean and sums them up.  Therefore an temporary copy of the input array is created. The function essentially computes the standard deviation as follows:</p>
  
  <pre><code>mu = sum(arr) / len(arr)
  tmp = arr - mu
  tmp = tmp * tmp
  sd = np.sum(tmp) / len(arr)
  </code></pre>
  
  <p>While this approach is OK for smaller input arrays it's definitely no way to go for larger ones. Since I'm using smaller blocks of memory as mentioned before this additional copy is not an game-breaking issue from the memory point of view in my program. However what bugs me is that for each block a new allocation is made and released before reading the next block.</p>
  
  <p>Is there some other function within NumPy or SciPy which utilizes an aproach with constant memory consumption like the Welford algorithm (<a href="https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm" rel="noreferrer">Wikipedia</a>) for one pass computation of the mean and standard deviation?</p>
  
  <p>Another way to go would be to implement a custom version of the <code>_var</code> function with an optional <code>out</code> argument for a preallocated buffer (like NumPy ufuncs). With this approach the additional copy would not be eliminated but at least the memory consumption would be constant and runtime for the allocations in each block is saved.</p>
  
  <p><strong>EDIT:</strong> Tested the Cython implementation of the Welford algorithm as suggested by kezzos.</p>
  
  <p>Cython implementation (modified from kezzos):</p>
  
  <pre><code>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] &lt; a_min:
              a_min = a[i]
          if a[i] &gt; a_max:
              a_max = a[i]
      return a_min, a_max, mean, sqrt(M2 / (n - 1))
  </code></pre>
  
  <p>NumPy implementation (mean and std can possibly be computed in one function):</p>
  
  <pre><code>def vector_approach(a):
      return np.min(a), np.max(a), np.mean(a), np.std(a, ddof=1)
  </code></pre>
  
  <p>Test results using a random data set (times in milliseconds, best of 25):</p>
  
  <pre><code>----------------------------------
  | 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 |
  ----------------------------------
  </code></pre>
  
  <p>It seems like the iterative approach using Cython is faster with smaller data sets and the NumPy vector (possibly SIMD accelerated) approach for larger data sets with 10000+ elements. All tests were conducted with Python 2.7.9 and NumPy version 1.9.2.</p>
  
  <p>Note that in the real case to upper functions would be used to compute the statistics for a single block of the raster. The standard deviations and means for all block are to be combined with methodology suggested in Wikipedia (<a href="https://en.wikipedia.org/wiki/Standard_deviation#Combining_standard_deviations" rel="noreferrer">here</a>). It has the advantage that not all elements of the raster need to be summed up and thereby avoids the float overflow problem (at least to some point).</p>