I am trying to use the one-sided Simultaneous Perturbation (SP) method to get a gradient approximation for multi-variable function.

The equations are very simple: $g_{i}=\frac{f(x+c\delta)-f(x)}{c{\delta}_i}$. which gives the $i^{th}$ component of the gradient. $\delta$ vector is a Bernoulli distribution ±1 with probability of 0.5 for each ±1 outcome. The step size c is fixed for an iteration of the approach. In the original paper, Spall uses this as one step in the approximate calculation of a minima, hence there are more that one iteration. There is just one step here (?).

What bothers me is that, in the above equation the magnitude of each gradient component is the same, with $\delta_i$ causing a possible difference in sign. Have any you used this method? Or, would you know of a better gradient approximation method? I really want to avoid 2N or even N function evaluations for a N variable function.