SPSA: A Stochastic Optimization Method

• SPSA is a stochastic optimization approach that uses random Bernoulli perturbations to estimate gradients with only two function evaluations per iteration.
• It employs decaying gain schedules to balance bias and variance, achieving asymptotically unbiased estimates even in noisy, high-dimensional settings.
• Variants like one-measurement, second-order, and parallel SPSA extend its applicability to constrained, distributed, and simulation-based optimization with strong convergence guarantees.

Simultaneous Perturbation Stochastic Approximation (SPSA) is a stochastic optimization methodology that estimates gradients and, in extended forms, higher-order derivatives, using only a minimal number of function evaluations per iteration. It is uniquely efficient for high-dimensional, noisy, or simulation-based settings where derivative information is unavailable or expensive to compute. SPSA and its variants have substantial theoretical guarantees and a broad empirical track record across control, machine learning, engineering, and operations research.

1. SPSA Fundamentals and Gradient Estimation

SPSA was introduced to minimize objective functions of the form

$J(\theta) = E[f(\theta) + \text{noise}],$

where only noisy observations $y(\theta) = J(\theta) + \epsilon$ are available. The classical stochastic approximation method—finite-difference (Kiefer–Wolfowitz)—requires $2d$ function measurements per iteration in $d$ dimensions; SPSA achieves an unbiased gradient estimate with only two function evaluations per iterate, independent of the dimension.

At iteration $k$, the standard SPSA estimator generates a random perturbation vector $\Delta_k\in\{+1,-1\}^d$ (entries drawn i.i.d. Bernoulli), and forms the gradient estimate: $$\hat{g}_k = \frac{y(\theta_k + c_k \Delta_k) - y(\theta_k - c_k \Delta_k)}{2 c_k}\, \Delta_k^{-1},$$ where $c_k>0$ is a small, decaying perturbation gain, and $\Delta_k^{-1}$ denotes elementwise reciprocals. The parameter is updated via

$\theta_{k+1} = \theta_k - a_k \hat{g}_k,$

where $a_k$ is the step size, typically $a_k = a/(k+1+A)^\alpha$ with $\alpha \in (0.6,1]$ and $A\ge0$ (Wang, 2020, Li et al., 2022).

This approach does not require explicit calculation or observation of partial derivatives. The expected value of the SPSA estimator is asymptotically unbiased with a bias term of order $O(c_k^2)$, and the variance scales as $O(1/c_k^2)$. Convergence to a local optimum is guaranteed under standard smoothness, noise, and gain sequence conditions (Wang, 2020).

2. Algorithmic Variants and Recent Advances

One-Measurement-Per-Iteration SPSA

To further reduce function evaluation cost, one-measurement methods have been developed. The “SPSA1-A” variant introduces a two-step process:

• Step 1 uses the classic two-point difference to generate $\hat{g}_k$ and performs a "half-step."
• Step 2, without additional function evaluation, samples a Bernoulli direction $\hat{\xi}_k$ predicted to be a descent direction based on $\hat{g}_k$, and performs a second move.

The expected number of function measurements per iteration approaches one asymptotically, with convergence and normality guarantees matching standard SPSA. Experimental results on standard test problems demonstrate that SPSA1-A halves the evaluation cost to reach fixed accuracy thresholds compared to regular SPSA (Li et al., 2022).

Second-Order and Newton-Type SPSA

Second-order SPSA (2SPSA) simultaneously estimates both the gradient and Hessian using $O(1)$ function evaluations per iterate, supporting faster local convergence:

• The Hessian is estimated by combining two simultaneous perturbations and symmetrizing the finite-difference matrix.
• A quasi-Newton update is performed with the Hessian inverse.

Recent work has reduced the computational burden of inverting these matrices from $O(p^3)$ to $O(p^2)$ per iteration by maintaining an $L\,B\,L^{\top}$ symmetric indefinite factorization, improving practical applicability in high dimensions without compromising formal convergence rates (Zhu et al., 2019).

