version: "1.0.1" name: spsa description: Concise guide to the SPSA estimator and sampler gradient examples for parameterized quantum circuits.

SPSA Gradient Estimation

Purpose

Use this skill for the SPSA gradient examples in this folder:

• scripts/algorithm.py

They estimate gradients of:

• expectation values via SPSAEstimatorGradient
• sampled output distributions via SPSASamplerGradient

SPSA is useful when circuits have many parameters, because each batch needs only 2 * batch_size circuit evaluations, independent of parameter count.

Overview

For each random perturbation vector delta in {+1, -1}^d, the implementation evaluates:

• f(theta + epsilon * delta)
• f(theta - epsilon * delta)

It then forms a central-difference estimate and averages across batch_size perturbations.

Prerequisites

• NumPy
• Qiskit parameterized circuits
• Estimator or Sampler primitives
• Basic finite-difference gradient intuition

Using the Provided Implementation

The real entry classes are:

• SPSAEstimatorGradient
• SPSASamplerGradient

Public usage goes through inherited run(...); the algorithm itself is implemented in _run(...).

Minimal estimator example:

from qiskit.circuit import QuantumCircuit, ParameterVector
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp
from qiskit_algorithms.gradients import SPSAEstimatorGradient
theta = ParameterVector("theta", 2)
qc = QuantumCircuit(2)
qc.ry(theta[0], 0)
qc.ry(theta[1], 1)
qc.cx(0, 1)
estimator = StatevectorEstimator()
obs = SparsePauliOp("ZZ")
grad = SPSAEstimatorGradient(estimator=estimator, epsilon=0.01, batch_size=4, seed=123)
result = grad.run(
    circuits=[qc],
    observables=[obs],
    parameter_values=[[0.3, -0.2]],
    parameters=[[theta[0], theta[1]]],
    precision=None,
).result()
print(result.gradients)

Minimal sampler example:

from qiskit.circuit import QuantumCircuit, ParameterVector
from qiskit.primitives import StatevectorSampler
from qiskit_algorithms.gradients import SPSASamplerGradient
theta = ParameterVector("theta", 2)
qc = QuantumCircuit(2)
qc.ry(theta[0], 0)
qc.ry(theta[1], 1)
qc.cx(0, 1)
qc.measure_all()
sampler = StatevectorSampler()
grad = SPSASamplerGradient(sampler=sampler, epsilon=0.01, batch_size=4, seed=123)
result = grad.run(
    circuits=[qc],
    parameter_values=[[0.3, -0.2]],
    parameters=[[theta[0], theta[1]]],
    shots=None,
).result()
print(result.gradients)

Core Parameters Explained

SPSAEstimatorGradient

• estimator: expectation-value primitive
• epsilon: perturbation size, must be positive
• batch_size=1: number of SPSA samples to average
• seed=None: RNG seed
• precision=None: estimator precision override
• transpiler=None, transpiler_options=None: optional transpilation

Run inputs:

• circuits
• observables
• parameter_values
• parameters
• precision

Return fields:

• gradients
• metadata
• precision

SPSASamplerGradient

• sampler: sampling primitive
• epsilon: perturbation size, must be positive
• batch_size=1: number of SPSA samples to average
• seed=None: RNG seed
• shots=None: sampler shot override
• transpiler=None, transpiler_options=None: optional transpilation

Run inputs:

• circuits
• parameter_values
• parameters
• shots

Return fields:

• gradients
• metadata
• shots

Implementation Architecture

Estimator path

SPSAEstimatorGradient._run(...) does the following:

1. Normalize precision.
2. Optionally transpile circuits and remap observables with circuit layout.
3. For each circuit, sample batch_size random sign vectors.
4. Build theta + epsilon * delta and theta - epsilon * delta parameter sets.
5. Submit one batched estimator job.
6. Compute

$$ \frac{f_+ - f_-}{2\epsilon} $$ and divide by the perturbation signs to recover per-parameter gradients.

1. Average across the batch and keep only requested parameters.

Simplified reconstruction:

