`spsa` — Spall’s (1992) simultaneous perturbation stochastic approximation¶

What. spsa (T2, ADR-002) is a BLACK-BOX OD calibrator: it never builds a route-choice proportion matrix at all. Each iteration perturbs the current OD estimate along a single random (Rademacher) direction, queries the assignment ORACLE (an inner MSA/AON sweep) at g+ck*Delta and g-ck*Delta, and takes a two-point finite-difference gradient step — the whole method needs only forward simulation queries, exactly Spall’s (1992) motivation (multivariate finite differencing would cost 2*dim queries per step; SPSA costs 2 regardless of dimension).

Why it is in the benchmark. Per its own module docstring, spsa is “the only shipped estimator that never sees P” — every other static-T2 estimator (gls/spiess/vzw-entropy/od-congested) builds and consumes an explicit MSA route-choice proportion matrix; spsa treats the inner assignment purely as a black-box loss oracle, exactly how a microsimulator or a neural surrogate would be calibrated. It is also the methodological ancestor of spsa-sumo (batch-12), which swaps the inner oracle for a real microsimulator, and it is the ONLY static-T2 estimator registered as deterministic=False (seeded, not deterministic — macroreplicated like a stochastic model).

Scope. Seeded byte-reproducibility on the convex two-route corridor, the macroreplication-draws-distinct-perturbations regression (needs >=2 OD pairs — SPSA is invariant to the Rademacher SIGN in 1-D), and a certified Braess recovery.

Canon. [spall1992multivariate], docs/REFERENCES.md / docs/references.bib.

How this notebook is graded¶

A notebook never claims a number it does not compute in that cell. Every scored quantity below is recomputed live by the P1 ODCertifier from the emitted OD matrix against the harness’s own pinned BFW assignment of the truth — never from spsa’s self-report (README, Certified, not self-reported).

# Setup. `spsa` is a core estimator: a plain `pip install -e .` suffices —
# no optional extra, so no guard cell. The inline backend is Agg-based
# (headless CI renders into the notebook); NEVER matplotlib.use("Agg")
# in-kernel — it silently suppresses inline figure capture.
%matplotlib inline
import numpy as np

from tabench import (
    Budget,
    Demand,
    Network,
    ODCertifier,
    RngBundle,
    SPSAEstimator,
    Trace,
    braess_scenario,
    two_route_scenario,
    viz,
)
from tabench.core.rng import SOURCE_OBSERVATION
from tabench.estimation import ODTrace
from tabench.estimation.base import EstimationTask
from tabench.models.frank_wolfe import BiconjugateFrankWolfeModel
from tabench.observe.levels import LinkCounts

BRAESS_TRUTH = np.array([4.0, 2.0, 2.0, 2.0, 4.0])  # UE(D=6), recomputed below
TWOROUTE_TRUTH = np.array([2.5, 2.5, 1.5, 1.5])  # UE(D=4), recomputed below

# Analytic anchors recomputed, never quoted: pin both by an independent
# high-precision BFW solve before EITHER constant is used as a data-generating
# truth anywhere below.
def _pinned_ue(scenario):
    trace = Trace()
    BiconjugateFrankWolfeModel().solve(
        scenario, Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), trace
    )
    return trace.final.link_flows

np.testing.assert_allclose(_pinned_ue(braess_scenario(6.0)), BRAESS_TRUTH, atol=1e-6)
np.testing.assert_allclose(_pinned_ue(two_route_scenario(sue_theta=None)), TWOROUTE_TRUTH, atol=1e-6)


def _single_pair_prior(d):
    m = np.zeros((2, 2))
    m[0, 1] = d
    return m


def _task(scenario, truth, sensors, prior_matrix):
    ds = LinkCounts(np.asarray(sensors), 1, "none").observe(
        scenario, truth, RngBundle(0).generator(SOURCE_OBSERVATION)
    )
    return EstimationTask(
        name="t", network=scenario.network, prior=Demand(np.asarray(prior_matrix, dtype=np.float64)),
        dataset=ds, identifiability={}, scenario_hash=scenario.content_hash(), seed=0,
    )