Parallel and Distributed SPSA

The Parallel SPSA (PSPO) variant exploits multi-core or distributed architectures by evaluating multiple random perturbations in parallel, forming a least-squares problem to fit the gradient. This achieves significant wall-clock reduction in high-noise and high-dimensional simulation optimization problems (Alaeddini et al., 2017, Alaeddini et al., 2017).

Distributed or decentralized adaptations (DSPG) allow each agent in a network to carry out SPSA-type updates locally, exchanging information stochastically over unreliable networks. Convergence to an $O(c)$ neighborhood of the global minimizer is provable despite communication delays and asynchrony (Ramaswamy, 2019).

Discrete and Constrained SPSA

SPSA has been extended to discrete settings (DSPSA), where updates are performed in a relaxed continuous space and projected or randomized to a discrete set. Convergence and rate of decay of the probability of not attaining the optimal solution have been established, with performance comparable to stochastic ruler and comparison algorithms (Wang, 2013, Mandal et al., 2023).

For inequality-constrained problems, "switch updating" SPSA combines classical unconstrained updates with targeted, one-measurement feasibility-restoring steps along violated constraint gradients, ensuring final feasibility and strong convergence properties without resorting to costly projections or sensitive penalty tuning (Jia et al., 2023).

For equality constraints, derivative-free SQP (DF-SSQP) methods employ SPSA-based gradient and Hessian surrogates for both objective and constraints, with momentum-style aggregation to manage the bias-variance trade-off, achieving dimension-independent evaluation complexity and classical local inference rates (Na, 25 Oct 2025).

Modifications for Variance, Tuning, and Stability

Optimal choice of perturbation vectors is Bernoulli in the asymptotic regime; for small-sample problems, specially designed segmented uniform distributions may outperform Bernoulli under certain tuning regimes (Cao, 2014).

Variants such as state-dependent gains and negatively correlated ("zig-zag") exploration have been shown to improve global stability and variance properties for the one-sided SPSA estimator (1SPSA), with strong theoretical justifications for use in unconstrained and ill-conditioned problems (Lauand et al., 4 Sep 2025).

3. Implementation, Tuning, and Theoretical Guarantees

Gains and Perturbations

The canonical gain schedules are: $$a_k = \frac{a}{(k+1+A)^{\alpha}}, \qquad c_k = \frac{c}{(k+1)^{\gamma}},$$ with $\alpha \in (0.6,1.0]$, $\gamma \approx 0.1$–$0.2$, and $A$ a stabilizing offset. Bernoulli $\pm1$ perturbations are optimal asymptotically, but the segmented uniform strategy may be preferred for short runs (Wang, 2020, Cao, 2014). Larger $a$ can speed initial convergence but amplify noise; the decay of $c_k$ trades bias versus variance in the gradient estimate (Li et al., 2022).

Convergence

Under standard assumptions (three continuous derivatives of $J$, bounded noise, compactness of iterates, step-size properties), strong convergence and asymptotic normality are guaranteed. SPSA converges almost surely to local minima, with precise scaling of the asymptotic variance, mean shift, and conditions for asymptotic normality (with explicit formulas for the limiting covariance matrix) (Li et al., 2022, Wang, 2020).

In constrained settings, similar ODE-based arguments and momentum-style bias correction ensure global and local optimality, feasibility preservation, and statistical inference for solution accuracy (Jia et al., 2023, Na, 25 Oct 2025).

4. Applications and Empirical Performance

SPSA and its extensions are used in:

