Simultaneous-perturbation optimization

What it fills

• Minimise the objective over the parameters

  Search the parameters of a fixed circuit family for the ones that minimise a chosen objective, given that every evaluation of that objective is a noisy estimate someone paid shots for. The output is a routine that prepares one particular state — the family collapsed to a member.

When it applies

Spall states the saving as a count of measurements, against the finite-difference alternative: "In contrast to SA algorithms based on finite difference methods, which require 2p (noisy) measurements of L at each iteration, the 'simultaneous perturbation' algorithm here requires only 2q, q ≥ 1, measurements of L at each iteration, where for large p we typically have q ≪ p" (§I). The mechanism is named in the same paper: "this estimate differs from the usual finite difference approximation in that only two measurements (instead of 2p) are used. (The name 'simultaneous perturbation'... arises from the fact that all elements of the θ_k vector are being varied simultaneously.)" (§II). **This citation is a 1992 control-theory paper and predates variational quantum algorithms entirely.** It is cited here for what it contains — the optimizer and its measurement count — and for nothing about quantum circuits; the claim this node makes is a claim about SPSA.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  loss function  L : R^p -> R^1, differentiable, root theta* of
       g(theta) = dL(theta)/dtheta sought (the general multivariate SA
       setting -- L is any differentiable objective, nothing here is
       specific to quantum circuits or to variational algorithms)
       only NOISY measurements of L are available; L and g are not
       observed directly
       starting point  theta_hat_0
       gain sequences  a_k, c_k > 0  with  a_k -> 0, c_k -> 0  as k -> infinity,
           sum_k a_k = infinity,  sum_k (a_k / c_k)^2 < infinity            (A1)
       # the standard form used below for the asymptotic-normality result is
       #   a_k = a / k^alpha,  c_k = c / k^gamma,  a, c, alpha, gamma > 0
       #   -- this paper does NOT use a stabilizing "+A" offset in either
       #   sequence; that form is from later SPSA literature, not here

requires  each perturbation draw Delta_k = (Delta_k1, ..., Delta_kp) has its
    p components mutually independent, mean-zero, i.i.d., symmetric about 0,
    almost-surely bounded (|Delta_ki| <= alpha_0), and with a BOUNDED INVERSE
    MOMENT  E|Delta_ki^-1| <= alpha_1                                (Lemma 1)
    # this inverse-moment condition is the one the bias bound below needs,
    #   and it RULES OUT the uniform and the normal distribution for
    #   Delta_ki; symmetric Bernoulli +-1 satisfies it and is what the
    #   paper's own numerical study draws                     (Sec. II; V-A)
    L has continuous third partial derivatives, uniformly bounded near each
    theta_hat_k, on the relevant domain                              (Lemma 1)

for k = 0, 1, 2, ...
    # --- one gradient estimate from TWO measurements, at any dimension p ---
    draw  Delta_k  satisfying the condition above, independent of
        theta_hat_0, ..., theta_hat_k
    y_k(+) = L(theta_hat_k + c_k Delta_k) + eps_k(+)
    y_k(-) = L(theta_hat_k - c_k Delta_k) + eps_k(-)
        # eps_k(+), eps_k(-) are measurement noise; only their DIFFERENCE
        #   need satisfy  E[eps_k(+) - eps_k(-) | F_k, Delta_k] = 0 a.s.
        #   -- individually biased or dependent noise is not ruled out (Sec. II)
    for l = 1, ..., p
        ghat_kl(theta_hat_k) = (y_k(+) - y_k(-)) / (2 c_k Delta_kl)     (Eq. 2.2)
    # optionally average q >= 1 independent (Delta_k, y+, y-) triples drawn
    #   at this same theta_hat_k before taking the step:
    #   ghat_k = (1/q) sum_{j=1..q} ghat_k^(j)                          (Eq. 2.3)
    #   -- this costs 2q measurements of L this iteration, against 2p for
    #   the finite-difference (FDSA) alternative; for large p, q << p (Sec. I)

    # --- Robbins-Monro step -------------------------------------------------
    theta_hat_{k+1} = theta_hat_k - a_k * ghat_k(theta_hat_k)           (Eq. 2.1)
        # ghat_k is a BIASED estimator of g(theta_hat_k): under the two
        #   conditions above, the bias is b_k(theta_hat_k) = O(c_k^2) as
        #   c_k -> 0 -- proved, not assumed away                (Lemma 1, Eq. 3.2)

