How to find the local minima of a smooth multidimensional array in NumPy efficiently?

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.

I have the following list comprehension which works for a two-dimensional array, ignoring potential minima on the boundaries:

import numpy as N

def local_minima(array2d):
    local_minima = [ index 
                     for index in N.ndindex(array2d.shape)
                     if index[0] > 0
                     if index[1] > 0
                     if index[0] < array2d.shape[0] - 1
                     if index[1] < array2d.shape[1] - 1
                     if array2d[index] < array2d[index[0] - 1, index[1] - 1]
                     if array2d[index] < array2d[index[0] - 1, index[1]]
                     if array2d[index] < array2d[index[0] - 1, index[1] + 1]
                     if array2d[index] < array2d[index[0], index[1] - 1]
                     if array2d[index] < array2d[index[0], index[1] + 1]
                     if array2d[index] < array2d[index[0] + 1, index[1] - 1]
                     if array2d[index] < array2d[index[0] + 1, index[1]]
                     if array2d[index] < array2d[index[0] + 1, index[1] + 1]
                   ]
    return local_minima

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 numpy.gradient() instead?