SPSAEstimator.capabilities confirms the distinctive contract directly — compared against every other static-T2 estimator in this benchmark, spsa is the ONLY one registered as deterministic=False.

from tabench import GLSEstimator, SpiessEstimator, VZWEntropyEstimator
from tabench.estimation import Yang1992Estimator, DavisNihanKalmanEstimator

for cls in (GLSEstimator, SpiessEstimator, VZWEntropyEstimator, Yang1992Estimator,
            DavisNihanKalmanEstimator, SPSAEstimator):
    print(f"{cls.name:12s} deterministic={cls().capabilities.deterministic}")

caps = SPSAEstimator().capabilities
assert caps.deterministic is False  # stochastic (Rademacher perturbations), but seedable
assert caps.seedable is True
for cls in (GLSEstimator, SpiessEstimator, VZWEntropyEstimator, Yang1992Estimator,
            DavisNihanKalmanEstimator):
    assert cls().capabilities.deterministic is True

gls          deterministic=True
spiess       deterministic=True
vzw-entropy  deterministic=True
od-congested deterministic=True
od-kalman    deterministic=True
spsa         deterministic=False

Seeded smoke: byte-reproducible under a fixed macroreplication¶

Two-route, full sensors, a 200-sp-call budget: spsa recovers D=4 to |g-4|<0.2, and re-running with the SAME (seed, macrorep) reproduces the OD matrix BYTE-IDENTICALLY (P8 determinism contract).

sc2 = two_route_scenario(sue_theta=None)
task2 = _task(sc2, TWOROUTE_TRUTH, np.arange(4), _single_pair_prior(3.0))
budget2 = Budget(sp_calls=200, iterations=10**6)
t1 = ODTrace()
SPSAEstimator().estimate(task2, budget2, RngBundle(0, macrorep=0), t1)
print(f"recovered D : {t1.final.od_matrix[0, 1]:.4f}  (truth: 4.0)")
assert abs(t1.final.od_matrix[0, 1] - 4.0) < 0.2

t2 = ODTrace()
SPSAEstimator().estimate(task2, budget2, RngBundle(0, macrorep=0), t2)
print(f"byte-identical rerun : {np.array_equal(t1.final.od_matrix, t2.final.od_matrix)}")
assert np.array_equal(t1.final.od_matrix, t2.final.od_matrix)

recovered D : 4.0923  (truth: 4.0)
byte-identical rerun : True

Macroreplications draw distinct perturbations¶

SPSA’s Rademacher direction is invariant to sign in 1-D, so this needs a network with >=2 OD pairs: a 3-zone hub (zones 1,2,3 route through node 4). Two runs at macrorep=0 reproduce byte-identically; macrorep=1 draws a genuinely different perturbation sequence.

net = Network(
    name="hub2", n_nodes=4, n_zones=3, first_thru_node=4,
    init_node=np.array([1, 4, 4], dtype=np.int64),
    term_node=np.array([4, 2, 3], dtype=np.int64),
    capacity=np.array([4.0, 3.0, 3.0]), length=np.zeros(3),
    free_flow_time=np.array([1.0, 1.0, 1.0]), b=np.array([0.1, 0.1, 0.1]),
    power=np.ones(3), toll=np.zeros(3), link_type=np.ones(3, dtype=np.int64),
)
od = np.zeros((3, 3))
od[0, 1], od[0, 2] = 3.0, 2.0
from tabench import Scenario
sc_hub = Scenario("hub2", net, Demand(od), family="hub2")
truth_trace = Trace()
BiconjugateFrankWolfeModel().solve(
    sc_hub, Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), truth_trace
)
truth_hub = truth_trace.final.link_flows

prior = np.zeros((3, 3))
prior[0, 1], prior[0, 2] = 2.0, 3.0
ds = LinkCounts(np.array([1, 2]), 1, "none").observe(
    sc_hub, truth_hub, RngBundle(0).generator(SOURCE_OBSERVATION)
)
task_hub = EstimationTask("t", sc_hub.network, Demand(prior), ds, {}, sc_hub.content_hash(), seed=0)
budget_hub = Budget(sp_calls=400, iterations=10**6)