return  theta_hat_k  as k -> infinity
# under A1 above plus regularity conditions A2-A5 (boundedness of the
#   iterate, a well-behaved limiting ODE), theta_hat_k -> theta* almost
#   surely -- convergence to a root, not a proof of global optimality (Prop. 1)
# with the standard gains a_k = a/k^alpha, c_k = c/k^gamma and
#   beta = alpha - 2*gamma > 0, 3*gamma - alpha/2 >= 0, the SCALED error
#   k^(beta/2) * (theta_hat_k - theta*) converges IN DISTRIBUTION to
#   N(mu, P M P^T) -- an asymptotic statement about the tail of the
#   sequence, not a bound at any finite k                          (Prop. 2)
#   -- Prop. 2 is STATED FOR q = 1; for q >= 2 the same normality holds with
#   a, c held constant and sigma^2 in M replaced by sigma^2/q, invoking the
#   strengthened A2'                                (Note after Prop. 2, p. 337)

# this paper proves no bound on the number of iterations needed and no
#   shot/measurement cost in a quantum sense; the only proved measurement
#   saving is per-iteration -- 2q here against 2p for finite-difference SA
# the asymptotic measurement-ratio conclusion is CONDITIONAL, not general:
#   (4.1) compares SPSA and FDSA at equal measurement count only under
#   sigma^2 = sigma_tilde^2 and the uniform-integrability caveat of the note
#   after Prop. 2; the reading that FDSA needs at least p/q times more
#   measurements -- p times more when mu = 0 -- is (4.2) under the further
#   assumptions mu approx mu_tilde, M = M_tilde, rho^2 = xi^2 = 1, common a
#   and c                                                    (Sec. IV, Eq. 4.2)
#   Spall says outright it "does not seem possible to make more general
#   statements about the relative superiority of SPSA or FDSA", because
#   rho^2 and xi^2 depend on the Delta_ki distribution and the mixed third
#   partials of L may or may not reduce SPSA's mu^T mu term       (Sec. IV)
#   and footnote 6 gives a normalization -- |Delta_kl| = 1/sqrt(p), equalizing
#   the two difference intervals instead of setting rho^2 = 1 -- under which
#   the mu = 0 ratio becomes p^(1-beta) > 1, i.e. FDSA the cheaper of the two
#                                                        (Sec. IV, footnote 6)
# a single p=20 numerical example found SPSA reaching FDSA's terminal
#   accuracy on about 1/20 to 1/4 of the measurements, LARGER than the
#   1/20-to-under-1/1000 range (4.2) predicts asymptotically for that same
#   case -- an OBSERVED finite-sample result, not a proved one       (Sec. V-B)

Cost, as the source states it

