Simultaneous perturbation stochastic approximation (SPSA) is analgorithmic method for optimizing systems with multiple unknownparameters. It is a type ofstochastic approximation algorithm. As an optimization method, it is appropriately suited to large-scale population models, adaptive modeling, simulationoptimization, andatmospheric modeling. Many examples are presented at the SPSA websitehttp://www.jhuapl.edu/SPSA. A comprehensive 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 such assimulated 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${\displaystyle u^{\ast }}$ with lossfunction${\displaystyle J(u)}$:

${\displaystyle u^{\ast }=\arg \underset{u\in U}{min}J(u).}$

Both Finite Differences Stochastic Approximation (FDSA)and SPSA use the same iterative process:

${\displaystyle u_{n+1}=u_n-a_n{\hat{g}}_n(u_n),}$

where${\displaystyle u_n=((u_n{)}_1,(u_n{)}_2,\dots ,(u_n{)}_p{)}^T}$represents the${\displaystyle n^{th}}$ iterate,${\displaystyle {\hat{g}}_n(u_n)}$ is the estimate of the gradient of the objective function${\displaystyle g(u)=\frac{\mathrm{\partial }}{\mathrm{\partial }u}J(u)}$ evaluated at${\displaystyle u_n}$, and${\displaystyle \{a_n\}}$ is a positive number sequence converging to 0. If${\displaystyle u_n}$ is ap-dimensional vector, the${\displaystyle i^{th}}$ component of thesymmetric finite difference gradient estimator is:

FD:${\displaystyle (\widehat{g_n}(u_n){)}_i=\frac{J(u_n+c_n{e}_i)-J(u_n-c_n{e}_i)}{2c_n},}$

1 ≤i ≤p, where${\displaystyle e_i}$ is theunit vector with a 1 in the${\displaystyle i^{th}}$place, and${\displaystyle c_n}$ is a small positive number that decreases withn. With this method,2p evaluations ofJ for each${\displaystyle g_n}$ are needed. Whenp is large, this estimator loses efficiency.

Let now${\displaystyle {\mathrm{\Delta }}_n}$ be a random perturbation vector. The${\displaystyle i^{th}}$ component of the stochastic perturbation gradient estimator is:

SP:${\displaystyle (\widehat{g_n}(u_n){)}_i=\frac{J(u_n+c_n{\mathrm{\Delta }}_n)-J(u_n-c_n{\mathrm{\Delta }}_n)}{2c_n({\mathrm{\Delta }}_n{)}_i}.}$

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 allp components). The number of loss function measurements needed in the SPSA method for each${\displaystyle g_n}$ is always 2, independent of thedimensionp. Thus, SPSA usesp times fewer function evaluations than FDSA, which makes it a lot more efficient.

Simple experiments withp=2 showed that SPSA converges in the same number of iterations as FDSA. The latter followsapproximately thesteepest 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 almostunbiasedestimator of the gradient, as shown in the following lemma.

Denote by

${\displaystyle b_n=E[{\hat{g}}_n|u_n]-\mathrm{\nabla }J(u_n)}$

the bias in the estimator${\displaystyle {\hat{g}}_n}$. Assume that${\displaystyle \{({\mathrm{\Delta }}_n{)}_i\}}$ are all mutually independent with zero-mean, bounded secondmoments, and${\displaystyle E(|({\mathrm{\Delta }}_n{)}_i{|}^{-1})}$ uniformly bounded. Then${\displaystyle b_n}$→0 w.p. 1.

The mainidea is to use conditioning on${\displaystyle {\mathrm{\Delta }}_n}$ to express${\displaystyle E[({\hat{g}}_n{)}_i]}$ and then to use a second order Taylor expansion of${\displaystyle J(u_n+c_n{\mathrm{\Delta }}_n{)}_i}$ and${\displaystyle J(u_n-c_n{\mathrm{\Delta }}_n{)}_i}$. After algebraic manipulations using the zero mean and the independence of${\displaystyle \{({\mathrm{\Delta }}_n{)}_i\}}$, we get

${\displaystyle E[({\hat{g}}_n{)}_i]=(g_n{)}_i+O(c_n^2)}$

The result follows from thehypothesis that${\displaystyle c_n}$→0.

Next we resume some of the hypotheses under which${\displaystyle u_t}$ converges inprobability to the set of global minima of${\displaystyle J(u)}$. The efficiency ofthe method depends on the shape of${\displaystyle J(u)}$, the values of the parameters${\displaystyle a_n}$ and${\displaystyle c_n}$ and the distribution of the perturbation terms${\displaystyle {\mathrm{\Delta }}_{ni}}$. First, the algorithm parameters must satisfy thefollowing conditions:

• ${\displaystyle a_n}$ >0,${\displaystyle a_n}$→0 when n→∝ and${\displaystyle \sum _{n=1}^{\mathrm{\infty }}{a}_n=\mathrm{\infty }}$. A good choice would be${\displaystyle a_n=\frac{a}{n};}$ a>0;
• ${\displaystyle c_n=\frac{c}{n^{\gamma }}}$, where c>0,${\displaystyle \gamma \in \left[\frac{1}{6},\frac{1}{2}\right]}$;
• 
• ${\displaystyle {\mathrm{\Delta }}_{ni}}$ must be mutually independent zero-mean random variables, symmetrically distributed about zero, with. The inverse first and second moments of the${\displaystyle {\mathrm{\Delta }}_{ni}}$ must be finite.

A good choice for${\displaystyle {\mathrm{\Delta }}_{ni}}$ is theRademacher 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 functionJ(u) must be thrice continuouslydifferentiable and the individual elements of thethird derivative must be bounded:. Also,${\displaystyle |J(u)|\to \mathrm{\infty }}$ as${\displaystyle u\to \mathrm{\infty }}$.

In addition,${\displaystyle \mathrm{\nabla }J}$ must be Lipschitz continuous, bounded and the ODE${\displaystyle \dot{u}=g(u)}$ must have a unique solution for each initial condition.Under these conditions and a few others,${\displaystyle u_k}$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.^[1]

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 theHessian 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 dimensionp. See the brief discussion inStochastic gradient descent.

• Bhatnagar, S., Prasad, H. L., and Prashanth, L. A. (2013),Stochastic Recursive Algorithms for Optimization: Simultaneous Perturbation Methods, Springer[1].
• Hirokami, T., Maeda, Y., Tsukada, H. (2006) "Parameter estimation using simultaneous perturbation stochastic approximation", Electrical Engineering in Japan, 154 (2), 30–3[2]
• Maryak, J.L., and Chin, D.C. (2008), "Global Random Optimization by Simultaneous Perturbation Stochastic Approximation,"IEEE Transactions on Automatic Control, vol. 53, pp. 780-783.
• Spall, J. C. (1987), “A Stochastic Approximation Technique for Generating Maximum Likelihood Parameter Estimates,”Proceedings of the American Control Conference, Minneapolis, MN, June 1987, pp. 1161–1167.
• Spall, J. C. (1992), “Multivariate Stochastic Approximation Using a Simultaneous Perturbation Gradient Approximation,”IEEE Transactions on Automatic Control, vol. 37(3), pp. 332–341.
• Spall, J.C. (1998). "Overview of the Simultaneous Perturbation Method for Efficient Optimization"2.Johns Hopkins APL Technical Digest, 19(4), 482–492.
• Spall, J.C. (2003)Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, Wiley.ISBN 0-471-33052-3 (Chapter 7)

• Numerical climate and weather models
• Stochastic optimization
• Randomized algorithms
• Optimization algorithms and methods

• Articles with short description
• Short description matches Wikidata