gls — Generalized Least Squares OD estimation (covers prior)

What. gls (T2, ADR-002) recovers an OD demand matrix from a set of observed link counts by minimizing a weighted sum of the prior misfit and the count misfit: min_g (g-g_pr)'W^-1(g-g_pr) + (p g - c)'V^-1(p g - c), where p is the (MSA-averaged) route-choice proportion matrix. It is the classical anchor every other T2 estimator in this benchmark is compared against. This notebook also covers prior — the stale-prior baseline gls is scored against is not a strawman but a first-class registered estimator (PriorBaseline, which emits the prior unchanged).

Why it is in the benchmark. OD estimation from traffic counts is underdetermined (Willumsen 1981, Cascetta 1984): gls is the field’s best-known regularized solution, and it needs its own MSA-proportion observation model (LinkCounts) and its own T2 held-out certification contract, both new relative to T1/day-to-day. See the model compendium (Cascetta 1984) and docs/design/adr-002-t2-estimation-certificate.md (P1).

Scope. This notebook certifies the scalar closed form, then the public run_estimation_experiment API end-to-end on Braess, comparing gls against the prior baseline it must beat, and finally the harder global-basin Braess D=6 recovery.

Canon. [cascetta1984estimation], 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 the estimator’s self-report (README, Certified, not self-reported).

# Setup. `gls` 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,
    GLSEstimator,
    ODCertifier,
    PriorBaseline,
    RngBundle,
    braess_scenario,
    run_estimation_experiment,
    two_route_scenario,
    viz,
)
from tabench.estimation import ODTrace, gls_solve
from tabench.estimation.base import EstimationTask
from tabench.core.rng import SOURCE_OBSERVATION
from tabench.observe.levels import LinkCounts
from tabench import Demand
from tabench.models.frank_wolfe import BiconjugateFrankWolfeModel
from tabench import Trace

The closed form: a scalar anchor

Single pair, single sensor, W=V=1: the GLS solve reduces to g* = (g_pr + p*c) / (1 + p^2) — a weighted average of the prior and the count-implied demand c/p. Recomputed here from gls_solve directly, no trusted digits.

p, c, g_pr = 0.625, 2.5, 3.0
expected = (g_pr + p * c) / (1.0 + p * p)
got = gls_solve(
    np.array([[p]]), np.array([c]), np.array([g_pr]), np.array([1.0]), np.array([1.0])
)
print(f"gls_solve      : {got[0]:.6f}")
print(f"closed form    : {expected:.6f}")
assert np.isclose(got[0], expected, atol=1e-10)

gls_solve      : 3.280899
closed form    : 3.280899

An end-to-end run: Braess, four sensors + one held out, gls vs prior

braess_scenario(6.0) is the D=6 Braess network (link order 1->3, 1->4, 3->4, 3->2, 4->2), frozen and content-hashed (P2). We run BOTH registered estimators through the public run_estimation_experiment API with a noiseless, exact-prior card: four links are observed and the middle link 3->4 is reserved as the held-out set — heldout_count_rmse is THE ranking column, so the design must reserve one (an empty held-out set is rejected on both T2 tracks, ADR-002/ADR-023). The prior baseline just echoes its (exact) input, so gls must do at least as well, and — under a noiseless card where the prior IS the truth — essentially exactly.

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

cfg = {
    "sensors": {"kind": "explicit", "links": [0, 1, 3, 4]},
    "heldout": {"kind": "explicit", "links": [2]},
    "n_periods": 1,
    "noise": "none",
    "prior": {"kind": "stale", "cv": 0.0},
    "identifiability_k_inner": 30,
}
result = run_estimation_experiment(
    sc, [PriorBaseline(), GLSEstimator(k_inner=40, outer_iters=8)],
    Budget(sp_calls=2000), seed=0, macroreps=1, estimation=cfg,
)
required = {
    "task_hash", "od_feasible", "obs_count_rmse", "oracle_obs_count_rmse",
    "heldout_count_rmse", "od_rmse", "od_nrmse", "total_demand_error",
    "od_identifiable",
}
assert required <= set(result.rows[0])
assert result.manifest["identifiability"]["linear_identifiable"] is True