• Reinforcement learning and policy search, where gradient information is inherently unavailable and simulation cost is high. SPSA's zeroth-order nature enables efficient policy tuning even with thousands of parameters (Wang, 2020).
• Hyperparameter and weight optimization in deep learning and meta-learning, where zero-order SPSA-based tracking often outperforms gradient-based meta-updates on few-shot classification (e.g., MAML, ProtoNet) (Boiarov et al., 2021).
• Stochastic simulation-based calibration, including epidemiological models (e.g., epidemic fitting via PSPO), satellite beamforming design (successive sub-array selection), and adaptive labor staffing under high-dimensional Markov cost models (Alaeddini et al., 2017, Alaeddini et al., 2017, Chen et al., 2021, Prashanth et al., 2013).
• Large-scale Markov-chain optimization (SM-SPSA), including web-graph ranking, where the coordinate-wise parameterization and transformations make high-dimensional problems tractable directly in the probability simplex (Dieleman et al., 20 Jul 2024).

Empirical studies report dramatic reductions in function evaluation requirements and wall-clock time versus standard finite-difference or derivative-based methods. In practical scenarios with noisy or expensive function measurements, SPSA and its variants achieve high-accuracy solutions at a fraction of the cost, often enabling optimization at scales impractical for classical alternatives (Li et al., 2022, Alaeddini et al., 2017, Na, 25 Oct 2025, Prashanth et al., 2013).

5. Practical Guidelines and Comparative Insights

A condensed table of reference parameter choices for classical SPSA and prominent variants follows:

Variant	Step Size $a_k$	Perturb Size $c_k$	Pert. Dist.	Function Evals/Iter
Classic SPSA	$a/(k+A)^{\alpha}$, $\alpha\approx0.602-1$	$c/(k+1)^{\gamma}$, $\gamma\approx0.101-0.2$	Bernoulli $\pm1$	2
SPSA1-A	Same (usually $\alpha=1$)	Same	Same (plus sign-based $\hat\xi_k$)	1 (asymp. avg)
2SPSA	As above	As above	Bernoulli $\pm1$	4 (gradient+Hess)
PSPO	As above	As above	$M$ Bernoulli	$M$ (parallel)
DSPSA	As above	As above	Bernoulli $\pm1$	2
SM-SPSA	Fixed $\epsilon$, $a=0.1$–$0.2$	$1/(i+1)$	Rademacher$^*$	2 (per updates for each mask entry)

$^*$ For SM-SPSA, the perturbation is masked to adjustable entries and mapped via a logistic transformation.

Crucially, all SPSA-based methods rely on precise tuning of $a$, $c$, and, where applicable, momentum or regularization parameters. Empirical performance is robust to moderate miscalibration, but formal convergence rates and variance-minimization require attention to the gain schedules (Li et al., 2022, Wang, 2020, Zhu et al., 2019).

6. Theoretical and Empirical Comparison with Competing Methods

Compared with finite-difference and gradient-based methods:

• SPSA offers dimension-independent function evaluation complexity.
• In high dimensions, finite-difference requires $2d$ evaluations per step versus SPSA's 2.
• Zero-order nature of SPSA enables deployment in non-differentiable, noisy, and simulation-based settings.
• In meta-learning and hyperparameter tuning, SPSA-based tracking demonstrates greater robustness to estimator variance and better empirical accuracy in data-constrained regimes (Boiarov et al., 2021).
• Parallel and second-order variants achieve substantial speedup, approaching the efficiency of analytical-gradient optimizers when function evaluation parallelism or higher-order structure is exploitable (Alaeddini et al., 2017, Zhu et al., 2019).

In constraint-handling, switch-updating and projection methods via SPSA provide formal feasibility guarantees and maintain competitive mean-square convergence rates, outperforming penalty-based alternatives in both accuracy and constraint satisfaction (Jia et al., 2023, Na, 25 Oct 2025).

Taken together, the Simultaneous Perturbation Stochastic Approximation framework constitutes a foundational, broadly applicable methodology in stochastic optimization, with powerful scalability properties, robust theoretical analysis, and a wide range of highly effective algorithmic variants (Wang, 2020, Li et al., 2022, Cao, 2014, Na, 25 Oct 2025).