Simultaneous perturbation stochastic approximation
Simultaneous perturbation stochastic approximation (SPSA) is an algorithmic method for optimizing systems with multiple unknown parameters. It is a type of stochastic approximation algorithm. As an optimization method, it is appropriately suited to large-scale population models, adaptive modeling, simulation optimization, and atmospheric modeling. Many examples are presented at the SPSA website http://www.jhuapl.edu/SPSA. A comprehensive recent book on the subject is Bhatnagar et al. (2013). An early paper on the subject is Spall (1987) and the foundational paper providing the key theory and justification is Spall (1992). SPSA is a descent method capable of finding global  minima, sharing this property with other methods as simulated annealing. Its main feature is the gradient approximation that requires only two measurements of the objective function, regardless of the dimension of the optimization problem. Recall that we want to find the optimal control formula_1 with loss function formula_2: Both Finite Differences Stochastic Approximation (FDSA) and SPSA use the same iterative process: where formula_5 represents the formula_6 iterate, formula_7 is the estimate of the gradient of the objective function formula_8 evaluated at formula_9, and formula_10 is a positive number sequence converging to 0. If formula_11 is a "p"-dimensional vector, the formula_12  component of the symmetric finite difference gradient estimator is: "1 ≤i ≤p", where formula_14 is the unit vector with a 1 in the formula_12 place, and formula_16is a small positive number that decreases with "n". With this method, "2p" evaluations of "J" for each formula_17 are needed. Clearly, when "p" is large, this estimator loses efficiency. Let now formula_18 be a random perturbation vector. The formula_12 component of the stochastic perturbation gradient estimator is: Remark that FD perturbs only one direction at a time, while the SP estimator disturbs all directions at the same time (the numerator is identical in  all "p" components). The number of loss function measurements needed in the SPSA method for each formula_17 is always 2, independent of the dimension "p". Thus, SPSA uses "p" times fewer function evaluations than FDSA, which makes it a lot more efficient. Simple experiments with "p=2" showed that SPSA converges in the same number of iterations as FDSA. The latter follows approximately the steepest descent direction, behaving like the gradient method. On the other hand, SPSA, with the random search direction, does not follow exactly the gradient path. In average though, it tracks it nearly because the gradient approximation is  an almost unbiased estimator of the gradient, as shown in the following lemma. Convergence lemma. Denote by the bias in the estimator formula_23. Assume that formula_24 are all mutually independent with zero-mean, bounded second moments, and formula_25 uniformly bounded. Then formula_26→0 w.p. 1. Sketch of the proof. The main idea is to use conditioning on formula_18 to express formula_28 and then to use a second order Taylor expansion of formula_29 and formula_30. After algebraic manipulations using the zero mean and the independence of formula_24, we get The result follows from the hypothesis that formula_16→0. Next we resume some of  the hypotheses under which formula_34 converges in probability to the set of global minima of formula_2. The efficiency of the method depends on the shape of formula_2, the values of the parameters formula_37 and formula_16 and the distribution of the perturbation terms formula_39. First, the algorithm parameters must satisfy the following conditions: A good choice for formula_39 is the Rademacher distribution, i.e. Bernoulli +-1 with probability 0.5. Other choices are possible too, but note that the uniform and normal distributions cannot be used because they do not satisfy the finite inverse moment conditions. The loss function "J(u)" must  be thrice continuously differentiable and the individual elements of the third derivative must be bounded: formula_51. Also, formula_52 as formula_53. In addition, formula_54 must be Lipschitz continuous, bounded and the ODE formula_55 must have a unique solution for each initial condition. Under these conditions and a few others, formula_56 converges in probability to the set of global minima of J(u) (see Maryak and Chin, 2008). It has been shown that differentiability is not required: continuity and convexity are sufficient for convergence. Extension to second-order (Newton) methods. It is known that a stochastic version of the standard (deterministic) Newton-Raphson  algorithm (a “second-order” method) provides an asymptotically optimal or near-optimal form of stochastic approximation. SPSA can also be used to efficiently estimate the Hessian matrix of the loss function based on either noisy loss measurements or noisy gradient measurements (stochastic gradients). As with the basic SPSA method, only a small fixed number of loss measurements or gradient measurements are needed at each iteration, regardless of the problem dimension "p". See the brief discussion in Stochastic gradient descent.