prior_row = [r for r in result.rows if r["estimator"] == "prior"][-1]
gls_row = [r for r in result.rows if r["estimator"] == "gls"][-1]
print(f"prior  : od_feasible={prior_row['od_feasible']:.0f}  od_rmse={prior_row['od_rmse']:.4f}")
print(f"gls    : od_feasible={gls_row['od_feasible']:.0f}  od_rmse={gls_row['od_rmse']:.4f}")
assert prior_row["od_feasible"] == 1.0 and gls_row["od_feasible"] == 1.0
# cv=0 prior IS the truth: prior is exact, gls (an iterative fixed point
# around the same truth) stays within a small inner-assignment slack.
assert prior_row["od_rmse"] < 1e-9
assert gls_row["od_rmse"] < 1e-2

scenario      : braess
content hash  : cf00f411cdccec88…
prior  : od_feasible=1  od_rmse=0.0000
gls    : od_feasible=1  od_rmse=0.0007

The harder anchor: global-basin recovery from an off prior

From an off prior D=5.5 (inside gls’s global basin — see spiess for the D=4 SPURIOUS-basin caveat this network also has), full noiseless sensors let gls recover Braess’s true D=6 demand to a tight tolerance. Built by hand here (EstimationTask + LinkCounts, the same primitives run_estimation_experiment uses internally) so the certificate is scored against an explicit oracle.

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

trace_pin = Trace()
BiconjugateFrankWolfeModel().solve(
    sc, Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), trace_pin
)
oracle_flows = trace_pin.final.link_flows
np.testing.assert_allclose(oracle_flows, BRAESS_TRUTH, atol=1e-6)

prior_matrix = np.zeros((2, 2))
prior_matrix[0, 1] = 5.5
ds = LinkCounts(np.arange(5), 1, "none").observe(
    sc, oracle_flows, RngBundle(0).generator(SOURCE_OBSERVATION)
)
task = EstimationTask(
    name="t", network=sc.network, prior=Demand(prior_matrix), dataset=ds,
    identifiability={}, scenario_hash=sc.content_hash(), seed=0,
)
trace = ODTrace()
GLSEstimator(k_inner=120, outer_iters=80, cv_prior=50.0).estimate(
    task, Budget(sp_calls=10**9, iterations=200), RngBundle(0), trace
)
recovered = trace.final.od_matrix[0, 1]
print(f"recovered D    : {recovered:.4f}  (truth: 6.0, off prior: 5.5)")
assert abs(recovered - 6.0) < 1e-3

certifier = ODCertifier(
    sc, np.arange(5), np.array([], dtype=np.int64),
    oracle_flows[None, :], oracle_flows[[]][None, :], oracle_flows,
    {"linear_identifiable": True},
)
metrics = certifier.certify(trace.final.od_matrix)
print(f"certified od_rmse : {metrics['od_rmse']:.4e}")
assert metrics["od_feasible"] == 1.0
assert metrics["od_rmse"] < 1e-3

recovered D    : 6.0000  (truth: 6.0, off prior: 5.5)
certified od_rmse : 1.2068e-05

Visualize

tabench.viz draws the road Network’s link flows: the certified gls recovery vs the off prior, both re-derived from the OD matrices above by re-running the pinned assignment (the same one the certifier scores against).

prior_flows = trace_pin.final.link_flows  # truth flows for reference
prior_trace = Trace()
from tabench import Scenario
BiconjugateFrankWolfeModel().solve(
    Scenario("prior", sc.network, Demand(prior_matrix), family=sc.family),
    Budget(iterations=5000, target_relative_gap=1e-10), RngBundle(0), prior_trace,
)
display(viz.plot_network_flows(sc.network, trace_pin.final.link_flows))
display(viz.plot_flow_scatter(
    ("truth (D=6)", trace_pin.final.link_flows),
    {"off prior (D=5.5)": prior_trace.final.link_flows},
))

Takeaways & pointers

• Certified, not self-reported. Every od_rmse/od_feasible above came from ODCertifier’s own re-assignment of the emitted OD matrix, never from gls’s internal proportions.

• prior is a first-class competitor, not a strawman. It is the baseline every T2 estimator on this benchmark must beat; this notebook certifies it directly (ADR-002 covers: [prior]).

• Where next. the entropy-balancing alternative vzw-entropy; the spurious-local-minimum caveat on this SAME network in spiess; the spatial-covariance cousin od-kalman; the lineage in the model compendium.