for each circuit:
    deltas = random_pm_one_vectors(batch_size)
    plus = [theta + epsilon * d for d in deltas]
    minus = [theta - epsilon * d for d in deltas]
    pubs.append((circuit, observable, plus + minus, precision))
results = estimator.run(pubs).result()
gradient = mean(((f_plus - f_minus) / (2 * epsilon)) / delta over batch)

Sampler path

SPSASamplerGradient._run(...) is similar, but works on probability distributions:

1. Normalize shots.
2. Optionally transpile circuits.
3. Build plus/minus parameter sets from random sign vectors.
4. Submit one batched sampler job.
5. Convert counts to probabilities.
6. Compute per-bitstring distribution differences.
7. Reconstruct each requested parameter gradient and average over the batch.

Simplified reconstruction:

for each circuit:
    deltas = random_pm_one_vectors(batch_size)
    plus = [theta + epsilon * d for d in deltas]
    minus = [theta - epsilon * d for d in deltas]
    pubs.append((circuit, plus + minus, shots))
results = sampler.run(pubs).result()
for each target parameter j:
    gradient_j = average_k(diff_distribution_k * deltas[k][j])

Understanding the Key Quantum Components

• Parameterized circuit U(theta)
• Estimator objective: expectation value f(theta)
• Sampler objective: output distribution p_theta(z)
• Random perturbation vector delta in {+1, -1}^d
• Quantum evaluations at theta +/- epsilon * delta
• Classical post-processing for gradient recovery

Theory-to-Code Mapping

• theta: parameter_values
• delta: random +1/-1 vectors built with NumPy RNG
• epsilon: constructor argument epsilon
• f(theta + epsilon * delta), f(theta - epsilon * delta): plus and minus evaluations
• requested parameters: parameters
• batch average: NumPy mean or dict-value averaging

Mathematical Deep Dive

For scalar objective f(theta):

$$ \hat g_i(\theta) = \frac{f(\theta + \epsilon \delta) - f(\theta - \epsilon \delta)}{2\epsilon\,\delta_i} $$

For batch_size = B:

$$ \hat g(\theta) = \frac{1}{B}\sum_{k=1}^B \frac{f(\theta + \epsilon \delta^{(k)}) - f(\theta - \epsilon \delta^{(k)})}{2\epsilon} \odot \delta^{(k)} $$

Since delta_i in {+1, -1}, dividing by delta_i is equivalent to multiplying by delta_i.

For sampler gradients, the same idea is applied per output bitstring probability.

Hands-On Example

If you request only a parameter subset, the output keeps that order:

result = grad.run(
    circuits=[qc],
    observables=[obs],
    parameter_values=[[0.1, -0.4, 0.7]],
    parameters=[[theta[0], theta[2]]],
    precision=None,
).result()
print(result.gradients[0])
print(result.metadata[0]["parameters"])

Expected behavior:

• output length matches requested parameter count
• output order matches parameters=[[...]]

Implementing Your Own Version

import numpy as np
def spsa_gradient(eval_fn, theta, epsilon, batch_size, rng):
    grads = []
    for _ in range(batch_size):
        delta = (-1) ** rng.integers(0, 2, len(theta))
        f_plus = eval_fn(theta + epsilon * delta)
        f_minus = eval_fn(theta - epsilon * delta)
        grads.append(((f_plus - f_minus) / (2 * epsilon)) / delta)
    return np.mean(np.asarray(grads), axis=0)

This matches the core logic of the estimator implementation; the sampler version replaces scalar outputs with sparse probability dictionaries.

Debugging Tips

• epsilon <= 0 raises ValueError.
• In estimator mode, the input lists must have matching lengths.
• Requested parameters must exist in the circuit.
• Too-small epsilon can make noise dominate.
• Larger batch_size reduces variance but increases runtime.
• In estimator mode, transpilation changes layout, so observables are remapped.
• In sampler mode, gradients are built from normalized counts, so shot count matters.
• In some Qiskit versions, constructor arguments such as epsilon and batch_size are not exposed as public attributes on gradient objects; prefer tracking configured values in local variables rather than reading grad.epsilon or grad.batch_size.