t0 = ODTrace()
SPSAEstimator().estimate(task_hub, budget_hub, RngBundle(0, macrorep=0), t0)
t0b = ODTrace()
SPSAEstimator().estimate(task_hub, budget_hub, RngBundle(0, macrorep=0), t0b)
t1_hub = ODTrace()
SPSAEstimator().estimate(task_hub, budget_hub, RngBundle(0, macrorep=1), t1_hub)

print(f"macrorep=0 repeat byte-identical : {np.array_equal(t0.final.od_matrix, t0b.final.od_matrix)}")
print(f"macrorep=0 vs macrorep=1 differ  : {not np.array_equal(t0.final.od_matrix, t1_hub.final.od_matrix)}")
assert np.array_equal(t0.final.od_matrix, t0b.final.od_matrix)
assert not np.array_equal(t0.final.od_matrix, t1_hub.final.od_matrix)

macrorep=0 repeat byte-identical : True
macrorep=0 vs macrorep=1 differ  : True

Certify: a Braess recovery, scored by the pinned harness¶

braess_scenario(6.0) is frozen and content-hashed (P2). A larger sp-call budget lets spsa recover the full 5-link Braess network from an off prior.

sc = braess_scenario(6.0)
print(f"scenario      : {sc.name}")
print(f"content hash  : {sc.content_hash()[:16]}…")

task_b = _task(sc, BRAESS_TRUTH, np.arange(5), _single_pair_prior(5.5))
trace_b = ODTrace()
SPSAEstimator(iters=200, k_inner=80).estimate(
    task_b, Budget(sp_calls=10**6, iterations=10**6), RngBundle(0, macrorep=0), trace_b
)
certifier = ODCertifier(
    sc, np.arange(5), np.array([], dtype=np.int64),
    BRAESS_TRUTH[None, :], BRAESS_TRUTH[[]][None, :], BRAESS_TRUTH,
    {"linear_identifiable": True},
)
metrics = certifier.certify(trace_b.final.od_matrix)
print(f"recovered D        : {trace_b.final.od_matrix[0, 1]:.4f}  (truth: 6.0)")
print(f"certified od_rmse  : {metrics['od_rmse']:.4f}")
assert metrics["od_feasible"] == 1.0
assert metrics["od_rmse"] < 1.0  # SPSA is stochastic, not asserted to gls/spiess-level tightness

scenario      : braess
content hash  : cf00f411cdccec88…
recovered D        : 6.3901  (truth: 6.0)
certified od_rmse  : 0.2758

Visualize¶

tabench.viz draws the road Network’s link flows for the certified spsa Braess recovery vs the truth, re-assigned from the OD matrix above.

from tabench import Scenario as _Scenario
recovered_trace = Trace()
BiconjugateFrankWolfeModel().solve(
    _Scenario("spsa-recovered", sc.network, Demand(trace_b.final.od_matrix), family=sc.family),
    Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), recovered_trace,
)
display(viz.plot_network_flows(sc.network, BRAESS_TRUTH))
display(viz.plot_flow_scatter(
    ("truth (D=6)", BRAESS_TRUTH),
    {"spsa recovered": recovered_trace.final.link_flows},
))

../../_images/0cdb6ca11f9d6af793a9d64f8c3ff388f7e19e40302cbab0dc5b56e1dfca2342.png

../../_images/28b4cd63e0fa2e1b83ae5e23ffc039784bdcc4de362bd1212e68f3e4a769f41a.png

Takeaways & pointers¶

Certified, not self-reported. The Braess od_rmse above came from ODCertifier’s own re-assignment of the emitted OD matrix, never from spsa’s self-report.
The only proportion-free estimator on the benchmark. Confirmed directly from SPSAEstimator.capabilities.inputs_required — SPSA works from assignment-oracle queries alone.
Seeded, not deterministic. deterministic=False but the RNG stream is pinned by (seed, macrorep) — verified as byte-identical reproducibility above, the P8 contract every stochastic estimator on this benchmark obeys.
Where next. the SUMO-in-the-loop extension spsa-sumo (batch-12, swaps the inner oracle for a real microsimulator); the proportion-based estimators gls / spiess; the lineage in the model compendium.

spsa — Spall’s (1992) simultaneous perturbation stochastic approximation¶