title: How to find the local minima of a smooth multidimensional array in NumPy efficiently?
body: <p>Say I have an array in NumPy containing evaluations of a continuous differentiable function, and I want to find the local minima. There is no noise, so every point whose value is lower than the values of all its neighbors meets my criterion for a local minimum.</p>
  
  <p>I have the following list comprehension which works for a two-dimensional array, ignoring potential minima on the boundaries:</p>
  
  <pre><code>import numpy as N
  
  def local_minima(array2d):
      local_minima = [ index 
                       for index in N.ndindex(array2d.shape)
                       if index[0] &gt; 0
                       if index[1] &gt; 0
                       if index[0] &lt; array2d.shape[0] - 1
                       if index[1] &lt; array2d.shape[1] - 1
                       if array2d[index] &lt; array2d[index[0] - 1, index[1] - 1]
                       if array2d[index] &lt; array2d[index[0] - 1, index[1]]
                       if array2d[index] &lt; array2d[index[0] - 1, index[1] + 1]
                       if array2d[index] &lt; array2d[index[0], index[1] - 1]
                       if array2d[index] &lt; array2d[index[0], index[1] + 1]
                       if array2d[index] &lt; array2d[index[0] + 1, index[1] - 1]
                       if array2d[index] &lt; array2d[index[0] + 1, index[1]]
                       if array2d[index] &lt; array2d[index[0] + 1, index[1] + 1]
                     ]
      return local_minima
  </code></pre>
  
  <p>However, this is quite slow. I would also like to get this to work for any number of dimensions. For example, is there an easy way to get all the neighbors of a point in an array of any dimensions? Or am I approaching this problem the wrong way altogether? Should I be using <code>numpy.gradient()</code> instead?</p>