Spall (1992): two noisy evaluations of the objective $L$ per iteration whatever the parameter dimension $p$ — $2q$, $q \ge 1$, if $q$ estimates are averaged (Eq. (2.3)) — against $2p$ per iteration for central-difference FDSA (§I–II). These count evaluations of a classical loss, not shots; Spall bounds no shot cost. Nothing bounds the iteration count: the estimate is $O(c_k^2)$-biased at finite $k$ (Lemma 1) and the guarantee is asymptotic, $O_p(k^{-\beta/2})$, $\beta = \alpha - 2\gamma$, under Propositions 1–2 (A2' strengthened, stated for $q = 1$). By Eq. (4.2) FDSA needs at least $p/q$ times more measurements for equal asymptotic MSE, $p$ times when $\mu = 0$ — but only assuming $\mu = \tilde\mu$, $M = \tilde M$, common $a$, $c$, and $\rho^2 = \xi^2 = 1$; under the normalisation $|\Delta_{kl}| = 1/\sqrt{p}$, footnote 6 makes the $\mu = 0$ ratio $p^{1-\beta} > 1$, SPSA worse. Measured on one $p = 20$ problem, SPSA reached FDSA's terminal accuracy on $1/20$ to $1/4$ of the measurements — larger, Spall notes, than (4.2)'s asymptotic $1/20$ to under $1/1000$ (§V-B), not proved.

Implementations

• Spall's own MATLAB listings on the SPSA website

  James C. Spall keeps an SPSA site at the Johns Hopkins University Applied Physics Laboratory, and its MATLAB page — headed "MATLAB Code", at https://www.jhuapl.edu/SPSA/Pages/MATLAB.htm — offers exactly two things: "Matlab—Basic SPSA Algorithm" and "Matlab—Automatic Gain Selection". Each is a slide PDF carrying exactly one embedded link, to a plain `.txt`, but only the first prints a listing: `Matlab-SPSA_Alg.pdf` runs to two pages and puts the code on the second, under "Matlab Code (cont'd) Click TXT File button for editable text.", while `Matlab-Auto_gain_sel-b-w.pdf` is a single page of four bullets that shows no code at all and closes "Click TXT File button for Matlab code." The first states its own scope, "The code below implements “Basic” SPSA for iterations k = 1,2,...,n", with "Global declarations and initialization for program variables theta, alpha, etc. not shown", the perturbation "elements are generated by Bernoulli ±1", and "Program calls an external function “loss” to obtain y(θ) values (noisy loss measurements)". What is *not* published is named as well: "Code for the adaptive SPSA algorithm (Spall (2000) sometimes called “2nd – order SPSA”) is available upon request from James Spall." This is the site Qiskit's own SPSA docstring sends its readers to — "Many examples are presented at the SPSA Web site".

  The basic listing is a single loop of nine lines, the two-measurement estimate and the Robbins-Monro step written out: `ak=a/(k+1+A)^alpha;` `ck=c/(k+1)^gamma;` `delta=2*round(rand(p,1))-1;` `thetaplus=theta+ck*delta;` `thetaminus=theta-ck*delta;` `yplus=loss(thetaplus);` `yminus=loss(thetaminus);` `ghat=(yplus-yminus)./(2*ck*delta);` `theta=theta-ak*ghat;`. Two calls to `loss` per pass whatever `p` is, one perturbation vector per iteration with no averaging of several gradient estimates, and an elementwise division `./` so that one scalar difference becomes all `p` components. The learning-rate gain carries a stability constant `A`, and the companion PDF attributes the guidelines for choosing `a`, `A` and `c` to a later paper rather than to the 1992 one this method records: "This code uses the guidelines in Spall, IEEE Transactions on Aerospace and Electronic Systems, 1998, pp. 817–823 to pick a, A, and c in gains of form ak = a/(k+1+A).602, ck = c/(k+1).101, k ≥ 0, for the basic SPSA algorithm". The code that bullet points at, `gainsSPSA.txt`, asks the user for five numbers at the prompt and returns the three constants: `A` is `.10` times the expected number of loss evaluations per run, divided by twice the number of averaged gradients; `c = max(c/gavg^0.5, .0001)` starting from "the standard deviation of the measurement noise at i.c."; and `a=step*((A+1)^alpha)/meangbar`, the desired first-step magnitude divided by the magnitude of the gradient estimate averaged over `NL/(2*gavg)` estimates. Constraints are offered only as an optional two-line addition below the update in the basic listing, `theta=min(theta,thetamax);` and `theta=max(theta,thetamin);`, and the gain-selection slide's own caveat about the values it produces is that "Refinement of these values based on numerical experimentation may be helpful".

  Nothing — neither listing carries a loss function or a problem. The basic code says so itself, "Algorithm initialization not shown above (see discussion in introduction to MATLAB code).", and the gain-selection code ships four placeholder assignments in their place, `sigma=TBD;`, `p=TBD;`, `loss='loss_TBD';` and `theta=TBD;`, so the noise level, the dimension, the objective and the starting point are all the reader's to supply. The numbers the files do fix are the exponent `alpha=0.602` (with the comment that it "is often used as a companion to gamma = 0.101"), the `.0001` floor on `c`, the `.10` factor in `A`, and the two seeds the gain-selection code sets for reproducibility, `rand('seed',31415927)` for the perturbations and `randn('seed',111113)` for the noise in the loss measurements. No run, no instance and no output is reported on the page.

  Two MATLAB fragments, both reached from https://www.jhuapl.edu/SPSA/Pages/MATLAB.htm. The basic algorithm is a two-page slide PDF at https://www.jhuapl.edu/spsa/PDF-SPSA/Matlab-SPSA_Alg.pdf whose editable text is https://www.jhuapl.edu/spsa/PDF-SPSA/MatlabCode-TXT.txt; the gain selection is the one-page https://www.jhuapl.edu/spsa/PDF-SPSA/Matlab-Auto_gain_sel-b-w.pdf with its code at https://www.jhuapl.edu/spsa/PDF-SPSA/gainsSPSA.txt. Both `.txt` files are fragments rather than functions — no `function` line, no arguments, and `global sigma p z` at the top of the gain-selection one — and neither carries a licence, a version number or a test, nor declares a dependency beyond MATLAB itself and the external `loss` routine the caller must write. There is no repository: the PDFs' own metadata records creation in November 2000 with a last modification in July 2004 for the basic listing, and December 2000 with a last modification in April 2001 for the gain selection, and the page offers a postal address and a `mailto:` for the author instead, "James C. Spall, Johns Hopkins University Applied Physics Laboratory, 11100 Johns Hopkins Rd., Laurel, MD 20723-6099 USA".

• The SPSA optimizer in qiskit-algorithms

  `SPSA` is shipped as a general-purpose optimizer class by `qiskit-algorithms`, a community package in the `qiskit-community` organisation installed with `pip install qiskit-algorithms` and not part of the Qiskit SDK, whose README carries the warning "**Qiskit Algorithms is no longer officially supported by IBM**. Like any other Apache 2 licensed code, you are free to use it or/and extend it, but please be aware that it is under your own risk." The module docstring is a title line and, after a blank line, one sentence — "Simultaneous Perturbation Stochastic Approximation (SPSA) optimizer." and "This implementation allows both standard first-order and second-order SPSA." — and the class docstring names the property this method exists for: "The main feature of SPSA is the stochastic gradient approximation, which requires only two measurements of the objective function, regardless of the dimension of the optimization problem." It also carries a recommendation, in the library's own voice and about quantum work specifically: "If you are executing a variational algorithm using a Quantum ASseMbly Language (QASM) simulator or a real device, SPSA would be the most recommended choice among the optimizers provided here." For the algorithm itself the docstring cites not the 1992 paper but Spall's 1998 overview, and it sends the reader on to Spall's own site — "Many examples are presented at the SPSA Web site".

  The estimate is the two-point difference along a Bernoulli $\pm 1$ direction: `bernoulli_perturbation` returns `1 - 2 * algorithm_globals.random.binomial(1, 0.5, size=dim)`, `_point_sample` evaluates the loss at `x + eps * delta1` and `x - eps * delta1` and charges 2 to the evaluation counter, and the gains are power series, `powerseries(eta, power, offset)` yielding `eta / ((n + offset) ** power)`. Around that core the class adds machinery the bare method does not have. `blocking` re-evaluates the loss at the proposed point and rejects the step unless it improves — `if fx + self.allowed_increase <= fx_next: ... continue` — at the price of a third evaluation per iteration, and with `allowed_increase=None` the tolerance is itself calibrated, from 25 further evaluations, to "twice the approximated standard deviation of the loss function". `trust_region=True` renormalises any update whose norm exceeds 1. `resamplings` averages several independent directions at the same point and accepts a dict read as `{iteration: number of resamplings per iteration}`. `perturbation_dims` narrows the perturbation to a fixed number of coordinates "chosen uniformly at random", so perturbing every parameter at once is a default here rather than a fixed property. `second_order=True` selects 2-SPSA: the Hessian is estimated from two more evaluations per sample, smoothed as `k / (k + 1) * self._smoothed_hessian + 1 / (k + 1) * hessian`, forced symmetric positive definite with a `regularization` that the constructor sets to 0.01 when it is left `None`, and applied by solving `lse_solver(spd_hessian, gradient)` — deferrable by `hessian_delay` iterations. `last_avg > 1` returns the mean of the last few iterates instead of the last one. Two arguments in `minimize`'s own signature are dead: `jac` and `bounds` appear nowhere in the body, and `get_support_level` declares both `"gradient"` and `"bounds"` as `OptimizerSupportLevel.ignored`.

  Toy circuits, in simulation — there is no dataset and no hardware anywhere in the module or its test. The docstring example minimises the expectation of `SparsePauliOp("ZZ")` on `pauli_two_design(2, reps=1, seed=2)` through a `StatevectorEstimator`, from a random initial point, with `SPSA(maxiter=300)`. The repository's own `test/optimizers/test_spsa.py` runs the same shape one qubit larger: `pauli_two_design(3, reps=1, seed=1)`, observable `SparsePauliOp("ZZI")`, the loss read as an exact `Statevector(bound_circ).expectation_value(obs).real` rather than sampled, a hard-coded six-element initial point, `algorithm_globals.random_seed = 12`, and settings `{"maxiter": 100, "blocking": True, "allowed_increase": 0}`. That file holds ten test methods in all. Of the nine besides `test_pauli_two_design`, seven use analytic scalars — `-(x**2)`, `(np.linalg.norm(x) - 2) ** 2` in three of them, `np.linalg.norm(x) + self.random_generator.random(1)`, `sum(x)`, and an objective that is the identity, `return x` — while the other two are circuit tests for the `QNSPSA` subclass rather than for `SPSA`: `test_qnspsa_fidelity_primitives` on `pauli_two_design(2, reps=1, seed=2)`, and `test_qnspsa_max_evals_grouped` on `pauli_two_design(3, reps=1, seed=1)` with `SparsePauliOp("ZZI")`, a `StatevectorEstimator(seed=12)` and a `StatevectorSampler(seed=12, default_shots=10_000)`.

  `qiskit_algorithms/optimizers/spsa.py` in https://github.com/qiskit-community/qiskit-algorithms — Python, Apache License 2.0, header "(C) Copyright IBM 2018, 2026". It is used as `from qiskit_algorithms.optimizers import SPSA` and then `SPSA(maxiter=300).minimize(loss, x0=initial_point)`; the defaults are `maxiter=100`, `blocking=False`, `trust_region=False`, `resamplings=1`, `last_avg=1`, `second_order=False`, with `learning_rate` and `perturbation` both `None`, which is what triggers calibration ("Note that either both or none must be set"). Calibration is a static method whose full signature is `calibrate(loss, initial_point, c=0.2, stability_constant=0, target_magnitude=None, alpha=0.602, gamma=0.101, modelspace=False, max_evals_grouped=1)`: it spends `steps = 25` two-point samples — 50 loss evaluations — and averages `|delta / (2 * c)|` over them. What it does with that average depends on the parameter `minimize` never passes: on the default branch `a = target_magnitude / avg_magnitudes`, and under `modelspace=True`, `a = target_magnitude / (avg_magnitudes**2)` — "modelspace: Whether the target magnitude is the difference of parameter values or function values (= model space)". `target_magnitude` defaults to `2 * np.pi / 10`, and "Calibration failed, using {target_magnitude} for `a`" is warned when `a < 1e-10`. The two exponents are the ones on Spall's own gain-selection slide and its `.txt`, `alpha=0.602` and `gamma=0.101`, but the recipe for `a` is credited elsewhere: "For further details on the automatic calibration, please refer to the supplementary information section IV. of [3]", where [3] is Kandala et al. (2017). **Three blocks of loss evaluations sit outside the number the optimizer reports.** `minimize` calls `calibrate` and only afterwards sets `self._nfev = 0`, so calibration's 50 evaluations never reach the counter; then, when `blocking` is on with `allowed_increase` left `None`, `self.allowed_increase = 2 * self.estimate_stddev(fun, x, ...)` runs `_batch_evaluate(loss, avg * [initial_point], ...)` with `avg` defaulting to 25 and never touches `_nfev` either; and the closing `result.fun = fun(x)` does not increment it. The first two are conditional and the last is not, which is why the repository's own test — passing `allowed_increase: 0`, so the second block never runs — still gets its `maxiter * 3 + 1` count exactly while the 50 calibration evaluations it also spends stay invisible. The quantum-natural variant is a separate `class QNSPSA(SPSA)` in the sibling file `qiskit_algorithms/optimizers/qnspsa.py`, and is not this method. The module is in the released package: `qiskit-algorithms` 0.4.0 on PyPI, uploaded 2025-08-29, Apache-2.0, Python >= 3.9, and the file at tag `0.4.0` differs from `main` only in the copyright year, one blank line, and whether `Callable` is imported from `typing` or `collections.abc`.

  The module reports no accuracy figures of its own; the numbers attached to it are its test's assertions, and they are about cost as much as accuracy. `test_pauli_two_design` runs the three-qubit case for 100 iterations and asserts `result.fun < -0.95` — the observable is a Pauli string, so the floor is $-1$ — and then names the evaluation count exactly: `maxiter * 3 + 1` for plain SPSA, `maxiter * 5 + 1` with `second_order=True`, `maxiter * 7 + 1` for `QNSPSA`. Those are the per-iteration prices the implementation's own arithmetic predicts — two evaluations for the gradient and one for `blocking`, four and one when the Hessian is estimated as well — with the single evaluation `blocking` makes before the loop as the trailing `+ 1`. The loss is exact in all three arms, but the test is not free of shot noise: the `qnspsa` arm passes `settings["fidelity"] = QNSPSA.get_fidelity(circuit, sampler=StatevectorSampler(seed=123))`, and `StatevectorSampler.__init__` is `def __init__(self, *, default_shots: int = 1024, seed=...)`, so that arm's fidelity is sampled at 1024 shots. Nothing here runs on hardware, and no assertion compares the result against an optimizer from outside the SPSA family.

What it needs

Nobody has taken this apart yet. That is a gap in this graph, not a claim that the method has no parts.

Other ways to fill the same slot

Different approaches

• Conditional-value-at-risk objective

  Change what the classical loop is minimising rather than how it minimises. Instead of averaging every measurement outcome into an expectation value, keep only the best tail of them and average that — which is defensible precisely when the answer is a single good bitstring rather than a physical average.

• Analytic-gradient parameter search

  Get the gradient of the objective exactly, rather than by finite differences, by running the same circuit again at shifted parameter values. The direction is then not an estimate of a slope taken from two noisy numbers; it is the slope, estimated to whatever precision the shots allow.

• Minimise the energy variance

  Minimise how much the energy fluctuates rather than the energy itself. Any eigenstate has zero variance, so the objective's own value tells you whether you have arrived — which the energy never does, since a low number is only low relative to a minimum nobody knows.

• Follow the steepest descent in the state's own geometry

  Take the step that moves the state fastest, not the one that moves the parameters fastest. The two differ because equal changes in parameters do not make equal changes in the state, and the metric measuring that difference has to be estimated before every step.

• Grow the circuit a layer at a time while training it

  Do not settle the circuit before optimising it. Start shallow, train what is there, then hold most of it fixed and add the next layer on top — so every step of the search runs on a shallow circuit with few free parameters, which is where a gradient is still large enough to follow.

In the Atlas

• SPSA-optimized VQE

  Simultaneous perturbation estimates a stochastic gradient with two objective evaluations per iteration.

Sources

• Multivariate stochastic approximation using a simultaneous perturbation gradient approximation↗J. C. Spall · 1992theory yes · simulation yes · hardware no